前言：
转载自http://blog.sina.com.cn/s/blog_45eac6860102uyni.html
信号与系统大作业需要滤除高斯白噪声，苦于MATLAB使用不熟练。
这篇文章质量很高，不仅是对不同滤波方式的简介和感性认识，也是对于MATLAB基本使用的熟悉。

文章目录

• 信号创建
• 传统滤波
• • Butterworth滤波
  • FIR
  • 移动平均滤波
  • 中值滤波
• 现代滤波
• • 维纳滤波
  • 自适应滤波
  • 小波变换

信号创建

%****************************************************************************************%% 创建两个信号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; %构造的混合信号

传统滤波

Butterworth滤波

巴特沃斯滤波器适合用于信号和噪声没有重叠的情况下。

%****************************************************************************************% % 信号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('原始信号1 ');subplot(4,1,2); %Mix_Signal_1 低通滤波滤波后信号 plot(Signal_Filter);axis([0,1000,-4,4]);title('巴特沃斯低通滤波后信号1');%混合信号 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('原始信号2 ');subplot(4,1,4); %Mix_Signal_2 低通滤波滤波后信号 plot(Signal_Filter);axis([0,1000,-10,30]);title('巴特沃斯低通滤波后信号2');

下图是巴特沃斯对两个信号的滤波效果。

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

FIR

情况同巴特沃斯低通滤波器相似。

%****************************************************************************************% % 信号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('原始信号1 ');subplot(4,1,2); %Mix_Signal_1 FIR低通滤波滤波后信号 plot(Signal_Filter);axis([0,1000,-5,5]);title('FIR低通滤波后的信号1');%混合信号 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('原始信号2 ');subplot(4,1,4); %Mix_Signal_2 FIR低通滤波滤波后信号 plot(Signal_Filter);axis([0,1000,-10,30]);title('FIR低通滤波后的信号2');

滤波后的效果如下：

移动平均滤波

移动平均的思路是将离散点附近的点取平均从而实现平滑化。但这种粗糙的平均，对于变换较为平缓的信号进行趋势预测非常有益。

%****************************************************************************************% % 信号Mix_Signal_1 和 Mix_Signal_2 分别作移动平均滤波%%***************************************************************************************%混合信号 Mix_Signal_1 移动平均滤波figure(3);b = [1 1 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 = [1 1 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('中值滤波后的信号');

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

现代滤波

维纳滤波

维纳滤波器是我们这里介绍的第一个现代滤波器。

传统的滤波器只能滤除信号和干扰频带没有重叠的情况，当信号和干扰频带有重叠的时候传统滤波器将无能为力，这时就需要用到现代滤波器。

现代滤波器利用信号和干扰的统计特征（如自相关函数、功率谱等）导出一套**估值算法，然后用硬件或软件予以实现。

维纳滤波是以均方误差最小（Least Mean Square，或简称LMS）为准则的，它根据过去观测值和当前观测值来估计信号的当前值，因此它的解形式是系统的传递函数或单位脉冲响应。

均方误差为：
E [ e 2 ( n ) ] = E [ ( s ( n ) , s ^ ( n ) ) 2 ] E\left[e^2(n)\right] = E\left[\left(s(n), \hat{s}(n)\right)^2\right] E[e2(n)]=E[(s(n),s^(n))2]

MATLAB代码：

%****************************************************************************************% % 信号Mix_Signal_1 和 Mix_Signal_2 分别作维纳滤波%%***************************************************************************************%混合信号 Mix_Signal_1 维纳滤波figure(5);Rxx=xcorr(Mix_Signal_1,Mix_Signal_1); %得到混合信号的自相关函数M=100; %维纳滤波器阶数for i=1:M %得到混合信号的自相关矩阵for j=1:Mrxx(i,j)=Rxx(abs(j-i)+N);endendRxy=xcorr(Mix_Signal_1,Signal_Original_1); %得到混合信号和原信号的互相关函数for i=1:Mrxy(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; %维纳滤波器阶数for i=1:M %得到混合信号的自相关矩阵for j=1:Mrxx(i,j)=Rxx(abs(j-i)+N);endendRxy=xcorr(Mix_Signal_2,Signal_Original_2); %得到混合信号和原信号的互相关函数for i=1:Mrxy(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('维纳滤波后的信号');

多次进行，可以看到统计特征的不利影响，即不稳定。时而效果是并不好的。事实上作者在源程序中对于信号2的处理也并不佳。

自适应滤波

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

%****************************************************************************************% % 信号Mix_Signal_1 和 Mix_Signal_2 分别作自适应滤波%%***************************************************************************************%混合信号 Mix_Signal_1 自适应滤波figure(6);N=1000; %输入信号抽样点数Nk=100; %时域抽头LMS算法滤波器阶数u=0.001; %步长因子%设置初值yn_1=zeros(1,N); %output signalyn_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):NXN=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;endsubplot(4,1,1);plot(Mix_Signal_1); %Mix_Signal_1 原始信号axis([k+1,1000,-4,4]);title('原始信号1');subplot(4,1,2);plot(yn_1); %Mix_Signal_1 自适应滤波后信号axis([k+1,1000,-4,4]);title('自适应滤波后信号1');%混合信号 Mix_Signal_2 自适应滤波N=1000; %输入信号抽样点数Nk=500; %时域抽头LMS算法滤波器阶数u=0.000011; %步长因子%设置初值yn_1=zeros(1,N); %output signalyn_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):NXN=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;endsubplot(4,1,3);plot(Mix_Signal_2); %Mix_Signal_1 原始信号axis([k+1,1000,-10,30]);title('原始信号2');subplot(4,1,4);plot(yn_1); %Mix_Signal_1 自适应滤波后信号axis([k+1,1000,-10,30]);title('自适应滤波后信号2');

自适应滤波的滤波效果如下：

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

小波变换

%****************************************************************************************% % 信号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('小波滤波后信号 ');

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