小波去噪的基本原理及其实现方法（Matlab）

weixin

小波去噪基本原理

小波去噪方法就是一种建立在小波变换多分辨分析基础上的算法，其基本思想是根据噪声与信号在不同频带上的小波分解系数具有不同强度分布的特点，将各频带上的噪声对应的小波系数去除，保留原始信号的小波分解系数，然后对处理后的系数进行小波重构，得到纯净信号。

相比以往的其他去噪方法，小波变换在低信噪比情况下的去噪效果较好，去噪后的信号识别率较高，同时小波去噪方法对时变信号和突变信号的去噪效果尤其明显。

图1 小波去噪基本原理图

Matlab实现小波去噪

利用小波分析对监测采集的信号进行去噪处理，可以恢复原始信号。利用小波分析进行去噪，包括以下3种方法：

  · 默认阈值去噪处理：该方法利用函数ddencmp() 生成信号的默认阈值，然后利用函数wdencmp() 进行去噪处理；

  · 给定阈值去噪处理：在实际的去噪处理过程中，阈值往往可通过经验公式获得，且这种阈值比默认阈值的可信度高。在进行阈值量化处理时可利用函数wthresh()；

  · 强制去噪处理：该方法是将小波分解结构中的高频系数全部置0，即滤掉所有高频部分，然后对信号进行小波重构。这种方法比较简单，且去噪后的信号比较平滑，但是容易丢失信号中的有用成分。

3种去噪方法原代码

[color=rgba(0, 0, 0, 0.5)]clear all;

[color=rgba(0, 0, 0, 0.5)]clc;

[color=rgba(0, 0, 0, 0.5)]load default.txt;

[color=rgba(0, 0, 0, 0.5)]%装载采集的信号

[color=rgba(0, 0, 0, 0.5)]x= default;

[color=rgba(0, 0, 0, 0.5)]lx=length(x);

[color=rgba(0, 0, 0, 0.5)]t=[0:1:length(x)-1]';

[color=rgba(0, 0, 0, 0.5)]%% 绘制监测所得信号%%

[color=rgba(0, 0, 0, 0.5)]subplot(2,2,1);

[color=rgba(0, 0, 0, 0.5)]plot(t,x);

[color=rgba(0, 0, 0, 0.5)]title('原始信号');

[color=rgba(0, 0, 0, 0.5)]grid on

[color=rgba(0, 0, 0, 0.5)]set(gcf,'color','w')

[color=rgba(0, 0, 0, 0.5)]set(gca,'fontname','times New Roman')

[color=rgba(0, 0, 0, 0.5)]set(gca,'fontsize',14.0)

[color=rgba(0, 0, 0, 0.5)]%% 用db1小波对原始信号进行3层分解并提取小波系数%%

[color=rgba(0, 0, 0, 0.5)][c,l]=wavedec(x,3,'db1');

[color=rgba(0, 0, 0, 0.5)]%sym8

[color=rgba(0, 0, 0, 0.5)]ca3=appcoef(c,l,'db1',3);

[color=rgba(0, 0, 0, 0.5)]%低频部分

[color=rgba(0, 0, 0, 0.5)]cd3=detcoef(c,l,3);

[color=rgba(0, 0, 0, 0.5)]%高频部分

[color=rgba(0, 0, 0, 0.5)]cd2=detcoef(c,l,2);

[color=rgba(0, 0, 0, 0.5)]%高频部分

[color=rgba(0, 0, 0, 0.5)]cd1=detcoef(c,l,1);

[color=rgba(0, 0, 0, 0.5)]%高频部分

[color=rgba(0, 0, 0, 0.5)]%% 对信号进行强制去噪处理并图示%%

[color=rgba(0, 0, 0, 0.5)]cdd3=zeros(1,length(cd3));

[color=rgba(0, 0, 0, 0.5)]cdd2=zeros(1,length(cd2));

[color=rgba(0, 0, 0, 0.5)]cdd1=zeros(1,length(cd1));

[color=rgba(0, 0, 0, 0.5)]c1=[ca3,cdd3,cdd2,cdd1];

[color=rgba(0, 0, 0, 0.5)]x1=waverec(c1,1,'db1');

[color=rgba(0, 0, 0, 0.5)]subplot(2,2,2);

[color=rgba(0, 0, 0, 0.5)]plot(x1);

[color=rgba(0, 0, 0, 0.5)]title('强制去噪后信号');

[color=rgba(0, 0, 0, 0.5)]grid on

[color=rgba(0, 0, 0, 0.5)]set(gcf,'color','w')

[color=rgba(0, 0, 0, 0.5)]set(gca,'fontname','times New Roman')

[color=rgba(0, 0, 0, 0.5)]set(gca,'fontsize',14.0)

[color=rgba(0, 0, 0, 0.5)]%% 默认阈值对信号去噪并图示%%

[color=rgba(0, 0, 0, 0.5)]%用ddencmp( )函数获得信号的默认阈值，使用wdencmp( )函数实现去噪过程

[color=rgba(0, 0, 0, 0.5)][thr,sorh,keepapp]=ddencmp('den','wv',x);

[color=rgba(0, 0, 0, 0.5)]x2=wdencmp('gbl',c,l,'db1',3,thr,sorh,keepapp);

