当前位置:文档之家› FIR带通滤波器的设计成熟版

FIR带通滤波器的设计成熟版

FIR带通滤波器的设计成熟版
FIR带通滤波器的设计成熟版

目录

1 技术要求 (1)

2 基本原理 (1)

2.1 FIR滤波器简介 (1)

2.2 窗函数法原理 (2)

3 建立模型描述 (4)

3.1 MATLAB常用函数 (4)

3.1.1 矩形窗函数 (4)

3.1.2 三角窗函数 (5)

3.1.3 广义余弦窗 (5)

3.1.4 汉宁窗(升余弦窗) (6)

3.1.5 fir1函数 (7)

3.1.6 freqz函数 (7)

3.1.7 其他函数与命令 (8)

3.2 方案设计与论证 (8)

3.3 程序流程图 (9)

4 模块功能分析或源程序代码 (9)

5 调试过程及结论 (12)

5.1 实验结果 (12)

5.2 结果分析 (15)

6 思考题 (15)

7 心得体会 (16)

8 参考文献 (17)

FIR 带通滤波器的设计

1 技术要求

用窗函数法设计FIR 带通滤波器。要求低端阻带截止频率ωls =0.2π,低端通带截止频率ωlp =0.35π,高端通带截止频率ωhp =0.65π,高端阻带截止频率ωhs =0.8π。绘出h(n)及其幅频响应特性曲线。

2 基本原理

2.1 FIR 滤波器简介

数字滤波器包括FIR (有限单位脉冲响应)滤波器与IIR (无限单位脉冲响应)滤波器两种。在现代信号处理技术中,例如数据传输、雷达接收以及一些要求较高的电子系统,都越来越多地要求信道具有线性的相位特性。在这方面,FIR 滤波器具有独到的优点,它可以在幅度特性随意设计的同时,保证精确、严格的线性相位特性。

FIR 滤波器的单位脉冲响应h (n )是有限长的(0≤n ≤N-1),其z 变换为1-z 的(N-1)阶多项式:

∑-=-==1

)()()

()(N n n

z n h z X z Y z H

可得FIR 滤波器的系统差分方程为:

∑-=?=-=+--++-+=1

0)

()()()()

