七种滤波方法测试matlab实现（转）

本文源自：http://blog.163.com/xiaheng0804@126/blog/static/1205282120132129471816/
创建两个混合信号，便于更好测试滤波器效果。同时用七中滤波方法测试。

混合信号ix_Signal_1 = 信号Signal_Original_1+白噪声。

混合信号Mix_Signal_2 = 信号Signal_Original_2+白噪声。

1.巴特沃斯低通滤波器去噪
巴特沃斯滤波器适合用于信号和噪声没有重叠的情况下。下图是巴特沃斯对两个信号的滤波效果。

从图上可以看出巴特沃斯低通滤波器对信号一的滤波效果还是可以的，主要是因为有效的信号最高频率才30Hz,本程序将50Hz以上的信号全部滤除，通过的频率成分中仍然是有白噪声的。
对于信号二，滤波后的信号与没有加噪声的信号相比就有失真了，上升沿和下降沿的高频信号被滤除了。
2.FIR低通滤波器去噪
情况同巴特沃斯低通滤波器相似。滤波后的效果如下：

3. 移动平均滤波去噪
滤波效果如下：

4. 中值滤波去噪

从上图可以看出，无论是对信号一还是对信号二，中值滤波的滤波效果都是很不错，特备是对于信号二，上升沿和下降失真比较的小。5. 维纳滤波去噪

维纳滤波器属于现代滤波器，传统的滤波器只能滤除信号和干扰频带没有重叠的情况，当信号和干扰频带有重叠的时候传统滤波器将无能为力，这时就需要用到现代滤波器，现代滤波器利用信号和干扰的统计特征（如自相关函数、功率谱等）导出一套最佳估值算法，然后用硬件或软件予以实现。
维纳滤波是以均方误差最小（LMS（LeastMeanSquare）为准则的，它根据过去观测值和当前观测值来估计信号的当前值，因此它的解形式是系统的传递函数或单位脉冲响应。
均方误差为：

维纳－霍夫（Wiener－Hopf）方程最小均方误差下的解为：

由于理解不深，对于信号二，没有什么滤波效果

6. 自适应滤波去噪
维纳滤波器参数是固定的，适合于平稳随机信号。卡尔曼滤波器参数是时变的，适合于非平稳随机信号。然而，只有在信号和噪声的统计特性先验已知的情况下，这两种滤波技术才能获得最优滤波。在实际应用中，常常无法得到信号和噪声统计特性的先验知识。在这种情况下，自适应滤波技术能够获得极佳的滤波性能，因而具有很好的应用价值。
自适应滤波的滤波效果如下：

本程序是基于LMS算法的自适应滤波，从上图可以看出，滤波效果也是很不错的，特别是对于信号二，上升沿有失真，下降沿保持还可以，最要的是得到的波形十分的平滑。由此可见自适应滤波极具使用价值。7. 小波去噪

首先看一下小波的去噪效果。

对于信号二，小波的去噪效果非常不错，虽然得到波形不是很平滑，但是上升沿和下降沿保持的非常高，基本可以看到棱角.
%****************************************************************************************
%  
%                      创建两个信号Mix_Signal_1和信号 Mix_Signal_2 
%
%***************************************************************************************

Fs =1000;                                                                        %采样率
N  =1000;                                                                        %采样点数
n  = 0:N-1;
t   =0:1/Fs:1-1/Fs;                                                            %时间序列 
Signal_Original_1=sin(2*pi*10*t)+sin(2*pi*20*t)+sin(2*pi*30*t); 
Noise_White_1    =[0.3*randn(1,500),rand(1,500)];          %前500点高斯分部白噪声，后500点均匀分布白噪声
Mix_Signal_1   =Signal_Original_1 +Noise_White_1;        %构造的混合信号

