# How effective is the signal denoising using the MATLAB based wavelet analysis

We see how to download seismic waveforms, convert them into mat format from mini-seed and then perform denoising using wavelet analysis. We first performed wavelet denoising on the high-frequency seismic time-series but after reconstruction. We also plot the wavelet coefficients at different scales. Further, I have shown on a synthetic time-series, how wavelet denoising works.

Wavelets analysis can be thought of as a general form of Fourier Analysis. Fourier Transform is often used in denoising the signals. Still, the biggest downside of this approach is that the signal needs to be stationary. Most of our real-world measurements are not stationary. Also, in Fourier-based denoising, we apply a lowpass filter to remove the noise. However, when the data has high-frequency features such as spikes in a signal or edges in an image, the lowpass filter smooths these out.

Here, the wavelet-based approach might have some advantages. Wavelets look at the signals in the multi-resolution window. It localizes features in the signal to a different scale. We can take advantage of that and preserve important signals, and removing nonuniform noise. I have covered the basics of the multi-resolution analysis using wavelets in previous post.

## Denoising using wavelets

When we take the wavelet transform of a time series, it concentrates the signal features in a few large-magnitude wavelet coefficients. In general, the coefficients that are smaller in value are considered noise. Theoretically, we can shrink those coefficients or simply remove them.

## Obtain the real-world signal

In this case, I will download the seismic time series recording a major earthquake. I will use the Obspy module to download the data from IRIS. For details on how to download waveforms using Obspy, see my previous post. I will download the waveforms for the arbitrarily selected event “2020-05-15 Mww 6.5 Nevada” located at (38.1689°N, 117.8497° W).

from obspy import read
from obspy.clients.fdsn import Client
from obspy import UTCDateTime
import matplotlib.pyplot as plt

client = Client("IRIS")

origin_time = UTCDateTime("2020-05-15T11:03:27Z")
net = "II"
stn = "PFO"
channel = "BHZ"
st = client.get_waveforms(net, stn, "00", channel,
origin_time+40, origin_time+300, attach_response=True)
st_copy = st.copy()
st_copy.remove_response(output="VEL")

# Filtering with a lowpass on a copy of the original Trace
st_copy.filter('highpass', freq=5.0, corners=2, zerophase=True)

st.plot(
outfile="myStream.png",
)

print(st_copy)
st.write(f'{net}-{stn}-{channel}.mseed', format="MSEED")


1 Trace(s) in Stream:
II.PFO.00.BHZ | 2020-05-15T11:04:07.019539Z - 2020-05-15T11:08:26.994539Z | 40.0 Hz, 10400 samples


This saves the mseed data into the script directory as well the plot of the “stream”. I also highpass filtered the signal to obtain the high-frequency part of the time series.

## Convert the mseed to mat format

Now, we want to make the data MATLAB readable. I have written a utility to convert any mseed format data to mat format. I will use that in this case.

python convert_mseed_mat.py -inp II-PFO-BHZ.mseed


## Denoise the signal using undecimated wavelet transform

I will use the Matlab function wdenoise to denoise the signal down to level 9 using the sym4 and db1 wavelets. wdenoise denoises the signal using an empirical Bayesian method with a Cauchy prior.

%% Using sym4
[XDEN,DENOISEDCFS] = wdenoise(data_double,9,'Wavelet','sym4','DenoisingMethod','BlockJS');
close;

fig2=figure('Renderer', 'painters', 'Position', [100 100 1000 400], 'color','w');
plot(data_double,'r')
hold on;
plot(XDEN, 'b')
legend('Original Signal','Denoised Signal')
axis tight;
hold off;
axis tight;
print(fig2,['denoise_signal_sym49','.jpg'],'-djpeg')

snrsym = -20*log10(norm(abs(data_double-XDEN))/norm(data_double))


In the above code, I denoised the seismic time-series using level 9 using block thresholding by setting the name-value pair, ‘DenoisingMethod’,’BlockJS’. I used the wavelet sym4.

The resultant difference in the signal strength compared with the original time-series is 82.5786.

Further, I try the same procedure using a different wavelet - db1.

%% Using db1
[XDEN_db1,DENOISEDCFS2] = wdenoise(data_double,9,'Wavelet','db1','DenoisingMethod','BlockJS');
close;

fig3=figure('Renderer', 'painters', 'Position', [100 100 1000 400], 'color','w');
plot(data_double,'r')
hold on;
plot(XDEN_db1, 'b')
legend('Original Signal','Denoised Signal')
axis tight;
hold off;
axis tight;
print(fig3,['denoise_signal_db1','.jpg'],'-djpeg')

snrdb1 = -20*log10(norm(abs(data_double-XDEN_db1))/norm(data_double))


This gives the resultant difference in the signal strength compared with the original time-series to be 34.5294. Hence, we find that the difference is significant in the wavelet, sym4. We can also play with several other parameters, including the levels.

I tried to change the level to 2, 5, 10 and got the SNR to be 82.7493, 82.5789, 82.5786, respectively, for the “sym4” wavelet. The insignificant difference between the SNR for different levels may show that the amount of high-frequency noise is not excessive in this case. Let us try to verify this hypothesis using the noisy synthetic data (in this case, we will have complete control).

