【数字信号去噪】基于matlab小波软阈值+硬阈值+改进阈值轴承故障仿真信号去噪【含Matlab源码 1024期】

一、小波语音降噪简介

对于噪声频谱遍布于语音信号频谱之中的宽带噪声，如果噪声振幅比大部分的语音信号振幅低，则削去低幅度成分也就削去了宽带噪声。基于这种思路，可以在频域中采取中心限幅的方法，即让带噪语音信号通过一限幅滤波器，高幅度频谱可以通过而低幅成分不允许通过，从而实现噪声抑制。需要注意的是中心削波不可避免地要损害语音质量，通常只在频域中进行，而一般不在时域中实施。
小波降噪的原理类似于中心削波法。小波降噪最初是由Donoho和Johnstone提出的， 其主要理论依据是，小波变换具有很强的去数据相关性，它能够使信号的能量在小波域集中在一些大的小波系数中；而噪声的能量却分布于整个小波域内。因此，经小波分解后，信号的小波系数幅值要大于噪声的系数幅值。因此，幅值比较大的小波系数一般以信号为主，而幅值比较小的系数在很大程度上是噪声。于是，采用阈值的办法可以把信号系数保留，而使大部分噪声系数减小至0。小波降噪的具体处理过程为：将含噪信号在各尺度上进行小波分解，设定一个阈值，幅值低于该阈值的小波系数置为0，高于该阈值的小波系数或者完全保
留， 或者做相应的“收缩”(shrinkage) 处理。最后； 将处理后获得的小波系数用逆小波变换进行重构，得到去噪后的信号。
阈值去噪中，阈值函数体现了对超过和低于阈值的小波系数的不同处理策略，是阈值去噪中关键的一步。设w表示小波系数， T为给定阈值， sgn(*) 为符号函数， 常见的阈值函数主要有：

二、部分源代码

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)

  

 

  
1
  
2
  
3
  
4
  
5
  
6
  
7
  
8
  
9
  
10
  
11
  
12
  
13
  
14
  
15
  
16
  
17
  
18
  
19
  
20
  
21
  
22
  
23
  
24
  
25
  
26
  
27
  
28
  
29
  
30
  
31
  
32
  
33
  
34
  
35
  
36
  
37
  
38
  
39
  
40
  
41
  
42
  
43
  
44
  
45
  
46
  
47
  
48
  
49
  
50
  
51
  
52
  
53
  
54
  
55
  
56
  
57
  
58
  
59
  
60
  
61
  
62
  
63
  
64
  
65
  
66
  
67
  
68
  
69
  
70
  
71
  
72
  
73
  
74
  
75
  
76
  
77
  
78
  
79
  
80
  
81
  
82
  
83
  
84
  
85
  
86
  
87
  
88
  
89
  
90
  
91
  
92
  
93
  
94
  
95
  
96
  
97
  
98
  
99
  
100
  
101
  
102
  
103
  
104
  
105
  
106
  
107
  
108
  
109
  
110
  
111
  
112
  
113
  
114
  
115
  
116
  
117
  
118
  
119
  
120
  
121
  
122
  
123
  
124
  
125
  
126
  
127
  
128
  
129
  
130
  
131
  
132
  
133
 

三、运行结果

四、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 沈再阳.精通MATLAB信号处理[M].清华大学出版社，2015.
[2]高宝建,彭进业,王琳,潘建寿.信号与系统——使用MATLAB分析与实现[M].清华大学出版社，2020.
[3]王文光,魏少明,任欣.信号处理与系统分析的MATLAB实现[M].电子工业出版社，2018.

文章来源: qq912100926.blog.csdn.net，作者：海神之光，版权归原作者所有，如需转载，请联系作者。

原文链接：qq912100926.blog.csdn.net/article/details/118034096