数字信号处理课程设计窗函数法线性相位FIR数字滤波器设计

1、汕 头 大 学 工 学 院三级项目报告课 程 名 称 ： 数字信号处理 三级项目题目：窗函数法线性相位fir数字滤波器设计指 导 教 师 ： * 系 别： 电子工程系 专 业： 姓名学号完成时间： 2010年 12月12 日 成绩： 评阅人：* 1 内容与要求1、线性相位fir数字滤波器应满足的条件。2、频率选择性线性相位fir数字滤波器的设计算法。2 报告正文一、线性相移fir 数字滤波器的条件：第一类线性相位：由偶对称h(n)=h(n-1-n),可得n=(n1)/2为h(n)的偶对称中心幅度函数：1）h(n)偶对称，n为奇数2）h(n)偶对称，n为偶数第二类线性相位由奇对称h(n)=-h(

2、n-1-n),可得n =(n1)/2为h(n)的奇对称中心幅度函数：3）h(n)奇对称，n为奇数4）h(n)奇对称，n为偶数总结四种线性相位fir的特点：当h(n)为实数且偶对称时，fir为恒相时延，相位曲线是一条过原点、以-(n-1)/2为斜率的直线。信号通过这类滤波器后，各种频率分量的时延都是(n-1)/2。当n为奇数时，时延(n-1)/2是整数，是采样间隔的整数倍，采样点时延后仍是采样点。但当n为偶数时，时延(n-1)/2不是整数，采样点时延后不在采样点位置上。同时，n为偶数时，点是幅度的零点，不能做高通、带阻滤波器。一般情况下，第一类fir特别适合做各种滤波器。当h(n)为实数且奇对称

