Subscribe to our weekly newsletter

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

Utpal Kumar   6 minute read      

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.

Similar posts

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[0].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.

Vertical component waveforms at II-PFO for 2020-05-15 Mww 6.5 Nevada
Vertical component waveforms at II-PFO for 2020-05-15 Mww 6.5 Nevada

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
Vertical component waveforms at II-PFO for 2020-05-15 Mww 6.5 Nevada using MATLAB
Vertical component waveforms at II-PFO for 2020-05-15 Mww 6.5 Nevada using MATLAB

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))
Denoised seismic time-series using sym4
Denoised seismic time-series using sym4

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))
Denoised seismic time-series using db1
Denoised seismic time-series using db1

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')
Wavelet reconstructions based on the level-4 approximation
Wavelet reconstructions based on the level-4 approximation

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]);
    load(fileloc);
    
    all_stats = fieldnames(stats);
    all_data = fieldnames(data);
    
    for id=1
        %% read data and header information
        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;
linkaxes([ax1 ax2],'x')
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.

Wavelet denoising for a synthetic data generated by MATLAB
Wavelet denoising for a synthetic data generated by MATLAB

Conclusions

We have seen how to download seismic waveforms, convert them into mat format from mseed and then perform denoising using wavelet analysis. We first performed wavelet denoising on the high-frequency seismic time series, but after reconstruction, the difference is not apparent visually. But from the SNR, we can see the differences. We have also plotted the wavelet coefficients at different scales. Further, I have shown on a synthetic time-series how wavelet denoising works.

References

  1. Wu, W., & Cai, P. (n.d.). Simulation of Wavelet Denoising Based on MATLAB–《Information and Electronic Engineering》2008年03期. Retrieved May 18, 2021, from https://en.cnki.com.cn/Article_en/CJFDTotal-XXYD200803017.htm
  2. MATLAB

Disclaimer of liability

The information provided by the Earth Inversion is made available for educational purposes only.

Whilst we endeavor to keep the information up-to-date and correct. Earth Inversion makes no representations or warranties of any kind, express or implied about the completeness, accuracy, reliability, suitability or availability with respect to the website or the information, products, services or related graphics content on the website for any purpose.

UNDER NO CIRCUMSTANCE SHALL WE HAVE ANY LIABILITY TO YOU FOR ANY LOSS OR DAMAGE OF ANY KIND INCURRED AS A RESULT OF THE USE OF THE SITE OR RELIANCE ON ANY INFORMATION PROVIDED ON THE SITE. ANY RELIANCE YOU PLACED ON SUCH MATERIAL IS THEREFORE STRICTLY AT YOUR OWN RISK.

Leave a comment