【數字訊號去噪】基於matlab小波軟閾值+硬閾值+改進閾值軸承故障模擬訊號去噪【含Matlab原始碼 1024期】

阿新 • • 發佈：2021-06-19

一、簡介

基於matlab小波軟閾值+硬閾值+改進閾值軸承故障模擬訊號去噪

二、原始碼

clc
clear all
close all
fs = 20e3;                  % 取樣頻率
fn = 3e3;                   % 固有頻率
y0 = 5;                      % 位移常數
g = 0.1;                     % 阻尼係數
T = 0.01;                   % 重複週期
N = 4096;                  % 取樣點數
NT = round(fs*T);      % 單週期取樣點數
t = 0:1/fs:(N-1)/fs;      % 取樣時刻
t0 = 0:1/fs:(NT-1)/fs;  % 單週期取樣時刻
K = ceil(N/NT)+1;       % 重複次數
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);                % 頻譜
figure(1);subplot(231);
plot(t,y);
axis([0,inf,-4,5])
title('軸承故障模擬訊號時域波形圖')
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('sqtwolog-去噪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-去噪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-去噪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-去噪db5');xlabel('Time(s)');ylabel('Amplitude');
mse4=MSE(y5,xd54)
PSNR4=PSNR(y5,xd54)
[thr,sorh,keepapp]=ddencmp('den','wv',y5);%ddencmp用於獲取訊號在消噪或壓縮過程中的預設閾值
xd55=den_gaijin(y5,'db5',5,thr);
subplot(235);plot(t,xd55);axis([0,inf,-4,5]);title('改進的方法-去噪db5');xlabel('Time(s)');ylabel('Amplitude');
mse5=MSE(y5,xd55)
PSNR5=PSNR(y5,xd55)

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

xdy8=denh(y5,'db3',5,thr);%硬閾值
subplot(331);plot(t,xdy8);axis([0,inf,-4,5]);title('硬閾值db3-去噪');xlabel('Time(s)');ylabel('Amplitude');
mse81=MSE(y5,xdy8)
PSNR81=PSNR(y5,xdy8)
xdr8=dens(y5,'db3',5,thr);%軟閾值
subplot(332);plot(t,xdr8);axis([0,inf,-4,5]);title('軟閾值db3-去噪');xlabel('Time(s)');ylabel('Amplitude');
mse82=MSE(y5,xdr8)
PSNR82=PSNR(y5,xdr8)
xdA8=den_gaijin(y5,'db3',5,thr);%改進
subplot(333);plot(t,xdA8);axis([0,inf,-4,5]);title('改進db3-去噪');xlabel('Time(s)');ylabel('Amplitude');
mse83=MSE(y5,xdA8)
PSNR83=PSNR(y5,xdA8)
%%——————————————————————
xdy9=denh(y5,'db4',5,thr);%硬閾值
subplot(334);plot(t,xdy9);axis([0,inf,-4,5]);title('硬閾值db4-去噪');xlabel('Time(s)');ylabel('Amplitude');
mse91=MSE(y5,xdy9)
PSNR91=PSNR(y5,xdy9)
xdr9=dens(y5,'db4',5,thr);%軟閾值
subplot(335);plot(t,xdr9);axis([0,inf,-4,5]);title('軟閾值db4-去噪');xlabel('Time(s)');ylabel('Amplitude');
mse92=MSE(y5,xdr9)
PSNR92=PSNR(y5,xdr9)
xdA9=den_gaijin(y5,'db4',5,thr);%改進
subplot(336);plot(t,xdA9);axis([0,inf,-4,5]);title('改進db4-去噪');xlabel('Time(s)');ylabel('Amplitude');
mse93=MSE(y5,xdA9)
PSNR93=PSNR(y5,xdA9)
%%-----------------------------------
xdy10=denh(y5,'db5',5,thr);%硬閾值
subplot(337);plot(t,xdy10);axis([0,inf,-4,5]);title('硬閾值db5-去噪');xlabel('Time(s)');ylabel('Amplitude');
mse101=MSE(y5,xdy10)
PSNR101=PSNR(y5,xdy10)
xdr10=dens(y5,'db5',5,thr);%軟閾值
subplot(338);plot(t,xdr9);axis([0,inf,-4,5]);title('軟閾值db5-去噪');xlabel('Time(s)');ylabel('Amplitude');
mse102=MSE(y5,xdr10)
PSNR102=PSNR(y5,xdr10)
xdA10=den_gaijin(y5,'db5',5,thr);%改進
subplot(339);plot(t,xdA10);axis([0,inf,-4,5]);title('改進db5-去噪');xlabel('Time(s)');ylabel('Amplitude');
mse103=MSE(y5,xdA10)
PSNR103=PSNR(y5,xdA10)
%%--------------------------------