1()1()1()1()()0()(N m n x n b m n x m b N n x N b n x b n x b n y

因此,FIR 滤波器又称为卷积滤波器。FIR 滤波器的频率响应表达式为:

∑-=-=1

)()(N n n

j j e n h e H ωω

信号通过FIR 滤波器不失真条件是在通带内具有恒定的幅频特性和线性相位特性。理论上可以证明:当FIR 滤波器的系数满足下列中心对称条件:

)1()(n N h n h --= 或者 )1()(n N h n h ---=

时,滤波器设计在逼近平直幅频特性的同时,还能获得严格的线性相位特性。线性相位FIR

滤波器的相位滞后和群延迟在整个频带上是相等且不变的。对于一个 N 阶的线性相位FIR 滤波器,群延迟为常数,即滤波后的信号简单地延迟常数个时间步长。这一特性使通带频率内信号通过滤波器后仍保持原有波形形状而无相位失真。

FIR 滤波器的设计任务是选择有限长度的h(n),使传输函数)(ω

j e H 满足技术要求。FIR

滤波器的设计方法有多种,如窗函数法、频率采样法及其它各种优化设计方法,本次设计使用窗函数法设计FIR 带通滤波器。

2.2 窗函数法原理

数字滤波器可以理解为是一个计算程序或算法,将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列,并在转化过程中,使信号按预定的形式变化。数字滤波器有多种分类,根据数字滤波器冲激响应的时域特征,可将数字滤波器分为两种,即无限长冲激响应(IIR )滤波器和有限长冲激响应(FIR )滤波器。IIR 数字滤波器具有无限宽的冲激响应,与模拟滤波器相匹配。所以IIR 滤波器的设计可以采取在模拟滤波器设计的基础上进一步变换的方法。FIR 数字滤波器的单位脉冲响应是有限长序列。它的设计问题实质上是确定能满足所要求的转移序列或脉冲响应的常数问题,设计方法主要有窗函数法、频率采样法和等波纹最佳逼近法等。FIR 滤波器具有严格的相位特性,这对于语音信号处理和数据传输是和重要的。目前FIR 滤波器的设计方法主要有三种:窗函数法、频率取样法和切比雪夫等波纹逼近的最优化设计方法。常用的是窗函数法和切比雪夫等波纹逼近的最优化设计方法。

因此设计FIR 滤波器的方法之一可以从时域出发,截取有限长的一段冲击响应作为H(z)的系数,冲击响应长度N 就是系统函数H(z)的阶数。只要N 足够长,截取的方法合理,总能满足频域的要求。一般这种时域设计、频域检验的方法要反复几个回合才能成功。

要设计一个线性相位的FIR 数字滤波器,首先要求理想频率响应)(jw

d e

H 。

)(jw d e H 是w 的周期函数,周期为π2,可以展开成傅氏级数:

∑∞

-∞

=-=

n jwn

d

jw

d e n h e H )()( (公式1-1) 使用上述的传递函数去逼近)(jw

d e H ,一个理想的频率响应)(jw d e H 的傅立叶

反变换:

ω

π

ωπ

d d

e H n h n

j jw

d d )(21

)(20

?=

(公式1-2)

其中)(n h d 是与理想频响对应的理想单位抽样响应序列。但不能用来作为设计FIR DF 用的h(n),因为)(n h d 一般都是无限长、非因果的,物理上无法实现。为了设计出频响类似于理想频响的滤波器,可以考虑用)(n h 来近似)(n h d 。

窗函数的基本思想:先选取一个理想滤波器(它的单位抽样响应是非因果、无限长的),再截取(或加窗)它的单位抽样响应得到线性相位因果FIR 滤波器。这种方法的重点是选择一个合适的窗函数和理想滤波器。

设x(n)是一个长序列,w(n)是长度为N 的窗函数,用w(n)截断x(n),得到N 点序列(n)x n ,即

)()()(n w n x n x n = (公式1-3)

在频域上则有

()

()()()

?--?=

ππ

j j j d e π21e θθωθωW e X X N (公式1-4) 由此可见,窗函数w(n)不仅仅会影响原信号x(n)在时域上的波形,而且也会影响到频域内的形状。

MATLAB 信号工具箱主要提供了以下几种窗函数,如表1

加矩形窗后的频谱和理想频谱可得到以下结论:

加窗使过渡带变宽,过渡带的带宽取决于窗谱的主瓣宽度。矩形窗情况下的过渡带宽是N /4π

。N 越大,过渡带越窄、越陡;

过渡带两旁产生肩峰,肩峰的两侧形成起伏振荡。肩峰幅度取决于窗谱主瓣和旁瓣面积之比。矩形窗情况下是8.95%,与N 无关。工程上习惯用相对衰耗来描述滤波器,相对衰耗定义为:

])0(/)(lg[20])(/)(lg[20)(0H w H e H e H w A j jw == (公式1-5)

这样两个肩峰点的相对衰耗分别是0.74dB 和-21dB 。其中(-0.0895)对应的点的值定义为阻带最小衰耗。

以上的分析可见,滤波器的各种重要指标都是由窗函数决定,因此改进滤波器的关键在于改进窗函数。窗函数谱的两个最重要的指标是:主瓣宽度和旁瓣峰值衰耗。旁瓣峰值衰耗定义为:

旁瓣峰值衰耗=20lg(第一旁瓣峰值/主瓣峰值) (公式1-6) 为了改善滤波器的性能,需使窗函数谱满足:

(1)主瓣尽可能窄,以使设计出来的滤波器有较陡的过渡带。

(2)尽量减少最大旁瓣的相对幅度,也就是能量集中于主瓣,以减小带内、带外波动的最大幅度,增大阻带衰减。

一般来说,以上两点很难同时满足。当选取主瓣宽度很窄时,旁瓣的分量势必增加,从而带内、带外的波动也增加了;当选取最小的旁瓣幅度时,降低了带内、带外的波动,但是过渡带的陡度减小了。所以实际采用的窗函数其特性往往是它们的折中,在保证主瓣宽度达到一定要求的前提下,适当牺牲主瓣宽度来换取旁瓣波动的减小

3 建立模型描述

3.1 MATLAB 常用函数

3.1.1 矩形窗函数

矩形窗(Rectangular Window)函数的时域形式可以表示为:

??

?-≤≤==

,

01

0,1)()(N n n R n w N (公式2-1)

它的频域特性为

()

?

?

? ???

?

? ??=??

?

?

?

--2sin 2sin e

e 21j j ωωωωN W N R (公式2-2) Boxcar 函数:生成矩形窗

调用方式w = boxcar (n):输入参数n 是窗函数的长度;输出参数w 是由窗函数的值组成的n 阶向量。从功能上讲,该函数又等价于w = ones(n,1)。

3.1.2 三角窗函数

三角窗(Bartlett Window )函数时域形式可表示为:

??

??

?---≤≤-=-<<-1)1(2

1

,

122)

1(2

1

01

2)(N n N N n N n N n n w R B (公式2-3)

窗谱为:

()

ωωω

ωωωω??

? ??--??? ??--?????? ????? ????? ??≈??????

? ????? ??????????? ??--=212

212

2s i n 4s i n 22s i n 41s i n 12N j N j j e N N e N N e W (公式2-4)

式中,当N 远大于1时,此时,窗谱主瓣宽度为8π/N 。

3.1.3 广义余弦窗

汉宁窗、海明窗和布莱克曼窗,都可以用一种通用的形式表示,这就是广义余弦窗。这些窗都是广义余弦窗的特例,汉宁窗又被称为余弦平方窗或升余弦窗,海明窗又被称为改进的升余弦窗,而布莱克曼窗又被称为二阶升余弦窗。采用这些窗可以有效地降低旁瓣的高度,但是同时会增加主瓣的宽度。

这些窗都是频率为0、2π/(N –1)和4π/(N –1)的余弦曲线的合成,其中N 为窗的长

度。通常采用下面的命令来生成这些窗:

()()1N π21N :0Ind '

-**-= (公式2-5)

ind)*cos(2C cos(ind)B A Window *+*-= (公式2-6)

其中,A 、B 、C 适用于自己定义的常数。根据它们取值的不同,可以形成不同的窗函数,分别是:

汉宁窗 A=0.5,B=0.5,C=0; 海明窗 A=0.54,B=0.54,C=0; 布莱克曼窗 A=0.5,B=0.5,C=0.08;

3.1.4 汉宁窗(升余弦窗)

汉宁窗(Hanning )函数时域形式可表示为: ())(12cos 121n R N n n W N ???

??

???? ??--=

π

(公式2-7) 利用傅利叶变换的调制特性,由上式可得汉宁窗的平谱函数为: ()

()ωω

πωπωω??

?

?

?--????

??????????? ??-++??? ??--+=21121225.05.0N j R R R j e N W N W W e

W

()ωω??

?

??--=21N j e

W (公式2-8)

式中, ()

()ωωω??

?

??--=21N j R j R e

W e W (公式2-9)

当N 远大于1时,上式可近似表示为:

()()???

??++????????? ?

?-+≈N W N W W W R R R πωπωωω2225.05.0 (公式2-10) 这三部分之和使旁瓣互相抵消,能量更集中在主瓣,汉宁窗函数的最大旁瓣值比主瓣值低31dB ,但是主瓣宽度比矩形窗函数的主瓣宽度增加了1倍,为8π/N 。

hanning 函数:生成汉宁窗 调用方式:

(1) w = hanning(n):输入参数n 是窗函数的长度;输出参数w 是由窗函数的值组成的n 阶向量。

注意:此函数不返回是零点的窗函数的首尾两个元素。 (2) w = hanning(n,'symmetric'):与上面相类似。

(3) w = hanning(n,'periodic'):此函数返回包括为零点的窗函数的首尾两个元素。

3.1.5 fir1函数

设计标准响应FIR 滤波器可使用firl 函数。fir1函数以经典方法实现加窗线性相位FIR 滤波器设计,它可以设计出标准的低通,带通,高通和带阻滤波器。形式为:

b=fir1 (n,Wc ,’ftype’,Window)

各个参数的含义如下:

b —滤波器系数。对于一个n 阶的FIR 滤波器,其n+1个滤波器系数可表示为:

()()()()n z n b z b b z b --++++=1211

n —滤波器阶数;

Wc —截止频率,0≤W c ≤1,Wc=1对应于采样频率的一半。当设计带通滤波器时,Wc=[Wc1 Wc2],Wc 1≤ω≤W c2;

ftype —当指定ftype 时,可设计高通和带阻滤波器。Ftype=high 时,设计高通FIR 滤波器;ftype=stop 时设计带阻FIR 滤波器。低通和带通FIR 滤波器无需输入ftype 参数; Window —窗函数。窗函数的长度应等于FIR 滤波器系数个数,即n+1。

3.1.6 freqz 函数

该函数基于FFT 算法计算数字滤波器Z 变换频率响应。形式为 [h , w] = freqz ( b , a , n ) 返回数字滤波器的n 点复频响应

()

()()()()()()a

b

n j a j n j b j j e n a e a a e n b e b b e H ωωωωω

----++++++=

2121

在简单形式中,b ,a 为滤波器系数,freqz 可得到数字滤波器的n 点复频响应,并将这n 点保存在w 中,相应的频率记录在h 中。

3.1.7 其他函数与命令

设计所用其他函数及命令如表2所示

表2 函数命令及其功能

3.2 方案设计与论证

用窗函数法设计一个FIR带通滤波器。指示如下:

低端阻带截止频率 wls = 0.2*pi;低端通带截止频率 wlp = 0.35*pi;

高端通带截止频率 whp = 0.65*pi;高端阻带截止频率 whs = 0.8*pi;

六种窗函数的基本参数如表3:

表3 窗函数基本参数

以上表格里的参数设置是最佳窗函数设计,根据设计方案的要求,选择一个合适的窗函数进行滤波器的设计,从上表可以看出:最小带阻衰减仅有窗函数决定,不受N的影响,而过渡带的宽度则随窗函数的增加而减小。

3.3 程序流程图

4 模块功能分析或源程序代码

clear; %清除工作空间

close all; %关闭所有打开的窗口

wls=0.2*pi;wlp=0.35*pi; %参数设置

whp=0.65*pi;whs=0.8*pi;

delta_w=min((wlp-wls),(whs-whp)); %求两个过渡带的较小者

wc1=(wls+wlp)/2;wc2=(whp+whs)/2; 截止频率取通带阻带边界频率的均值%矩形窗

N1=ceil(1.8*pi/delta_w); %根据矩形窗精确过渡带宽1.8∏/N计算窗宽

hn1=fir1(N1-1,[wc1,wc2]/pi,boxcar(N1)); %检验设计的滤波器单位脉冲响应[h1,w1]=freqz(hn1,1);

%Hamming窗

N2=ceil(6.6*pi/delta_w); %根据Hamming窗精确过渡带宽6.6∏/N计算窗宽hn2=fir1(N2-1,[wc1,wc2]/pi,hamming(N2));

[h2,w2]=freqz(hn2,1);

%Blackman 窗

N3=ceil(11*pi/delta_w); %根据Blackman窗精确过渡带宽11∏/N计算窗宽hn3=fir1(N3-1,[wc1,wc2]/pi,blackman(N3));

[h3,w3]=freqz(hn3,1);

%Kaiser 窗

N4=ceil(10*pi/delta_w); %根据Kaiser窗技术精确过渡带宽10∏/N计算窗宽hn4=fir1(N4-1,[wc1,wc2]/pi,kaiser(N4));

[h4,w4]=freqz(hn4,1);

%绘图

figure(1) %建立图形窗口

subplot(3,1,1); %把窗口分割成3行1列

n=0:N1-1;stem(n,hn1,'.'); %绘制矩形窗的单位脉冲响应

axis([0,N1-1,-0.4,0.4]); %设置显示范围

xlabel('n');ylabel('h(n)');grid on; %确定x,y轴坐标名称,加网格

title('矩形窗单位冲击响应h(n)'); %添加图形的标题

subplot(3,1,2);

plot(w1/pi,20*log10(abs(h1))); %绘制矩形窗的幅频特性曲线

axis([0,1,-150,5]); %设置显示范围

xlabel('归一化角频率');ylabel('幅度(单位:分贝)');grid on; %确定x,y坐标title('矩形窗幅频响应'); %添加图形的标题

subplot(3,1,3);

plot(w1/pi,180/pi*unwrap(angle(h1))); %绘制矩形窗的相频特性曲线

xlabel('归一化角频率');ylabel('单位:度');grid on;

title('矩形窗相频相应');

figure(2) %建立图形窗口

subplot(3,1,1);

n=0:N2-1;stem(n,hn2,'.'); %绘制Hamming窗单位脉冲响应

axis([0,N2-1,-0.4,0.4]); %确定显示范围

xlabel('n');ylabel('h(n)');grid on;

title('Hamming窗单位脉冲响应h(n)');

subplot(3,1,2);

plot(w2/pi,20*log10(abs(h2))); %绘制Hamming窗幅频响应

axis([0,1,-150,5]);

xlabel('归一化角频率');ylabel('幅度(单位:分贝)');grid on;

title('Hamming窗幅频响应');

subplot(3,1,3);

plot(w2/pi,180/pi*unwrap(angle(h2))); %绘制Hamming窗相频响应

xlabel('归一化角频率');ylabel('单位:度');grid on;

title('Hamming窗相频相应');

figure(3)

subplot(3,1,1);

n=0:N3-1;stem(n,hn3,'.');

axis([0,N3-1,-0.4,0.4]);

xlabel('n');ylabel('h(n)');grid on;

title('Blackman窗单位冲击响应h(n)');

subplot(3,1,2);

plot(w3/pi,20*log10(abs(h3)));

axis([0,1,-150,5]);

xlabel('归一化角频率');ylabel('幅度(单位:分贝)');grid on; title('Blackman窗幅频响应');

subplot(3,1,3);

plot(w3/pi,180/pi*unwrap(angle(h3)));

xlabel('归一化角频率');ylabel('单位:度');grid on;

title('Blackman窗相频相应');

figure(4) %建立图形窗口

subplot(3,1,1);

n=0:N4-1;stem(n,hn4,'.');

axis([0,N4-1,-0.4,0.4]);

xlabel('n');ylabel('h(n)');grid on;

title('Kaiser窗单位脉冲响应h(n)');

subplot(3,1,2);

plot(w4/pi,20*log10(abs(h4)));

axis([0,1,-150,5]);

xlabel('归一化角频率');ylabel('幅度(单位:分贝)');grid on; title('Kaiser窗幅频响应');

subplot(3,1,3);

plot(w4/pi,180/pi*unwrap(angle(h4)));

xlabel('归一化角频率');ylabel('单位:度');grid on;

title('Kaiser窗相频相应');

5 调试过程及结论

5.1 实验结果

新建m文件然后运行程序,检查无误后,运行结果如下:

图1 矩形窗滤波器

图2 Hamming窗滤波器

图3 Blackman窗滤波器

图4 Kaiser窗滤波器

5.2 结果分析

设计结果中,矩形窗和Hamming 窗的情况结果如图1和2所示,Blackman窗和Kaiser 窗的运行结果如图3和4所示。根据理论分析,幅频响应图中,在3db的对应点应该为通带截止频率,在60db处对应的是阻带截止频率。

(1)对于矩形窗:窗宽N=12,h(n)为偶对称,对称中心为n=5.5,由于n为整数,故在n=5和n=6处存在两个极大值;在幅频响应图中,实际设计的低端,高端通带截止频率为0.3262pi和0.6758pi,与技术指标相差6.8%和3.9%。而低端和高端的阻带截止频率为0.167pi和0.833pi,与技术指标相差16.5%和4.13%。其阻带的纹波较大,第一阻带最小衰减27db;由相频特性可以看出,次滤波器是线性相位的。

(2)对于Hamming窗:窗宽N=44,h(n)为偶对称,对称中心为n=21.5,由于n为整数,故在n=21和n=22处存在两个极大值;在幅频响应图中,实际设计的低端,高端通带截止频率为0.29pi和0.705pi,与技术指标相差17.17%和8.47%。而低端和高端的阻带截止频率为0.1953pi和0.8047pi,与技术指标相差2.35%和0.59%。第一阻带最小衰减50db;由相频特性可以看出,次滤波器是线性相位的。

(3)对于Blackman窗:窗宽N=80,h(n)为偶对称,对称中心为n=39.5,由于n为整数,故在n=39和n=40处存在两个极大值;在幅频响应图中,实际设计低端,高端通带截止频率为0.289pi和0.7109pi,与技术指标相差17.4%和9.37%。而低端和高端的阻带截止频率为0.207pi和0.793pi,与技术指标相差3.5%和0.875%。第一阻带最小衰减75db;由相频特性可以看出,滤波器为线性相位的。

(4)对于Kaiser窗:窗宽为N=67,h(n)偶对称,对称中心n=33,有用n为整数,故在n=33处存在一个极大值;在幅频响应图中,实际设计的低端和高端通带的截止频率为0.289pi和0.7109pi,与技术指标相差17.4%和9.37%。而低端和高端的阻带截止频率为0.205pi和0.7969pi,与技术指标相差2.55%和0.39%。第一阻带最小衰减80db;有相频特性可以看出,次滤波器为线性相位的。

6 思考题

窗函数设计法中,选择窗函数的类型与滤波器阶数对滤波器设计的影响?

答:首先明确一个关系式,如果窗宽为N ,则滤波器阶数为N-1,所以滤波器的阶数的增大或减小,相当于窗宽的增大或减小。

FIR 滤波器的设计目标就是尽可能的接近理想滤波器。加窗处理后滤波器的幅度函数为()()()θθωθπωπ

π

d W H H R d -=?-21

。可见只有当窗谱()ωR W 逼近冲激响应时,()ωH 才会

逼近()ωd H 。但这时,窗宽N 趋近于无穷大,相当于没有加窗处理。因此,我们希望是:①窗谱主瓣尽可能窄,以便获得较陡过渡带,而主瓣宽度只跟N 呈反比,主瓣窄要求窗宽N 要大;②尽量减小最大旁瓣的相对幅度,以便增大阻带的衰减,而阻带衰减大则要求窗宽N 要小。

7 心得体会

这已经是我们第三次做课程设计了,可以说是轻车熟路了。这次带通滤波器的设计,从原理上来说,难度并不大。只要把基本参数如通带截止频率,阻带截止频率,通带波动,阻带衰减等掌握其含义,就能很好的完成这次试验。

碰到的难点其一就在于Matlab 的使用。最初的学习使用MATLAB 软件阶段,由于操作不熟练,经常出现函数或者命令输入错误,中英文标点输入没有区分,漏掉棒引号或者分号等情况。虽然都是小错,但是极其容易被忽视而产生错误。通常按照错误要找半天。而慢慢熟练之后只需要调用少量几个函数就能实现设计功能,因此后面的调试过程基本上不存在问题。

当我学会使用Matlab 后,发现一直让我感到枯燥的编程原来并不困难,起码设计带通滤波器是这样,只要调用几个函数就可以轻松完成。在这之后的难点就是结论分析了。虽然作出的图形效果让我有成就感,但这并不是最终结果。一次实验结束后,通过得到的实验数据可以提出或者验证什么理论,才是实验的目的所在。而这恰恰是我目前所缺乏的能力。一直以来,我们所做的实验设计基本都是验证性的,对与不对,一看结果就知道,这直接导致了我忽视了后续的分析,认为只要结果正确就行。就拿这次来说,虽然用不同的窗实现了同一个功能,但用这么多窗的意义何在?当然是为了找到最好的。那如何找到?任务就落在了后续分析上。而经过分析发现,实际上并没有最好的,只有根据不同情况选择最适合的。

尽管现在只是初步学会了简单数字滤波器的设计,离真正掌握还有一定距离,但这段

日子确实令我收益匪浅,这将对我今后的学产生积极的影响。

8 参考文献

[1] 程佩青.数字信号处理教程(第3版).清华大学出版社,2007

[2] 陈亚勇等.MATLAB信号处理详解.人民邮电出版社,2001

[3] 万永革.数字信号处理的MATLAB实现.科学出版社,2007

[4] 王力宁.MATLAB与通信仿真.人民邮电出版社,1999

相关主题
文本预览
相关文档 最新文档