Digital signal denoising based on MATLAB wavelet soft threshold + hard threshold + improved threshold bearing fault simulation signal denoising

clc
clear all
close all
fs = 20e3; % Sampling frequency fn =3e3; % natural frequency y0 =5; % displacement constant g =0.1; % damping coefficient T =0.01; % repeat period N =4096; NT = round(fs*T); % Number of sampling points in one period t =0:1/fs:(N- 1)/fs; % sampling time t0 =0:1/fs:(NT- 1)/fs; % Single-period sampling time K =ceil(N/NT)+1; % number of repetitions y = [];for i = 1:K
    y = [y,y0*exp(-g*2*pi*fn*t0).*sin(2*pi*fn*sqrt(1-g^2)*t0)];
end
y = y(1:N); Yf = fft(y); % spectrum figure (1); subplot(231);
plot(t,y);
axis([0,inf,4 -.5])
title('Bearing Fault Simulation Signal Time domain Waveform')
xlabel('Time(s)')
ylabel('Amplitude')
y5 = awgn(y,5.'measured'); % Add white Gaussian noise
y10 = awgn(y,10.'measured'); % Add white Gaussian noise
y15 = awgn(y,15.'measured'); % Add white Gaussian noise
figure(2);
xd51 = wden(y5,'sqtwolog'.'s'.'one'.5.'db5');
subplot(231); plot(t,xd51); axis([0,inf,2 -.2]); title('the sqtwolog - denoising db5'); xlabel('Time(s)'); ylabel('Amplitude');
mse1=MSE(y5,xd51)
PSNR1=PSNR(y5,xd51)
xd52 = wden(y5,'rigrsure'.'s'.'one'.5.'db5');
subplot(232); plot(t,xd52); axis([0,inf,- 3.4]); title('rigrsure - denoising db5'); xlabel('Time(s)'); ylabel('Amplitude');
mse2=MSE(y5,xd52)
PSNR2=PSNR(y5,xd52)
xd53 = wden(y5,'heursure'.'s'.'one'.5.'db5');
subplot(233); plot(t,xd53); axis([0,inf,2 -.2]); title('heursure - denoising db5'); xlabel('Time(s)'); ylabel('Amplitude');
mse3=MSE(y5,xd53)
PSNR3=PSNR(y5,xd53)
xd54 = wden(y5,'minimaxi'.'s'.'one'.5.'db5');
subplot(234); plot(t,xd54); axis([0,inf,2 -.2]); title('minimaxi - denoising db5'); xlabel('Time(s)'); ylabel('Amplitude');
mse4=MSE(y5,xd54)
PSNR4=PSNR(y5,xd54)
[thr,sorh,keepapp]=ddencmp('den'.'wv',y5); % ddencMP is used to obtain the default threshold of signal during denoising or compression xd55=den_gaijin(y5,'db5'.5,thr);
subplot(235); plot(t,xd55); axis([0,inf,4 -.5]); title('Improved method - denoising DB5'); xlabel('Time(s)'); ylabel('Amplitude');
mse5=MSE(y5,xd55)
PSNR5=PSNR(y5,xd55)

%%--------------------------------------------
figure(3);

xdy8=denh(y5,'db3'.5,thr); % hard threshold subplot(331); plot(t,xdy8); axis([0,inf,4 -.5]); title('Hard threshold DB3 - Denoising'); xlabel('Time(s)'); ylabel('Amplitude');
mse81=MSE(y5,xdy8)
PSNR81=PSNR(y5,xdy8)
xdr8=dens(y5,'db3'.5,thr); % soft threshold subplot(332); plot(t,xdr8); axis([0,inf,4 -.5]); title('Soft threshold DB3 - Denoising'); xlabel('Time(s)'); ylabel('Amplitude');
mse82=MSE(y5,xdr8)
PSNR82=PSNR(y5,xdr8)
xdA8=den_gaijin(y5,'db3'.5,thr); % improvement subplot (333); plot(t,xdA8); axis([0,inf,4 -.5]); title('Improved DB3 - Denoising'); xlabel('Time(s)'); ylabel('Amplitude'); Mse83 = MSE (y5, xdA8) PSNR83 = PSNR (y5, xdA8) % % -- -- -- -- -- -- -- -- -- -- -- xdy9 = denh (y5,'db4'.5,thr); % hard threshold subplot(334); plot(t,xdy9); axis([0,inf,4 -.5]); title('Hard threshold DB4 - Denoising'); xlabel('Time(s)'); ylabel('Amplitude');
mse91=MSE(y5,xdy9)
PSNR91=PSNR(y5,xdy9)
xdr9=dens(y5,'db4'.5,thr); % soft threshold subplot(335); plot(t,xdr9); axis([0,inf,4 -.5]); title('Soft threshold DB4 - Denoising'); xlabel('Time(s)'); ylabel('Amplitude');
mse92=MSE(y5,xdr9)
PSNR92=PSNR(y5,xdr9)
xdA9=den_gaijin(y5,'db4'.5,thr); % improvement subplot (336); plot(t,xdA9); axis([0,inf,4 -.5]); title('Improved DB4 - Denoising'); xlabel('Time(s)'); ylabel('Amplitude');
mse93=MSE(y5,xdA9)
PSNR93=PSNR(y5,xdA9)
%%-----------------------------------
xdy10=denh(y5,'db5'.5,thr); % hard threshold subplot(337); plot(t,xdy10); axis([0,inf,4 -.5]); title('Hard threshold DB5 - Denoising'); xlabel('Time(s)'); ylabel('Amplitude');
mse101=MSE(y5,xdy10)
PSNR101=PSNR(y5,xdy10)
xdr10=dens(y5,'db5'.5,thr); % soft threshold subplot(338); plot(t,xdr9); axis([0,inf,4 -.5]); title('Soft threshold DB5 - Denoising'); xlabel('Time(s)'); ylabel('Amplitude');
mse102=MSE(y5,xdr10)
PSNR102=PSNR(y5,xdr10)
xdA10=den_gaijin(y5,'db5'.5,thr); % improvement subplot (339); plot(t,xdA10); axis([0,inf,4 -.5]); title('Improved DB5 - Denoising'); xlabel('Time(s)'); ylabel('Amplitude');
mse103=MSE(y5,xdA10)
PSNR103=PSNR(y5,xdA10)
%%--------------------------------