[color=rgba(0, 0, 0, 0.5)]subplot(2,2,3);

[color=rgba(0, 0, 0, 0.5)]plot(x2);

[color=rgba(0, 0, 0, 0.5)]title('默认阈值去噪后信号');

[color=rgba(0, 0, 0, 0.5)]grid on

[color=rgba(0, 0, 0, 0.5)]set(gcf,'color','w')

[color=rgba(0, 0, 0, 0.5)]set(gca,'fontname','times New Roman')

[color=rgba(0, 0, 0, 0.5)]set(gca,'fontsize',14.0)

[color=rgba(0, 0, 0, 0.5)]%% 给定的软阈值进行去噪处理并图示%%

[color=rgba(0, 0, 0, 0.5)]cd1soft=wthresh(cd1,'x',1.465);

[color=rgba(0, 0, 0, 0.5)]%经验给出软阈值数

[color=rgba(0, 0, 0, 0.5)]cd2soft=wthresh(cd2,'x',1.823);

[color=rgba(0, 0, 0, 0.5)]%经验给出软阈值数

[color=rgba(0, 0, 0, 0.5)]cd3soft=wthresh(cd3,'x',2.768);

[color=rgba(0, 0, 0, 0.5)]%经验给出软阈值数

[color=rgba(0, 0, 0, 0.5)]c2=[ca3,cd3soft,cd2soft,cd1soft];

[color=rgba(0, 0, 0, 0.5)]x3=waverec(c2,1,'db1');

[color=rgba(0, 0, 0, 0.5)]subplot(2,2,4);

[color=rgba(0, 0, 0, 0.5)]plot(x3);

[color=rgba(0, 0, 0, 0.5)]title('给定软阈值去噪后信号');

[color=rgba(0, 0, 0, 0.5)]grid on

[color=rgba(0, 0, 0, 0.5)]set(gcf,'color','w')

[color=rgba(0, 0, 0, 0.5)]set(gca,'fontname','times New Roman')

[color=rgba(0, 0, 0, 0.5)]set(gca,'fontsize',14.0)

小波去噪：实例分享

[color=rgba(0, 0, 0, 0.5)]%在马达负载为1，转速为1771，内圈的故障信号

[color=rgba(0, 0, 0, 0.5)]load D:\我的文档\MATLAB\106.mat

[color=rgba(0, 0, 0, 0.5)]%读数据序列

[color=rgba(0, 0, 0, 0.5)]x=X106_DE_time(1:2000);

[color=rgba(0, 0, 0, 0.5)]%将信号1到2000个采样点给x

[color=rgba(0, 0, 0, 0.5)]%x=X106_DE_time;           %赋予变量

[color=rgba(0, 0, 0, 0.5)]fs=12000;                            %采样率

[color=rgba(0, 0, 0, 0.5)]dt=1/fs;                                %采样间隔时间1/fs

[color=rgba(0, 0, 0, 0.5)]n=length(x);                        %得到序列的长度

[color=rgba(0, 0, 0, 0.5)]t=[0:n-1]*dt;

[color=rgba(0, 0, 0, 0.5)]y=fft(x);                            %对信号进行快速FFT变换

[color=rgba(0, 0, 0, 0.5)]figure(1)

[color=rgba(0, 0, 0, 0.5)]subplot(2,2,1),

[color=rgba(0, 0, 0, 0.5)]plot(t,x);

[color=rgba(0, 0, 0, 0.5)]%绘制原始信号序列

[color=rgba(0, 0, 0, 0.5)]xlabel('时间/s'),

[color=rgba(0, 0, 0, 0.5)]title('时间域')

[color=rgba(0, 0, 0, 0.5)]grid on

[color=rgba(0, 0, 0, 0.5)]%subplot(2,2,2),

[color=rgba(0, 0, 0, 0.5)]%plot([0:n-1]/(n*dt),abs(y)*2/n)

[color=rgba(0, 0, 0, 0.5)]%绘制信号的振幅谱

[color=rgba(0, 0, 0, 0.5)]%xlabel('频率/HZ'),title('幅频图')

[color=rgba(0, 0, 0, 0.5)]%ylabel('振幅')

[color=rgba(0, 0, 0, 0.5)]%grid on

[color=rgba(0, 0, 0, 0.5)]%小波去噪1

[color=rgba(0, 0, 0, 0.5)]%将信号用小波函数SYM5分解到第5层，用minimaxi阈值选择对系数进行处理，消除噪声

[color=rgba(0, 0, 0, 0.5)]lev=5;

[color=rgba(0, 0, 0, 0.5)]xd=wden(x,'minimaxi','s','mln',lev,'sym5');

[color=rgba(0, 0, 0, 0.5)]subplot(2,2,2),

[color=rgba(0, 0, 0, 0.5)]plot(t,xd)

[color=rgba(0, 0, 0, 0.5)]xlabel('时间/s'),

[color=rgba(0, 0, 0, 0.5)]title('消噪后的信号')

[color=rgba(0, 0, 0, 0.5)]grid on

程序运行结果如下：

图2 去噪前

图3 去噪后

来源：本文内容整理自lxppopo新浪博客和gufengyueban的MATLAB论坛。