基于MATLAB的FIR滤波器的设计及应用（图像去噪）

一、实现功能

1．学会MATLAB的使用，掌握MATLAB的程序设计方法；

2．掌握数字信号处理的基本概念、基本理论和基本方法；

3. 在MATLAB环境下产生噪声信号；

4．掌握MATLAB设计FIR数字滤波器的窗函数设计方法；

5．学会用MATLAB对图像进行分析和处理。

二、设计原理

1.利用窗函数法设计FIR滤波器的原理
如果所希望的滤波器的理想的频率响应函数为，则其对应的单位脉冲响应为

（2.1）

窗函数设计法的基本原理是用有限长单位脉冲响应序列逼近。由于往往是无限长序列，而且是非因果的，所以用窗函数将截断，并进行加权处理，得到：

（2.2）

就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数为

（2.3）

式(2.3)中，N为所选窗函数的长度。

我们知道，用窗函数法设计的滤波器性能取决于窗函数的类型及窗口长度N的取值。设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N 。各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见表2.1。

表2.1 各种窗函数的基本参数

这样选定窗函数类型和长度N之后，求出单位脉冲响应，并按照式（2.3）求出。是否满足要求，要进行演算。一般在尾部加零使长度满足2的整数次幂，以便用FFT计算。如果要观察细节，补零点数增多即可。必要时检查频率响应

开始

读入窗口长度N

计算hd(n)

调用窗函数子程序求w(n)

调用子程序（函数）计算H(k)=DFT[h(n)]

调用绘图子程序（函数）绘制H(k)幅度相位曲线

结束

图2.1主程序框图

计算h(n)= hd(n) w(n)

2．等波纹最佳逼近法设计FIR滤波器的原理

用表示希望逼近的幅度特性函数，要求设计线性相位FIR数字滤波器时，必须满足线性相位约束条件。用表示实际设计的滤波器幅度特性函数。定义加权误差函数为式中，成为误差加权函数，用来控制不同频段（一般指通带和祖带）的逼近精度。等波纹最佳逼近基于切比雪夫逼近，在通带和阻带以的最大值和最小值为准则，采用Remez多重交代迭代算法求解滤波器系数。在Remez多重交换迭代过程中是确知函数。

利用等波纹最佳逼近准则设计线性相位FIR数字滤波器数学模型的建立及其求解算法的推导复杂，所以求解必须借助MATLAB信号处理工具箱函数remezord和remez，只要简单的调用这两个函数就可以完成线性相位FIR数字滤波器的等波纹最佳逼近设计。

如果要求线性相位特性，则还必须满足

根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类，例如，要设计线性相位低通特性，可以选择这一类，而不能选择这一类。

三、主要实验仪器及材料

计算机、MATLAB7.1教学版。

四、实验内容

1.噪声信号的频谱分析。

2.设计数字滤波器和画出频率响应。

利用窗函数和最佳逼近的方法设计FIR滤波器；最后，利用MATLAB中的函数freqz画出各滤波器的频率响应。

3.用滤波器对噪声图像进行滤波。

4.比较几种滤波器去滤波前后噪声图像的波形及频谱。

五、程序设计

I = imread('eight.tif');
J = imnoise(I,'gaussian',0.1);
figure(1) ;imshow(I);title('原始图像');
figure(2);imshow(uint8(J));title('加噪图像');
%窗函数法
fp=1000;fs=2000;Fs=10000;
wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;
B=ws-wp;N=ceil(12*pi/B);wc=(wp+ws)/2/pi;
hn=fir1(N-1,wc);
n=0:N-1;
figure(3)
subplot(1,2,1);stem(n,hn,'.');title('h(n)波形'); grid on;
[hw,w]=freqz(hn,1);
subplot(1,2,2);
plot(w/pi,20*log10(abs(hw)));title('频率响应hw');grid on;
axis([0 1 -100 4]);
yt=filter2(hn,J);
figure(4);imshow(yt,[]);title('滤波后图像1');
%用等波纹最佳逼近方法绘制滤波器损耗函数以及输出信号波形
fp=120;fs=150;Ap=0.2;As=60;Fs=1000;
wc=(fs+fp)/Fs;
B=2*pi*(fs-fp)/Fs;
Nb=ceil(11*pi/B);
fb=[fp,fs];m=[1,0];
dev=[(10^(Ap/20)-1)/(10^(Ap/20)+1),10^(-As/20)];
[Ne,fo,mo,W]=remezord(fb,m,dev,Fs); hn=remez(Ne,fo,mo,W);
yet=filter2(hn,J);
[H,W]=freqz(hn);
m=abs(H);
figure(5);subplot(3,1,1);
plot(W/pi,20*log10(m/max(m)));grid on;
xlabel('\omega/\pi');ylabel('幅度(dB)');
axis([0,1,-100,0]);title('损耗函数曲线');
subplot(3,1,2);plot(hn);title('输出信号波形');
figure(6);imshow(yet,[]);title('滤波后图像2')

六、仿真结果及分析

图6.1　原始图像　　　　　　　　　　　　　　　

图6.2 加噪图像

图6.3 波形和频率响应

图6.4 滤波后图像1

图6.5 损耗函数曲线和输出信号波形

图6.6 滤波后图像2

七、设计总结

1.该课题设计是先画出原始图像，然后画出加噪图像，然后由各项参数求出单位脉冲响应和其频率响应，再分别用布莱克曼窗函数法和等波纹最佳逼近法画出信号波形和去噪图像，并分析和比较，在分别求出损耗函数。

2.FIR滤波器的设计方法方法包括窗函数法、频率采样法和等波纹最佳逼近。本文设计主要采用的是布莱克曼窗函数法和等波纹最佳逼近法设计FIR滤波器对图像进行滤波去噪，通过布莱克曼窗函数这一过程，对滤波前后波形进行对比分析得到结论。等波纹最佳逼近法是一种优化设计方法，它克服了窗函数设计法和频率采样法的缺点，使最大误差最小化，并在整个逼近频段上均匀分布。设计的滤波器的幅频响应在通带和阻带都是等波纹的，而且可以分别控制通带和阻带波纹幅度。

3.经过用窗函数法和等波纹最佳逼近法设计FIR滤波器所得的图形可以看出，窗函数法设计出的滤波器效果更好。

八、参考文献

[1] 陈怀深，吴大正，高西全. MATLAB及在电子信息课程中的应用．3版．北京：电子工业出版社，2006.

[2] 王世一．数字信号处理．北京：北京工业学院出版社，1987.
[3] 胡广书．数字信号处理——理论、算法与实现．北京：清华大学出版社， 1998.

[4] 陈怀深．数字信号处理教程——MATLAB释疑与实现．北京：电子工业出版社，2004.