figure(4);

xdy8=denh(y5,'sym3',5,thr);%硬閾值
subplot(331);plot(t,xdy8);axis([0,inf,-4,5]);title('硬閾值sym3-去噪');xlabel('Time(s)');ylabel('Amplitude');
mse81=MSE(y5,xdy8)
PSNR81=PSNR(y5,xdy8)
xdr8=dens(y5,'sym3',5,thr);%軟閾值
subplot(332);plot(t,xdr8);axis([0,inf,-4,5]);title('軟閾值sym3-去噪');xlabel('Time(s)');ylabel('Amplitude');
mse82=MSE(y5,xdr8)
PSNR82=PSNR(y5,xdr8)
xdA8=den_gaijin(y5,'sym3',5,thr);%改進
subplot(333);plot(t,xdA8);axis([0,inf,-4,5]);title('改進sym3-去噪');xlabel('Time(s)');ylabel('Amplitude');
mse83=MSE(y5,xdA8)
PSNR83=PSNR(y5,xdA8)
%%——————————————————————
xdy9=denh(y5,'sym4',5,thr);%硬閾值
subplot(334);plot(t,xdy9);axis([0,inf,-4,5]);title('硬閾值sym4-去噪');xlabel('Time(s)');ylabel('Amplitude');
mse91=MSE(y5,xdy9)
PSNR91=PSNR(y5,xdy9)
xdr9=dens(y5,'sym4',5,thr);%軟閾值
subplot(335);plot(t,xdr9);axis([0,inf,-4,5]);title('軟閾值sym4-去噪');xlabel('Time(s)');ylabel('Amplitude');
mse92=MSE(y5,xdr9)
PSNR92=PSNR(y5,xdr9)
xdA9=den_gaijin(y5,'sym4',5,thr);%改進
subplot(336);plot(t,xdA9);axis([0,inf,-4,5]);title('改進sym4-去噪');xlabel('Time(s)');ylabel('Amplitude');
mse93=MSE(y5,xdA9)
PSNR93=PSNR(y5,xdA9)
%%-----------------------------------
xdy10=denh(y5,'sym5',5,thr);%硬閾值
subplot(337);plot(t,xdy10);axis([0,inf,-4,5]);title('硬閾值sym5-去噪');xlabel('Time(s)');ylabel('Amplitude');
mse101=MSE(y5,xdy10)
PSNR101=PSNR(y5,xdy10)
xdr10=dens(y5,'sym5',5,thr);%軟閾值
subplot(338);plot(t,xdr9);axis([0,inf,-4,5]);title('軟閾值sym5-去噪');xlabel('Time(s)');ylabel('Amplitude');
mse102=MSE(y5,xdr10)
PSNR102=PSNR(y5,xdr10)
xdA10=den_gaijin(y5,'sym5',5,thr);%改進
subplot(339);plot(t,xdA10);axis([0,inf,-4,5]);title('改進sym5-去噪');xlabel('Time(s)');ylabel('Amplitude');
mse103=MSE(y5,xdA10)