窗函数法线性相位FIR数字滤波器设计.docx

窗函数法线性相位FIR数字滤波器设计.docx

《窗函数法线性相位FIR数字滤波器设计.docx》由会员分享，可在线阅读，更多相关《窗函数法线性相位FIR数字滤波器设计.docx（21页珍藏版）》请在冰豆网上搜索。

窗函数法线性相位FIR数字滤波器设计

汕头大学工学院

三级项目报告

课程名称：

数字信号处理

三级项目题目：

窗函数法线性相位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=[-41-1-2565-2-11-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=[-41-1-25665-2-11-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（[02-2040]）;

subplot（2,2,4）;

pzplotz（h,1）;%画极零图

图三：

h=[-41-1-250-521-14]

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=[-41-1-256-6-521-14]

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.';

functionpzplotz（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）;

end

pz=[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'）;hold

plot（zzr,zzi,'bo',pzr,pzi,'rx'）;hold

text（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程序代码

clearall

clf

y=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（'frequencyresponseofFIRFilter'）

clf

y=klh（-1,20000,5000,4000,0.1,80）;

plot（y）

title（'h（n）'）

pause

clf

[H,w]=freqz（y,l,512）;

f=w*pi;

plot（f,20*log10（abs（H）））

xlabel（'f/kHz'）

title（'frequencyresponseofFIRFilter'）

clf

y=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）））

xlabel（'f/kHz'）

title（'frequencyresponseofFIRFilter'）

%kaiser窗

%alpha（即α）和N分别是kaiser窗的控制参数

functionw=kwind（alpha,N）

M=（N-1）/2;

den=IO（alpha）;

forn=0:

N-1,

w（n+1）=IO（alpha*sqrt（n*（N-1-n））/M）/den;

end

%%%%%%%%%%%%%%%%%%%%%%%

%IO为bessel方程

functionS=IO（x）

eps=10^（-9）;

n=1;S=1;D=1;

whileD>（eps*S）,

T=x/（2*n）;

n=n+1;

D=D*T^2;

S=S+D;

end

%%%%%%%%%%%%%%%%%%%%%%%

%计算kaiser窗控制参数

functionw=kwind（alpha,N）

M=（N-1）/2;

den=IO（alpha）;

forn=0:

N-1,

w（n+1）=IO（alpha*sqrt（n*（N-1-n））/M）/den;

end

%%%%%%%%%%%%%%%%%%%%%%%%%

%移位并截断理想低/高通滤波器单位冲激响应

%wc为截止频率

%N为滤波器长度

%s=1为低通，s=-1为高通；

functionh=dlh（s,wc,N）

M=（N-1）/2;

fork=-M:

M,

ifk==0,

h（k+M+1）=（1-s）/2+s*wc/pi;

else

h（k+M+1）=s*sin（wc*k）/（pi*k）;

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%

%低/高通基于kaiser窗的FIR滤波器单位冲激响应

%fs为采样频率，fpass为通带截止频率，fstop为阻带截止频率

%

functionh=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

%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%

%带通滤波器原理为高低通滤波器的级联

%fs为采样频率，fpa为高通的通带截止频率

%fpb为低通的通带截止频率

%fsa,fsb分别为高/低阻带截止频率

functionh=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（dpass,dstop）;

A=-20*log10（d）;

[alpha,N]=kparm（DF,A）;

h=dbp（wa,wb,N）.*kwind（alpha,N）;

%%%%%%%%%%%%%%%%%%%%%%

%理想带通滤波器

functionh=dbp（wa,wb,N）

M=（N-1）/2;

fork=-M:

M,

ifk==0,

h（k+M+1）=（wb-wa）/pi;

else

h（k+M+1）=sin（wb*k）/（pi*k）-sin（wa*k）/（pi*k）;

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

3结果与分析

各个技术指标如下

fs=20kHzfpass=4kHzfstop=5kHzApass=0.1dBAstop=80dB

通过仿真可以得到基于kaiser窗的FIR滤波器的单位冲激响应为如下:

系统的对数频率响应为：

此时N的阶数为103

由于

阶数N与Δf成反比，改变Δf可以实现N的控制；

当参数截止频率改为：

fpass=4kHzfstop=4.5kHz此时计算得N=203

滤波器的对数频率响应为：

当参数截止频率改为：

fpass=4kHzfstop=6kHz此时计算得N=203

滤波器的对数频率响应为：

结果分析

从以上基于kaiser窗的低通FIR滤波器的结果可以得知：

相对于矩形窗和Hamming窗，kaiser窗的阶数及形状可以实现控制，而前者两个窗其δ、Apass、Astop和D是一定的，实现不了了任意需求的滤波器，而凯瑟窗的各个参数是是可以变化的，可以实现任意技术指标需求的滤波器。

kaiser窗是性能取决于窗的阶数N（窗口长度）和其形状控制参数α，而阶数N和形状控制参数α可以由各个技术指标计算得出，从而实现了满足任意需求的滤波器。

通过改变通带截止频率和阻带截止频率，可以实现不改变kaiser窗的形状控制参数而改滤波器的阶数N，上面分别做了fpass=4kHzfstop=4.5kHz和fpass=4kHzfstop=6kHz两张情况下滤波器的对数频率特性，实现了对滤波器的控制。

从结果上看过渡带越短，滤波器的阶数越高，滤波器的效果越好，越接近理想滤波器。

基于kaiser窗的FIR高通滤波器

技术指标为:

fs=20kHzfpass=5kHzfstop=4kHzApass=0.1dBAstop=80dB

单位冲激响应为：

其对数频率响应为：

基于kaiser窗的FIR带通滤波器

技术指标为:

fs=20kHzfpa=5kHzfpb=8kHzfsa=4kHz

fsb=9kHzApass=0.1dBAstop=80dB

单位冲激响应为：

其对数频率响应为：

4总结

本次三级项目，小组成员都积极配合。

每个人都负责好自己的部分。

但由于我分配任务的失策，总以为把必做题做完大家再做附加题，结果是必做题做完了，但没时间再做附加题了。

大家都很忙，附加题也不知道怎么分配。

当然，小组部分成员还是有对其思考的。

通过此次练习，每个人对课本的认识都更加深刻，对matlab软件的使用也有一定的提高。

5参考文献

[1]A.V.Oppenheim等著，刘树棠译.信号与系统（第二版），西安交通大学出版社，1998.

[2]陈怀琛著，数字信号处理教程——MATLAB释义与实现，电子工业出版社2004。

[3]陈怀琛吴大正高西全著，MATLAB及在电子信息课程中的应用，电子工业出版社2003。

[4]李海涛邓樱著，MATLAB程序设计教程，高等教育出版社，2002

[5]程佩青著，数字信号处理教程，清华大学出版社2001.