※※※※※※※※※ ※2009级数字信号理 ※
※※
课程设计
※ ※
※※※※
※
※※※※
数字信号处理课程设计报告书
基于MATLAB的FIR滤波器的设计及应用
课题名称 姓 名
学 号
（图像去噪）
院、系、部 电气工程系 专 业 电子信息工程 指导教师
2012年 6 月20日
基于MATLAB的FIR滤波器的设计及应用（图像去噪）
一、设计目的
1．学会MATLAB的使用，掌握MATLAB的程序设计方法
2．掌握数字信号处理的基本概念、基本理论和基本方法 3．在MATLAB环境下产生噪声信号
4．掌握MATLAB设计FIR数字滤波器的窗函数和等波纹最佳逼近设计法 5．学会用MATLAB对图像进行分析和处理
二、设计要求
1．噪声信号的频谱分析；
2．设计数字滤波器和画出频率响应；
利用窗函数和等波纹最佳逼近的方法设计FIR滤波器；最后利用MATLAB 中的函数freqz画出各滤波器的频率响应； 3．用滤波器对噪声图像进行滤波；
4．比较几种滤波器去滤波前后噪声图像的波形及频谱。
三、设计原理及思想
如果所希望的滤波器的理想的频率响应函数为Hd?ej??，则其对应的单位脉
冲响应为
1hd?n??2?????H?e?ed? （3.1） ?jj?d窗函数设计法的基本原理是用有限长单位脉冲响应序列h?n?逼近hd?n?。由于hd?n?往往是无限长序列，而且是非因果的，所以用窗函数??n?将hd?n?截断，并进行加权处理，得到：
h?n??hd?n???n?
（3.2）
h?n?就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数Hej?为
He??????h?n?ej?n?0N?1j?n （3.3）
式中，N为所选窗函数??n?的长度。
2
我们知道，用窗函数法设计的滤波器性能取决于窗函数??n?的类型及窗口长度N的取值。设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N 。各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见表1。
表1 各种窗函数的基本参数
窗函数 旁瓣峰值幅度/dB 过渡带宽 阻带最小衰减/dB 矩形窗 三角形窗 汉宁窗 哈明窗 不莱克曼窗 凯塞窗(α=7.865)
这样选定窗函数类型和长度N之后，求出单位脉冲响应h?n??hd?n????n?，并按照式（3.3）求出Hej?。H?ej??是否满足要求，要进行演算。一般在h?n?尾部加零使长度满足2的整数次幂，以便用FFT计算H?ej??。如果要观察细节，补零点数增多即可。如果H?ej??不满足要求，则要重新选择窗函数类型和长度N，再次验算，直至满足要求。
如果要求线性相位特性，则h?n?还必须满足
-13 -25 -31 -41 -57 -57 4π/N 8π/N 8π/N 8π/N 12π/N 10π/N -12 -25 -44 -53 -74 -80 ??h?n???h?N?1?n?
根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类，例如，要设计线性相位低通特性，可以选择h?n??h?N?1?n?这一类，而不能选择h?n???h?N?1?n?这一类。
四、主程序框图
3
开始 读入窗口长度N 计算hd(n) 调用窗函数子程序求w(n) 计算h(n)= hd(n) w(n) 调用子程序（函数）计算H(k)=DFT[h(n)] 调用绘图子程序（函数）绘制H(k)幅度相位结束 图1 主程序框图
五、程序设计
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;Rp=0.2;As=60;Fs=1000;
4
wc=(fs+fp)/Fs; B=2*pi*(fs-fp)/Fs; Nb=ceil(11*pi/B); fb=[fp,fs];m=[1,0];
dev=[(10^(Rp/20)-1)/(10^(Rp/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')
六、仿真结果及分析
5