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

汕 头 大 学 工 学 院三级项目报告课 程 名 称 ： 数字信号处理 三级项目题目：窗函数法线性相位FIR数字滤波器设计指 导 教 师 ： **** 系 别： 电子工程系 专 业： 姓名学号完成时间： 2010年 12月12 日 成绩： 评阅人：**** 1 内容与要求1、线性相位FIR数字滤波器应满足的条件2、频率选择性线性相位FIR数字滤波器的设计算法2 报告正文一、线性相移FIR 数字滤波器的条件：第一类线性相位：由偶对称h(n)=h(N-1-n),可得n=(N–1)/2为h(n)的偶对称中心幅度函数：1）h(n)偶对称，N为奇数2）h(n)偶对称，N为偶数第二类线性相位由奇对称h(n)=-h(N-1-n),可得n =(N–1)/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)为实数且奇对称时，FIR仅是恒群时延相位曲线是一条截距为π／2，以-(N-1)/2为斜率的直线信号通过该滤波器产生的时延也是(N-1)/2个采样周期，但另外对所有频率分量均有一个附加的90度的相移程序代码：图一：h=[-4 1 -1 -2 5 6 5 -2 -1 1 -4]M=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,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 -4]M=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);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 4]M=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)系数')subplot(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 4]M=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(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)%计算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:1000].'*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向量可从z的最高幂降幂排至z^0,也可由z^0开始,按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;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(zmm*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):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('frequency response of FIR Filter')clfy=klh(-1,20000,5000,4000,0.1,80);plot(y)title('h(n)')pauseclf[H,w]=freqz(y,l,512);f=w*pi;plot(f,20*log10(abs(H)))xlabel('f/kHz')title('frequency response of FIR Filter')clfy。