%% Plot reconstructions based on the level-4 approximation
[sigden,coefs] = cmddenoise(data_double,'sym4',4);
[C,L] = wavedec(data_double,4,'sym4');

fig5= figure('Renderer', 'painters', 'Position', [100 100 1000 400], 'color','w');
app = wrcoef('a',C,L,'sym4',4); %Coefficients to reconstruct, specified as 'a' or 'd', for approximation or detail coefficients, respectively.
subplot(5,1,1);
plot(data_double, 'r-')
title('Approximation Coefficients');
hold on;
plot(app, 'b');
legend('Original Signal','Reconstructed Signal')
for nn = 1:4
det = wrcoef('d',C,L,'sym4',nn);
subplot(5,1,nn+1)
plot(det); title(['Noisy Wavelet Coefficients - Level '...
num2str(nn)]);
end
print(fig5,['wavelet_reconstruction','.jpg'],'-djpeg')


### Complete MATLAB code

clear; close; clc;
wdir='./';

fileloc0=[wdir,'II-PFO-BHZ'];
fileloc_ext = '.mat';
fileloc = [fileloc0 fileloc_ext];

if exist(fileloc,'file')
disp(['File exists ', fileloc]);

all_stats = fieldnames(stats);
all_data = fieldnames(data);

for id=1
stats_0 = stats.(all_stats{id});
data_0 = data.(all_data{id});

sampling_rate = stats_0.('sampling_rate');
delta = stats_0.('delta');
starttime = stats_0.('starttime');
endtime = stats_0.('endtime');
t1 = datetime(starttime,'InputFormat',"yyyy-MM-dd'T'HH:mm:ss.SSS'Z'");
t2 = datetime(endtime,'InputFormat',"yyyy-MM-dd'T'HH:mm:ss.SSS'Z'");
datetime_array = t1:seconds(delta):t2;

%% Plot waveforms
fig = figure('Renderer', 'painters', 'Position', [100 100 1000 400], 'color','w');
plot(t1:seconds(delta):t2, data_0, 'k-')
title([stats_0.('network'),'-', stats_0.('station'), '-', stats_0.('channel')])
axis tight;
print(fig,[fileloc0, '_ts', num2str(id),'.jpg'],'-djpeg')
data_double = double(data_0);

%% Using sym4
[XDEN,DENOISEDCFS] = wdenoise(data_double,9,'Wavelet','sym4','DenoisingMethod','BlockJS');
close;

fig2=figure('Renderer', 'painters', 'Position', [100 100 1000 400], 'color','w');
plot(data_double,'r')
hold on;
plot(XDEN, 'b')
legend('Original Signal','Denoised Signal')
axis tight;
hold off;
axis tight;
print(fig2,['denoise_signal_sym49','.jpg'],'-djpeg')

snrsym = -20*log10(norm(abs(data_double-XDEN))/norm(data_double))

%% Using db1
[XDEN_db1,DENOISEDCFS2] = wdenoise(data_double,9,'Wavelet','db1','DenoisingMethod','BlockJS');
close;

fig3=figure('Renderer', 'painters', 'Position', [100 100 1000 400], 'color','w');
plot(data_double,'r')
hold on;
plot(XDEN_db1, 'b')
legend('Original Signal','Denoised Signal')
axis tight;
hold off;
axis tight;
print(fig3,['denoise_signal_db1','.jpg'],'-djpeg')

snrdb1 = -20*log10(norm(abs(data_double-XDEN_db1))/norm(data_double))

end
end


## Wavelet denoising the noisy synthetic data

clear; close; clc;

rng default;
[X,XN] = wnoise('bumps',12,sqrt(6)); %returns values x of the test signal fun evaluated at 2n linearly spaced points from 0 to 1
denoised_signal = wdenoise(XN,4,'Wavelet','sym4','DenoisingMethod','BlockJS');
denoised_signal_db1 = wdenoise(XN,4,'Wavelet','db1','DenoisingMethod','BlockJS');

fig3=figure('Renderer', 'painters', 'Position', [100 100 1000 400], 'color','w');
ax1 = subplot(211)
plot(XN,'r')
hold on;
plot(denoised_signal,'b', 'linewidth',2)
legend('Noisy Data', 'Denoised Signal (sym4)')
hold off;
xlim([0, 4000])
axis tight;
ax2 = subplot(212)
plot(XN,'r')
hold on;
plot(denoised_signal_db1,'g', 'linewidth',2)
legend('Noisy Data', 'Denoised Signal (db1)')
hold off;
xlim([0, 4000])
axis tight;
print(fig3,['noisy_test_signal','.jpg'],'-djpeg')

snrsym = -20*log10(norm(abs(denoised_signal-XN))/norm(XN))
snrdb1 = -20*log10(norm(abs(denoised_signal_db1-XN))/norm(XN))


This example has been derived from the MATLAB example. In the above case, for the same data, I got the SNR to be 9.6994 and 9.4954 for sym4 and db1 wavelets. This also shows that there is a significant difference for the sym4. Please note that “sym4” works well in most cases and probably is the reason why it has been adopted as the default.