figure(4);

xdy8=denh(y5,'sym3'.5,thr); % hard threshold subplot(331); plot(t,xdy8); axis([0,inf,4 -.5]); title('Hard threshold SYM3 - Denoising'); xlabel('Time(s)'); ylabel('Amplitude');
mse81=MSE(y5,xdy8)
PSNR81=PSNR(y5,xdy8)
xdr8=dens(y5,'sym3'.5,thr); % soft threshold subplot(332); plot(t,xdr8); axis([0,inf,4 -.5]); title('Soft Threshold SYM3 - Denoising'); xlabel('Time(s)'); ylabel('Amplitude');
mse82=MSE(y5,xdr8)
PSNR82=PSNR(y5,xdr8)
xdA8=den_gaijin(y5,'sym3'.5,thr); % improvement subplot (333); plot(t,xdA8); axis([0,inf,4 -.5]); title('Improved SYM3 - Denoising'); xlabel('Time(s)'); ylabel('Amplitude'); Mse83 = MSE (y5, xdA8) PSNR83 = PSNR (y5, xdA8) % % -- -- -- -- -- -- -- -- -- -- -- xdy9 = denh (y5,'sym4'.5,thr); % hard threshold subplot(334); plot(t,xdy9); axis([0,inf,4 -.5]); title('Hard threshold SYM4 - Denoising'); xlabel('Time(s)'); ylabel('Amplitude');
mse91=MSE(y5,xdy9)
PSNR91=PSNR(y5,xdy9)
xdr9=dens(y5,'sym4'.5,thr); % soft threshold subplot(335); plot(t,xdr9); axis([0,inf,4 -.5]); title('Soft Threshold SYM4 - Denoising'); xlabel('Time(s)'); ylabel('Amplitude');
mse92=MSE(y5,xdr9)
PSNR92=PSNR(y5,xdr9)
xdA9=den_gaijin(y5,'sym4'.5,thr); % improvement subplot (336); plot(t,xdA9); axis([0,inf,4 -.5]); title('Improved SYM4 - Denoising'); xlabel('Time(s)'); ylabel('Amplitude');
mse93=MSE(y5,xdA9)
PSNR93=PSNR(y5,xdA9)
%%-----------------------------------
xdy10=denh(y5,'sym5'.5,thr); % hard threshold subplot(337); plot(t,xdy10); axis([0,inf,4 -.5]); title('Hard threshold SYM5 - Denoising'); xlabel('Time(s)'); ylabel('Amplitude');
mse101=MSE(y5,xdy10)
PSNR101=PSNR(y5,xdy10)
xdr10=dens(y5,'sym5'.5,thr); % soft threshold subplot(338); plot(t,xdr9); axis([0,inf,4 -.5]); title('Soft Threshold SYM5 - Denoising'); xlabel('Time(s)'); ylabel('Amplitude');
mse102=MSE(y5,xdr10)
PSNR102=PSNR(y5,xdr10)
xdA10=den_gaijin(y5,'sym5'.5,thr); % improvement subplot (339); plot(t,xdA10); axis([0,inf,4 -.5]); title('Improved SYM5 - Denoising'); xlabel('Time(s)'); ylabel('Amplitude');
mse103=MSE(y5,xdA10)
Copy the code

Version: 2014 a