Signal_Original_2  =  [zeros(1,100),20*ones(1,20), -2*ones(1,30), 5*ones(1,80), -5*ones(1,30),9*ones(1,140), -4*ones(1,40), 3*ones(1,220), 12*ones(1,100),5*ones(1,20), 25*ones(1,30), 7*ones(1,190)]; 
Noise_White_2    =  0.5*randn(1,1000);                                %高斯白噪声
Mix_Signal_2        =  Signal_Original_2+Noise_White_2;      %构造的混合信号

%****************************************************************************************
%  
%                信号Mix_Signal_1和Mix_Signal_2  分别作巴特沃斯低通滤波。
%
%***************************************************************************************

%混合信号Mix_Signal_1  巴特沃斯低通滤波
figure(1);
Wc=2*50/Fs;                                          %截止频率50Hz
[b,a]=butter(4,Wc);
Signal_Filter=filter(b,a,Mix_Signal_1);

subplot(4,1,1);                                        %Mix_Signal_1原始信号                 
plot(Mix_Signal_1);
axis([0,1000,-4,4]);
title('原始信号 ');

subplot(4,1,2);                                        %Mix_Signal_1低通滤波滤波后信号  
plot(Signal_Filter);
axis([0,1000,-4,4]);
title('巴特沃斯低通滤波后信号');

%混合信号Mix_Signal_2  巴特沃斯低通滤波
Wc=2*100/Fs;                                          %截止频率100Hz
[b,a]=butter(4,Wc);
Signal_Filter=filter(b,a,Mix_Signal_2);

subplot(4,1,3);                                        %Mix_Signal_2原始信号                 
plot(Mix_Signal_2);
axis([0,1000,-10,30]);
title('原始信号 ');

subplot(4,1,4);                                      %Mix_Signal_2低通滤波滤波后信号  
plot(Signal_Filter);
axis([0,1000,-10,30]);
title('巴特沃斯低通滤波后信号');

%****************************************************************************************
%  
%                信号Mix_Signal_1和Mix_Signal_2  分别作FIR低通滤波。
%
%***************************************************************************************

%混合信号Mix_Signal_1  FIR低通滤波
figure(2);
F  =  [0:0.05:0.95]; 
A  =  [1    1      0    0    0    0      0    0    0    0    0    0    0    0    0    0    0  0  0   0] ;
b  =  firls(20,F,A);
Signal_Filter = filter(b,1,Mix_Signal_1);

subplot(4,1,1);                                          %Mix_Signal_1原始信号                 
plot(Mix_Signal_1);
axis([0,1000,-4,4]);
title('原始信号 ');

subplot(4,1,2);                                          %Mix_Signal_1FIR低通滤波滤波后信号  
plot(Signal_Filter);
axis([0,1000,-5,5]);
title('FIR低通滤波后的信号');

%混合信号Mix_Signal_2  FIR低通滤波
F  =  [0:0.05:0.95]; 
A  =  [1    1      1    1    1    0      0    0    0    0    0    0    0    0    0    0    0  0  0   0] ;
b  =  firls(20,F,A);
Signal_Filter = filter(b,1,Mix_Signal_2);
subplot(4,1,3);                                          %Mix_Signal_2原始信号                 
plot(Mix_Signal_2);
axis([0,1000,-10,30]);
title('原始信号 ');

subplot(4,1,4);                                          %Mix_Signal_2FIR低通滤波滤波后信号  
plot(Signal_Filter);
axis([0,1000,-10,30]);
title('FIR低通滤波后的信号');

%****************************************************************************************
%  
%                信号Mix_Signal_1和Mix_Signal_2  分别作移动平均滤波
%
%***************************************************************************************

%混合信号Mix_Signal_1  移动平均滤波
figure(3);
b  =  [11 1 1 1 1]/6;
Signal_Filter = filter(b,1,Mix_Signal_1);

subplot(4,1,1);                                          %Mix_Signal_1原始信号                 
plot(Mix_Signal_1);
axis([0,1000,-4,4]);
title('原始信号 ');

subplot(4,1,2);                                          %Mix_Signal_1移动平均滤波后信号  
plot(Signal_Filter);
axis([0,1000,-4,4]);
title('移动平均滤波后的信号');

