--基于Matlab的FIR数字滤波器的设计
专业年级： 2009级通信工程专业
指导老师： 日 期： 2011年12月26号
基于Matlab的FIR的数字滤波器的设计
1、选题依据
数字滤波器精确度高、使用灵活、可靠性高，具有模拟设备所没有的许多优点，已广泛应用于各个学科技术领域，例如数字电视、语音、通信、雷达、声纳、遥感、图像、生物医学以及许多工程应用领域。随着信息时代数字时代的到来，数字滤波技术已经成为一门极其重要的学科和技术领域。以往的滤波器大多采用模拟电路技术，但是，模拟电路技术存在很多难以解决的问题，例如，模拟电路元件对温度的敏感性，等等。而采用数字技术则避免很多类似的难题，当然数字滤波器在其他方面也有很多突出的优点，所以采用数字滤波器对信号进行处理是目前的发展方向。
2、设计要求
(1)通过课程设计把自己在大学中所学的知识应用到实践当中。
(2)深入了解利用Matlab设计FIR数字滤波器的基本方法。
(3)在课程设计的过程中掌握程序编译及软件设计的基本方法。
(4)提高自己对于新知识的学习能力及进行实际操作的能力。
(5)锻炼自己通过网络及各种资料解决实际问题的能力。
3、设计原理
FIR 滤波器具有严格的相位特性，对于信号处理和数据传输是很重要的。目
前 FIR滤波器的设计方法主要有三种：窗函数法、频率取样法和切比雪夫等波
纹逼近的最优化设计方法。常用的是窗函数法和切比雪夫等波纹逼近的最优化设
计方法。本课题中的窗函数法比较简单，可应用现成的窗函数公式，在技术指标
要求高的时候是比较灵活方便的。
如果 FIR 滤波器的为实数，而且满足以下任意条件，滤波器就具有
准确的线性相位
第一种：偶对称，
第二种：奇对称，
对称中心在处。
4、FIR滤波器的基本概念及特点
抽样响应是有限长的，因而滤波器是稳定的系统。因此，FIR滤波器在
通信、图像处理、模式识别等领域都有着广泛的应用。滤波器设计是根据给定滤波器的频率特性，求得满足该特性的传输函数。
有限长单位冲激响应（FIR）滤波器有以下特点：
(1)系统的单位冲激响应在有限个值处不为零；
(2)系统函数在处收敛，极点全部在处（因果系统）；
(3)结构上主要是非递归结构，没有输出到输入的反馈，但有些结构中
(例如频率抽样结构)也包含有反馈的递归部分。
设FIR滤波器的单位冲激响应为一个点序列，则滤波器的系统函数为:
就是说，它有阶极点在处，有个零点位于有限平面的任何位置。
优点 ：
(1)很容易获得严格的线性相位，避免被处理的信号产生相位失真，这一特点在宽频带信号处理、阵列信号处理、数据传输等系统中非常重要；
(2)可得到多带幅频特性；
(3)极点全部在原点（永远稳定），无稳定性问题；
(4)任何一个非因果的有限长序列，总可以通过一定的延时，转变为因果序列， 所以因果性总是满足；
(5)无反馈运算，运算误差小。
缺点：
(1)因为无极点，要获得好的过渡带特性，需以较高的阶数为代价；
(2)无法利用模拟滤波器的设计结果，一般无解析设计公式，要借助计算机辅助设计程序完成。
5、数字滤波器的设计方法
5.1 窗函数法
设计FIR数字滤波器的最简单的方法是窗函数法，通常也称之为傅立叶级数法。FIR数字滤波器的设计首先给出要求的理想滤波器的频率响应，设计一个FIR数字滤波器频率响应，去逼近理想的滤波响应。然而，窗函数法设计FIR数字滤波器是在时域进行的，因而必须由理想的频率响应推导出对应的单位取样响应，再设计一个FIR数字滤波器的单位取样响应去逼近。设计过程如下：
加窗的作用是通过把理想滤波器的无限长脉冲响应乘以窗函数来产生一个被截断的脉冲响应，即并且对频率响应进行平滑。Matlab工具箱提供的窗函数有：矩形窗(Rectangularwindow)、三角窗(Triangular window)、布拉克曼窗(Blackman window)、汉宁窗(Hanningwindow)、海明窗(Hamming window)、凯塞窗(Kaiser window)、切比雪夫窗(Chebyshev window)。
窗函数主要用来减少序列因截断而产生的Gibbs效应。但当这个窗函数为矩形时，得到的FIR滤波器幅频响应会有明显的Gibbs效应，并且任意增加窗函数的长度（即FIR滤波器的抽头数）Gibbs效应也不能得到改善。为了克服这种现象，窗函数应该使设计的滤波器：
(1)频率特性的主瓣宽度应尽量窄，且尽可能将能量集中在主瓣内；
(2)窗函数频率特性的旁瓣趋于的过程中，其能量迅速减小为零。
5.2频率采样法
频率采样法是从频域出发，根据频域采样定理，对给定的理想滤波器的频率响应加以等间隔的抽样 ，得到：
再利用可求得FIR滤波器的系统函数及频率响应。
而在各采样点间的频率响应则是其的加权内插函数延伸叠加的结果。但对于一个无限长的序列，用频率采样法必然有一定的逼近误差，误差的大小取决于理想频响曲线的形状， 理想频响特性变换越平缓， 则内插函数值越接近理想值，误差越小。为了提高逼近的质量，可以通过在频率相应的过渡带内插入比较连续的采样点，扩展过渡带使其比较连续，从而使得通带和阻带之间变换比较缓慢，以达到减少逼近误差的目的。
选取内个采样点的约束条件为：
增大阻带衰减三种方法：
(1)加宽过渡带宽，以牺牲过渡带换取阻带衰减的增加。
(2)过渡带的优化设计。
利用线性最优化的方法确定过渡带采样点的值，得到要求的滤波器的最佳逼近（而不是盲目地设定一个过渡带值）。
(3)增大N。如果要进一步增加阻带衰减，但又不增加过渡带宽，可增加采样点数N。代价是滤波器阶数增加，运算量增加。
直接从频域进行设计，物理概念清楚，直观方便；适合于窄带滤波器设计，这时频率响应只有少数几个非零值，但是截止频率难以控制。
典型应用：用一串窄带滤波器组成多卜勒雷达接收机，覆盖不同的频段，多卜勒频偏可反映被测目标的运动速度。
5.3最优化设计
最优化设计方法是指采用最优化准则来设计的方法。在 FIR DF的最优化设计中，最优化准则有均方误差最小化准则和等波纹切比雪夫逼近(也称最大误差最小化)准则两种。实际设计中，只有采用窗函数法中的矩形窗才能满足前一种最优化准则，但由于吉布斯 (Gibbs )效应的存在，使其根本不能满足设计的要求。为了满足设计的要求，可以采用其它的窗函数来消除吉布斯效应，但此时的设计已经不能满足该最优化准则了。因此，要完成 FIR DF的最优化设计，只能采用后一种优化准则来实现。
5.3.1等波纹切比雪夫逼近准则
在滤波器的设计中，通常情况下通带和阻带的误差要求是不一样的。等波纹切比雪夫逼近准则就是通过对通带和阻带使用不同的加权函数，实现在不同频段(通常指的是通带和阻带)的加权误差最大值相同，从而实现其最大误差在满足性能指标的条件下达到最小值。
尽管窗函数法与频率采样法在FIR数滤波器的设计中有着广泛的应用，但两者不是最优化的设计 。通常线性相位滤波在不同的频带内逼近的最大容许误差要求不同。等波纹切比雪夫逼近准则就是通过通带和阻带使用不同的加权函数，实现在不同频段（通常指的是通带和阻带）的加权误差最大值相同，从而实现其最大误差在满足性能指标的条件下达到最小值，即使得和之间的最大绝对误差最小。
等波纹切比雪夫逼近是采用加权逼近误差，它可以表示为：
(3-1)
其中，为逼近误差加权函数在误差要求高的频段上，可以取较大的加权值，否则，应当取较小的加权值。
尽管按照 FIR 数字滤波器单位取样响应的对称性和的奇、偶性，FIR 数字滤波器可以分为 4 种类型，但滤波器的频率响应可以写成统一的形式：
(3-2)
其中，为幅度函数，且是一个纯实数，表达式也可以写成统一的形式：
(3-3)
其中，为的固定函数，为M个余弦函数的线性组合。
5.3.2仿真函数
利用数字信号处理工具箱中的 remezord 和 remez函数可以实现 FIRDF的最优化设计。在此先介绍这两个函数：
功能:利用 remezord 函数可以通过估算得到滤波器的近似阶数，归一化频率带边界fo ，频带内幅值ao 及各个频带内的加权系数weights。输入参数f为频带边缘频率，a 为各个频带所期望的幅度值，dev是各个频带允许的最大波动。
功能:利用 remez 函数可以得到最优化设计的FIR DF的系数 ，输入参数是滤波器的阶数，fo，ao，weights参数含义说明同(1)。ftype 是所设计的滤波器类型，它除了可以设计普通的滤波器外，它还可以设计数字希尔钞特变换器以及数字微分器。实际设计中，由于 remezord 函数可跑高估或低估滤波器的阶数，因此在得到滤波器的系数后，必须检查其阻带最小衰减是否满足设计要求。如果此时的技术指标不能满足设计要求，则必须提高滤波器的阶数到等。故等波纹切比雪夫逼近法设计FIR数字滤波器的步骤是：
(1)给出所需的频率响应，加权函数和滤波器的单位取样响应
的长度。
(2)由(1)中给定的参数来形成所需的、和的表达式。
(3)根据Remez算法，求解逼近问题。
(4)利用傅立叶逆变换计算出单位取样响应。
6、结果分析
窗函数法
图1 窗函数设计的FIR低通滤波器频率响应
由图1可以看出，采用特殊的窗函数如Hamming窗，可以减小Gibbs效应，但同时也会使滤波器的过度带变宽。波动幅度取决于窗函数幅度频谱旁瓣的相对幅度，而波纹的多少取决于窗函数旁瓣的多少。以上两点是就是窗函数直接截断引起的截断效应在频域的反映，截断效应直接影响滤波器的性能，因为通带内的波动会影响滤波器痛带中的平稳性，阻带内的波动则影响阻带最小衰减，因此，减少截断效应也是FIR数字滤波器设计的关键之一。
频率采样法
图2 FIR的单位取样响应
图3 FIR的低通衰减幅频特性
图2为在间断点处增加一个过渡点后的情况。从图3中可以看出滤波器的带外衰减指标有了明显的改善。但这同时增加了滤波器的过渡带宽。所以，在带外衰减和过渡带宽这两个指标之间需要有一个折衷。
因频率取样点都局限在的整数倍点上，所以在指定通带和阻带截止频率时，这种方法受到，比较死板。充分加大N，可以接近任何给定的频率，但计算量和复杂性增加。频率采样法偏离设计指标明显，阻带衰减最小，只有适当选取过渡带样点值，才会取得较好的衰耗特性。
最优化设计
图4 滤波器输出的幅频及相频响应特性
由图形可以看出，在设计中，如果该滤波器的特性不满足要求，那么，原有参数必须作适当调整。这在程序中很容易实现，只需对参数进行重新设定，就可以得到新条件下滤波器的特性。采用最优化设计方法时大大减小了滤波器的阶数，从而减小了滤波器的体积，并最终降低了滤波器的成本。这样使得设计出来的滤波器更为简单经济。因而在实际的滤波器设计中，这种最优化方法是完全可行的。在实际应用中，如果需要对某一信号源进行特定的滤波，并要检验滤波效果 ，应用传统方法实施起来比较繁琐。在Matlab环境下，可先用软件模拟产生信号源 ，再设计滤波器对其进行滤波 。
同样是设计一个FIR低通数字滤波器，综合分析可以看出：
(1)窗函数法在阶数较低时，阻带特性不满足设计要求，只有当滤波器阶数较高时，使用海明窗和凯塞窗基本可以达到阻带衰耗要求；
(2)频率采样法偏离设计指标最明显，阻带衰减最小，而且设计比采用窗函数法复杂。只有适当选取过渡带样点值，才会取得较好的衰耗特性；
(3)利用等波纹切比雪夫逼近法则的设计可以获得最佳的频率特性和衰耗特性，具有通带和阻带平坦，过渡带窄等优点。
综上所述，FIR滤波器很容易实现具有严格线性相位的系统， 使信号经过处理后不产生相位失真，舍入误差小，而且稳定，因此越来越受到广泛的重视。Matlab软件的诞生，使数字信号处理系统的分析与设计得简单，它已经成为电子工程师必备的一个工具软件。
7、课程设计体会
通过本次设计，感受到了课程设计是一项很需要耐心于细心的工作，在不的
改正与检查中，自我感觉学到了不少，概括为以下几点：
(1)对于Matlab语句有了更加深刻的理解，也注意到了一些运算符号的使
用，例如数组的相乘需用(*)来表示，而一般数字相乘应用。还有当运用数组
的法时，必须保持数组是等长的，否则，不能相加。
(2)想要改变图形的尺寸，可调用AXIS函数。
(3)在编程过程中应该注意一些细节问题，例如中英文符号的区别，往往一些错误都是由于粗心而导致的。
(4)设计过程中，学习了许多数字信号处理课程中关于数字滤波器的设计的内容，再通过利用参考文献与网络，完成了用Matlab进行数字信号处理课程设计。
(5)通过课程设计，加深了对课堂抽象概念的理解，巩固了课堂上所学的理论知识，并能很好地理解与掌握数字信号处理中的基本概念、基本原理、基本分析方法。同时掌握编程方法和解决实际问题的技巧。
(6)与其他高级语言的程序设计相比，Matlab环境下可以更方便、快捷地设计出具有严格线性相位的FIR滤波器，节省大量的编程时间，提高编程效率，且参数的修改也十分方便，还可以进一步进行优化设计。
(7)随着版本的不断提高，Matlab在数字滤波器技术发挥着更大的作用。同时，用Matlab计算有关数字滤波器的设计参数，如等，对于数字滤波器的硬件实现也提供了一条简单而准确的途径和依据。
附 录
窗函数实现程序：
passrad=0.4*pi;
w1=boxcar(61);
w2=hamming(61)
n=1:1:61;
hd=sin(passrad*(n-31))./(pi*(n-31));
hd(31)=passrad/pi;
h1=hd.*rot90(w1);
h2=hd.*rot90(w2);
[mag1，rad]=freqz(h1);
[mag2，rad]=freqz(h2);
subplot(2，2，1);
plot(rad，20*log10(abs(mag1)));
grid on;
title('designed by Rectangular window');
subplot(2，2，2);
plot(rad，20*log10(abs(mag2)));
grid on;
title('designed by Hamming window');
[h1，w1]=freqz(h1，1，100，2);
subplot(2，2，3);
plot(w1，unwrap(angle(h1)));
grid on;
[h2，w2]=freqz(h2，1，100，2);
subplot(2，2，4);
plot(w2，unwrap(angle(h2)));
grid on;
频率采样法实现程序：
f=[0.2 0.3];
a=[1 0];
dev=[0.04 0.02]; %给出滤波器的参数
[n f0 a0 w]=remezord(f，a，dev);
N=n;
alpha=(N-1)/2;
k=0:N-1;
wp=0.2*pi; ws=0.3*pi; %计算理想低通滤波器的截止频率
wc=(wp+ws)/2;
m=fix(wc*N/(2*pi)+1); %在两边过渡带取值为 0.5 的采样点
T = 0.5;
Hrs=[ones(1，m)，T，zeros(1，N-2*m-1)，T，ones(1，m-1)];
k1 = 0:floor(alpha);
k2 = floor(alpha+1):N-1;
phai=[-alpha*(2*pi)/N*k1，alpha*(2*pi)/N*(N -k2)];
H =Hrs.*exp(j*phai); %计算单位冲激响应
h =ifft(H，N);
[h1，w1] = freqz(h，1，256，1);
hr = abs( h1);
h1 = 20* log10(hr); %画出 FIR DF的单位取样响应
figure(1);
k=0:N-1;
stem(k，h，'k.')
axis([0，N-1，1.1*min(real(h))，1.1*max(real(h))]);
xlabel('n'); ylabel('h(n)');
grid on; %画出 FIR DF的低通衰减幅频特性
figure(2);
plot(w1，h1);
xlabel('Normalized')
Frequency(×πrad/sample)');
ylabel('Magnitude(dB)');
grid on;
最优化设计实现程序
fs =2000; %采样频率
rp =3; %通带波纹
rs =40; %阻带波纹
f = [500 600]; %截止频率
a = [1 0]; %期望的幅度
dev= [(10^(rp/20) - 1)/ (10^(rp/20) + 1) 10^( - rs/20) ];
[n ，fo ，ao ，w] = remezord(f ，a ，dev，fs) ;
b = remez(n ，fo ，ao ，w) ; %调用最优设计法中 remez函数
freqz(b ，1 ，1024 ，fs);