3、时，fir仅是恒群时延。相位曲线是一条截距为2，以-(n-1)/2为斜率的直线。信号通过该滤波器产生的时延也是(n-1)/2个采样周期，但另外对所有频率分量均有一个附加的90度的相移。程序代码：图一：h=-4 1 -1 -2 5 6 5 -2 -1 1 -4m=length(h);n=0:m-1;hr,w,a,l=hr_type1(h);subplot(2,2,1);stem(n,h);xlabel(n);ylabel(h(n);title(脉冲响应)subplot(2,2,3);stem(0:l,a);xlabel(n);ylabel(a(n);title(a(n)系数)subplot(2,

4、2,2);plot(w/pi,hr);xlabel(频率 单位pi);ylabel(hr);title(1型幅度响应)subplot(2,2,4);pzplotz(h,1);%画极零图图二：h=-4 1 -1 -2 5 6 6 5 -2 -1 1 -4m=length(h);n=0:m-1;hr,w,b,l=hr_type2(h);subplot(2,2,1);stem(n,h);xlabel(n);ylabel(h(n);title(脉冲响应)subplot(2,2,3);stem(1:l,b);xlabel(n);ylabel(b(n);title(b(n)系数)subplot(2,2,2

5、);plot(w/pi,hr);xlabel(频率 单位pi);ylabel(h);title(2型幅度响应)axis(0 2 -20 40);subplot(2,2,4);pzplotz(h,1);%画极零图图三：h=-4 1 -1 -2 5 0 -5 2 1 -1 4m=length(h);n=0:m-1;hr,w,c,l=hr_type3(h);subplot(2,2,1);stem(n,h);xlabel(n);ylabel(h(n);title(脉冲响应)subplot(2,2,3);stem(0:l,c);xlabel(n);ylabel(c(n);title(c(n)系数)sub

6、plot(2,2,2);plot(w/pi,hr);xlabel(频率 单位pi);ylabel(h);title(3型幅度响应)subplot(2,2,4);pzplotz(h,1);%画极零图图四：h=-4 1 -1 -2 5 6 -6 -5 2 1 -1 4m=length(h);n=0:m-1;hr,w,d,l=hr_type4(h);subplot(2,2,1);stem(n,h);xlabel(n);ylabel(h(n);title(脉冲响应)subplot(2,2,3);stem(1:l,d);xlabel(n);ylabel(d(n);title(d(n)系数)subplot

7、(2,2,2);plot(w/pi,hr);xlabel(频率 单位pi);ylabel(h);title(4型幅度响应)subplot(2,2,4);pzplotz(h,1);%画极零图function hr,w,a,l=hr_type1(h)%计算1型滤波器的振幅响应hr（w）%hr=振幅响应%w=数字频率%a=1型低通滤波器系数%l=hr的阶次%h=1型低通滤波器的脉冲响应m=length(h);l=(m-1)/2;a=2*h(l+1:-1:1);n=0:1:l;w=0:1:1000.*pi/500;hr=cos(w*n)*a.;function hr,w,b,l=hr_type2(h)

8、%计算2型滤波器的振幅响应hr（w）%hr=振幅响应%w=数字频率%b=2型低通滤波器系数%l=hr的阶次%h=2型低通滤波器的脉冲响应m=length(h);l=m/2;b=2*h(l:-1:1);n=1:1:l;n=n-0.5;w=0:1:1000.*pi/500;hr=cos(w*n)*b.;function hr,w,c,l=hr_type3(h)%计算3型滤波器的振幅响应hr（w）%hr=振幅响应%w=数字频率%c=2型低通滤波器系数%l=hr的阶次%h=3型低通滤波器的脉冲响应m=length(h);l=(m-1)/2;c=2*h(l+1:-1:1);n=0:1:l;w=0:1:1

9、000.*pi/500;hr=sin(w*n)*c.;function hr,w,d,l=hr_type4(h)%计算2型滤波器的振幅响应hr（w）%hr=振幅响应%w=数字频率%d=4型低通滤波器系数%l=hr的阶次%h=4型低通滤波器的脉冲响应m=length(h);l=m/2;d=2*h(l:-1:1);n=1:1:l;n=n-0.5;w=0:1:1000.*pi/500;hr=sin(w*n)*d.;function pzplotz(b,a)% pzplotz(b,a) 按给定系数向量b,a在z平面上画出零极点分布图% b - 分子多项式系数向量% a - 分母多项式系数向量% a,b

10、向量可从z的最高幂降幂排至z0,也可由z0开始,按z-1的升幂排至z的最负幂.n = length(a); m = length(b); pz = ; zz = ; if (n m) zz = zeros(n-m),1);elseif (m n) pz = zeros(m-n),1); endpz = pz;roots(a); zz = zz;roots(b);pzr = real(pz); pzi = imag(pz);zzr = real(zz); zzi = imag(zz);rzmin = min(pzr,zzr,-1)-0.5; rzmax = max(pzr,zzr,1)+0.5;

11、izmin = min(pzi,zzi,-1)-0.5; izmax = max(pzi,zzi,1)+0.5;zmin = min(rzmin,izmin); zmax = max(rzmax,izmax); zmm = max(abs(zmin,zmax);%uc=exp(j*2*pi*0:1:500/500); % 单位圆plot(real(uc),imag(uc),b,-zmm,zmm,0,0,b,0,0,-zmm,zmm,b);axis(-zmm,zmm,-zmm,zmm);axis(square);holdplot(zzr,zzi,bo,pzr,pzi,rx);holdtext(z

12、mm*1.1,zmm*0.95,z-平面)xlabel(实轴);ylabel(虚轴)title(零极点图)二、频率选择性线性相位fir数字滤波器的设计算法及设计步骤fir数字滤波器设计的基本思想(以窗函数法为例)如果希望得到的滤波器的理想特性为 hd(ejw),那么fir 滤波器的设计就是寻求一个系统函数h(z),用其频响去逼近 hd(ejw)。理想特性的单位冲激响应可由傅里叶反变换求取, 即，为了使hd(n)变成有限长, 可将hd ( n) 利用有限长度的窗函数w(n)直接截短成有限项n , 即。基于kaiser窗的fir数字滤波器的设计步骤1.有理想低通滤波器的频率响应得到理想的hd(n)

13、:2.把hd(n)进行截断移位得到有限长、因果的物理可实现的单位抽样响应h(n)即3.根据技术指标计算窗函数的设计参数n和4.将h(n)乘上kaiser窗函数，既可得到基于kaiser窗的fir数字滤波器5.滤波器的技术指标频带指标依采样频率的归一化频率起伏指标简化为对于较小的进一步的工程简化二、matlab程序代码clear all clfy=klh(1,20000,4000,5000,0.1,80); plot(y)title(h(n) pause clf h,w=freqz(y,1,512); f=w*pi;plot(f,20*log10(abs(h)xlabel(f/khz)title

14、(frequency response of fir filter)clfy=klh(-1,20000,5000,4000,0.1,80);plot(y)title(h(n)pauseclfh,w=freqz(y,l,512);f=w*pi;plot(f,20*log10(abs(h)xlabel(f/khz)title(frequency response of fir filter)clfy=kbp(20000,5000,8000,4000,9000,0.1,80,1);plot(y)h,w=freqz(y,1,512);f=w*pi;plot(f,20*log10(abs(h) xlab

15、el(f/khz)title(frequency response of fir filter)%kaiser窗%alpha（即）和n分别是kaiser窗的控制参数function w=kwind(alpha,n)m=(n-1)/2;den=io(alpha);for n=0:n-1, w(n+1)=io(alpha*sqrt(n*(n-1-n)/m)/den;end%io为bessel方程function s=io(x)eps=10(-9);n=1;s=1;d=1;while d(eps*s), t=x/(2*n); n=n+1; d=d*t2; s=s+d;end%计算kaiser窗控制参

16、数function w=kwind(alpha,n)m=(n-1)/2;den=io(alpha);for n=0:n-1, w(n+1)=io(alpha*sqrt(n*(n-1-n)/m)/den;end%移位并截断理想低/高通滤波器单位冲激响应%wc为截止频率%n为滤波器长度%s=1为低通，s=-1为高通；function h=dlh(s,wc,n)m=(n-1)/2;for k=-m:m, if k=0, h(k+m+1)=(1-s)/2+s*wc/pi; else h(k+m+1)=s*sin(wc*k)/(pi*k); endend%低/高通基于kaiser窗的fir滤波器单位冲激

17、响应%fs为采样频率，fpass为通带截止频率，fstop为阻带截止频率%function h=klh(s,fs,fpass,fstop,apass,astop)fc=(fpass+fstop)/2;wc=2*pi*fc/fs;df=s*(fstop-fpass); df=df/fs;dpass=(10(apass/20)-1)/(10(apass/20)+1);dstop=10(-astop/20);d=min(dpass,dstop);a=-20*log10(d);alpha,n=kparm(df,a);h=dlh(s,wc,n).*kwind(alpha,n);end%带通滤波器原理为高

18、低通滤波器的级联%fs为采样频率，fpa为高通的通带截止频率%fpb为低通的通带截止频率%fsa,fsb分别为高/低阻带截止频率function h=kbp(fs,fpa,fpb,fsa,fsb,apass,astop,s)df=min(fpa-fsa,fsb-fpb);df=df/fs;fa=(1+s)*fpa+(1-s)*fsa-s*df)/2;wa=2*pi*fa/fs;fb=(1+s)*fpb+(1-s)*fsb+s*df)/2;wb=2*pi*fb/fs;dpass=(10(apass/20)-1)/(10(apass/20)+1);dstop=10(-astop/20);d=min

19、(dpass,dstop);a=-20*log10(d);alpha,n=kparm(df,a);h=dbp(wa,wb,n).*kwind(alpha,n);%理想带通滤波器function h=dbp(wa,wb,n)m=(n-1)/2;for k=-m:m, if k=0, h(k+m+1)=(wb-wa)/pi; else h(k+m+1)=sin(wb*k)/(pi*k)-sin(wa*k)/(pi*k); endend%3 结果与分析各个技术指标如下fs=20khz fpass=4khz fstop=5khz apass=0.1db astop=80db通过仿真可以得到基于kais

20、er窗的fir滤波器的单位冲激响应为如下:系统的对数频率响应为：此时n的阶数为103由于阶数n与f成反比，改变f可以实现n的控制；当参数截止频率改为：fpass=4khz fstop=4.5khz 此时计算得n=203滤波器的对数频率响应为：当参数截止频率改为：fpass=4khz fstop=6khz 此时计算得n=203滤波器的对数频率响应为：结果分析从以上基于kaiser窗的低通fir滤波器的结果可以得知：相对于矩形窗和hamming窗，kaiser窗的阶数及形状可以实现控制，而前者两个窗其、apass、astop和d是一定的，实现不了了任意需求的滤波器，而凯瑟窗的各个参数是是可以变化的，可以实现任意技术指标需求的滤波器。kaiser窗是性能取决于窗的阶数n（窗口长度）和其形状控制参数，而阶数n和形状控制参数可以由各个技术指标计算得出，从而实现了满足任意需求的滤波器。通过改变通带截止频率和阻带截止频率，可以实现不改变kaiser窗的形状控制参数而改滤波器的阶数n，上面分别做了fpass=4khz