%混合信号Mix_Signal_2  移动平均滤波
b  =  [11 1 1 1 1]/6;
Signal_Filter = filter(b,1,Mix_Signal_2);
subplot(4,1,3);                                          %Mix_Signal_2原始信号                 
plot(Mix_Signal_2);
axis([0,1000,-10,30]);
title('原始信号 ');

subplot(4,1,4);                                          %Mix_Signal_2移动平均滤波后信号  
plot(Signal_Filter);
axis([0,1000,-10,30]);
title('移动平均滤波后的信号');

%****************************************************************************************
%  
%                信号Mix_Signal_1和 Mix_Signal_2  分别作中值滤波
%
%***************************************************************************************

%混合信号Mix_Signal_1  中值滤波
figure(4);
Signal_Filter=medfilt1(Mix_Signal_1,10);

subplot(4,1,1);                                          %Mix_Signal_1原始信号                 
plot(Mix_Signal_1);
axis([0,1000,-5,5]);
title('原始信号 ');

subplot(4,1,2);                                          %Mix_Signal_1中值滤波后信号  
plot(Signal_Filter);
axis([0,1000,-5,5]);
title('中值滤波后的信号');

%混合信号Mix_Signal_2  中值滤波
Signal_Filter=medfilt1(Mix_Signal_2,10);
subplot(4,1,3);                                          %Mix_Signal_2原始信号                 
plot(Mix_Signal_2);
axis([0,1000,-10,30]);
title('原始信号 ');

subplot(4,1,4);                                          %Mix_Signal_2中值滤波后信号  
plot(Signal_Filter);
axis([0,1000,-10,30]);
title('中值滤波后的信号');

%****************************************************************************************
%  
%                信号Mix_Signal_1和 Mix_Signal_2  分别作维纳滤波
%
%***************************************************************************************

%混合信号Mix_Signal_1  维纳滤波
figure(5);
Rxx=xcorr(Mix_Signal_1,Mix_Signal_1);              %得到混合信号的自相关函数
M=100;                                                            %维纳滤波器阶数
fori=1:M                                                          %得到混合信号的自相关矩阵
    forj=1:M
        rxx(i,j)=Rxx(abs(j-i)+N);
    end
end
Rxy=xcorr(Mix_Signal_1,Signal_Original_1);      %得到混合信号和原信号的互相关函数
for i=1:M
    rxy(i)=Rxy(i+N-1);
end                                                                  %得到混合信号和原信号的互相关向量
h =inv(rxx)*rxy';                                              %得到所要涉及的wiener滤波器系数
Signal_Filter=filter(h,1,Mix_Signal_1);              %将输入信号通过维纳滤波器

subplot(4,1,1);                                                  %Mix_Signal_1原始信号                 
plot(Mix_Signal_1);
axis([0,1000,-5,5]);
title('原始信号 ');

subplot(4,1,2);                                                  %Mix_Signal_1 维纳滤波后信号  
plot(Signal_Filter);
axis([0,1000,-5,5]);
title('维纳滤波后的信号');

%混合信号Mix_Signal_2  维纳滤波
Rxx=xcorr(Mix_Signal_2,Mix_Signal_2);              %得到混合信号的自相关函数
M=500;                                                            %维纳滤波器阶数
fori=1:M                                                          %得到混合信号的自相关矩阵
    forj=1:M
        rxx(i,j)=Rxx(abs(j-i)+N);
    end
end
Rxy=xcorr(Mix_Signal_2,Signal_Original_2);      %得到混合信号和原信号的互相关函数
for i=1:M
    rxy(i)=Rxy(i+N-1);
end                                                                  %得到混合信号和原信号的互相关向量
h=inv(rxx)*rxy';                                              %得到所要涉及的wiener滤波器系数
Signal_Filter=filter(h,1,Mix_Signal_2);            %将输入信号通过维纳滤波器

