Digital signal denoising based on MATLAB wavelet soft threshold + hard threshold + improved threshold digital 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 y5 = awgn(y,0.5.'measured'); % Add white Gaussian noise
y10 = awgn(y,1.'measured'); % Add white Gaussian noise
y15 = awgn(y,5.'measured'); % Add White Gaussian noise %% SNR =0.5--------------------------------------------

mse12=MSE(y5,xdy12);
PSNR12=PSNR(y5,xdy12);
st=sprintf('Mean square error after denoising by soft threshold function =%.2f, SNR =%.2f',mse12,PSNR12); disp(st) st=sprintf('Mean square error after denoising by semi-soft threshold function =%.2f, SNR =%.2f',mse13,PSNR13); disp(st) figure(1);
subplot(511); plot(t,y); axis([0,inf,4 -.5]); title('Raw signal'); xlabel('time e (s)'); ylabel('幅度')
subplot(512); plot(t,y5); axis([0,inf,4 -.5]); title('Add 0.5 dB gaussian white noise signal'); xlabel('time e (s)'); ylabel('幅度')
subplot(513); plot(t,xdy11); axis([0,inf,4 -.5]); title('Signal processed by hard threshold function'); xlabel('time e (s)'); ylabel('幅度')
subplot(514); plot(t,xdy12); axis([0,inf,4 -.5]); title('Signal processed by soft threshold function'); xlabel('time e (s)'); ylabel('幅度')
subplot(515); plot(t,xdy13); axis([0,inf,4 -.5]); title('Signal processed by semi-soft threshold function'); xlabel('time e (s)'); ylabel('幅度') function X = denH (X, wNAME, n, THR) % hard threshold [C, S] = WaveDEC2 (X, n, wNAME); Dcoef = C(prod(S(1, :)) + 1: end); % ind = find(abs(dcoef) < thr) + prod(S(1:)); % is less than the threshold THR coefficient C(ind)=0; % is set to zero a=0.6;
ind = find( abs(dcoef) >= thr) + prod(S(1:)); % is greater than the threshold THR coefficient C(ind) = C(IND); Function X = denR (X, wname, n, THR) % soft threshold [C, S] = wavedec2(X, n, wNAME); Dcoef = C(prod(S(1, :)) + 1: end); % ind = find(abs(dcoef) < thr) + prod(S(1:)); % is less than the threshold THR coefficient C(ind) =0; % is set to zeroCopy the code

Version: 2014 a