subplot(4,1,3);                                                  %Mix_Signal_2原始信号                 
plot(Mix_Signal_2);
axis([0,1000,-10,30]);
title('原始信号 ');

subplot(4,1,4);                                                  %Mix_Signal_2维纳滤波后信号  
plot(Signal_Filter);
axis([0,1000,-10,30]);
title('维纳滤波后的信号');

%****************************************************************************************
%  
%                信号Mix_Signal_1和Mix_Signal_2  分别作自适应滤波
%
%***************************************************************************************

%混合信号 Mix_Signal_1 自适应滤波
figure(6);
N=1000;                                            %输入信号抽样点数N
k=100;                                                  %时域抽头LMS算法滤波器阶数
u=0.001;                                            %步长因子

%设置初值
yn_1=zeros(1,N);                                  %outputsignal
yn_1(1:k)=Mix_Signal_1(1:k);                %将输入信号SignalAddNoise的前k个值作为输出yn_1的前k个值
w=zeros(1,k);                                        %设置抽头加权初值
e=zeros(1,N);                                        %误差信号

%用LMS算法迭代滤波
for i=(k+1):N
        XN=Mix_Signal_1((i-k+1):(i));
        yn_1(i)=w*XN';
        e(i)=Signal_Original_1(i)-yn_1(i);
        w=w+2*u*e(i)*XN;
end

subplot(4,1,1);
plot(Mix_Signal_1);                              %Mix_Signal_1 原始信号
axis([k+1,1000,-4,4]);
title('原始信号');

subplot(4,1,2);
plot(yn_1);                                            %Mix_Signal_1自适应滤波后信号
axis([k+1,1000,-4,4]);
title('自适应滤波后信号');

%混合信号 Mix_Signal_2 自适应滤波
N=1000;                                            %输入信号抽样点数N
k=500;                                                %时域抽头LMS算法滤波器阶数
u=0.000011;                                        %步长因子

%设置初值
yn_1=zeros(1,N);                                  %output signal
yn_1(1:k)=Mix_Signal_2(1:k);                  %将输入信号SignalAddNoise的前k个值作为输出yn_1的前k个值
w=zeros(1,k);                                        %设置抽头加权初值
e=zeros(1,N);                                        %误差信号

%用LMS算法迭代滤波
for i=(k+1):N
        XN=Mix_Signal_2((i-k+1):(i));
        yn_1(i)=w*XN';
        e(i)=Signal_Original_2(i)-yn_1(i);
        w=w+2*u*e(i)*XN;
end

subplot(4,1,3);
plot(Mix_Signal_2);                              %Mix_Signal_1 原始信号
axis([k+1,1000,-10,30]);
title('原始信号');

subplot(4,1,4);
plot(yn_1);                                            %Mix_Signal_1自适应滤波后信号
axis([k+1,1000,-10,30]);
title('自适应滤波后信号');

%****************************************************************************************
%  
%                信号Mix_Signal_1和 Mix_Signal_2  分别作小波滤波
%
%***************************************************************************************

%混合信号Mix_Signal_1  小波滤波
figure(7);
subplot(4,1,1);
plot(Mix_Signal_1);                                %Mix_Signal_1 原始信号
axis([0,1000,-5,5]);
title('原始信号 ');

subplot(4,1,2);
[xd,cxd,lxd] =wden(Mix_Signal_1,'sqtwolog','s','one',2,'db3');
plot(xd);                                                %Mix_Signal_1 小波滤波后信号
axis([0,1000,-5,5]);
title('小波滤波后信号 ');

%混合信号Mix_Signal_2  小波滤波
subplot(4,1,3);
plot(Mix_Signal_2);                                %Mix_Signal_2 原始信号
axis([0,1000,-10,30]);
title('原始信号 ');

subplot(4,1,4);
[xd,cxd,lxd] =wden(Mix_Signal_2,'sqtwolog','h','sln',3,'db3');
plot(xd);                                                %Mix_Signal_2小波滤波后信号
axis([0,1000,-10,30]);
title('小波滤波后信号 ');