窗函数设计FIR滤波器.doc
《窗函数设计FIR滤波器.doc》由会员分享,可在线阅读,更多相关《窗函数设计FIR滤波器.doc(17页珍藏版)》请在冰点文库上搜索。
1.课题描述...................................................................................1
2.题目及要求................................................................................1
3.设计原理....................................................................................1
3.1滤波器的分类....................................................................1
3.2数字滤波器工作原理........................................................1
3.3FIR滤波器的设计指标...................................................33.4窗函数设计FIR滤波器的设计原理.................................5
3.5用窗函数设计滤波器的步骤............................................10
3.6实验所用MATLAB函数说数..........................................11
4设计内容.....................................................................................12
4.1用MATLAB编程实现.....................................................12
4.2结果分析...........................................................................15
5总结...........................................................................................17
6参考文献....................................................................................17
1.课题描述
数字滤波器是指输入、输出均为数字信号,通过数值运算处理改变输入信号所含频率成分的相对比例,或者滤除某些频率成分的数字器件或程序。
因此,数字滤波的概念和模拟滤波相同,只是信号的形成和实现滤波方法不同。
正因为数字滤波通过数值运算实现滤波,所以数字滤波处理精度高、稳定、体积小、质量轻、灵活、不存在阻抗匹配问题,可以实验模拟滤波器无法实现的特殊滤波功能。
本课题使用MATLAB信号处理箱和运用窗函数的FIR滤波器去除无用信号。
2.题目及要求
产生包含三个正弦成分(120hz,80hz,20hz)的信号,设计基于窗函数的FIR滤波器去除120hz,20hz成分,保留80hz信号。
通带允许的最大衰减为0.25dB,阻带应达到的最小衰减为20dB。
滤波器的采样频率为500Hz。
3.设计原理
3.1滤波器的分类
从功能上可以分为:
低通、高通、带通和带阻。
从处理信号分为:
经典滤波器和现代滤波器。
从设计方法上分为:
切比雪夫和巴特沃斯
从实现方法上分为:
FIR和IIR
3.2数字滤波器的工作原理
数字滤波器是一个离散时间系统,输入x(n)是一个时间序列,输出y(n)也是一个时间序列。
如数字滤波器的系统函数为H(Z),其脉冲响应为h(n),则在时间域内存在下列关系
y(n)=x(n)h(n)
在Z域内,输入输出存在下列关系
Y(Z)=H(Z)X(Z)
式中,X(Z),Y(Z)分别为输入x(n)和输出y(n)的Z变换。
同样在频率域内,输入和输出存在下列关系
Y(jw)=X(jw)H(jw)
式中,H(jw)为数字滤波器的频率特性,X(jw)和Y(jw)分别为x(n)和y(n)的频谱。
w为数字角频率,单位rad。
通常设计H(jw)在某些频段的响应值为1,在某些频段的响应为0.X(jw)和H(jw)的乘积在频率响应为1的那些频段的值仍为X(jw),即在这些频段的振幅可以无阻碍地通过滤波器,这些频带为通带。
X(jw)和H(jw)的乘积在频段响应为0的那些频段的值不管X(jw)大小如何均为零,即在这些频段里的振幅不能通过滤波器,这些频带称为阻带。
一个合适的数字滤波器系统函数H(Z)可以根据需要输入x(n)的频率特性,经数字滤波器处理后的信号y(n)保留信号x(n)中的有用频率成分,去除无用频率成分。
3.3FIR滤波器的设计指标
我们在进行滤波器设计时,需要确定其性能指标。
一般来说,滤波器的性能要求往往以频率响应的幅度特性的允许误差来表征。
以低通滤波器特性为例,频率响应有通带、过渡带及阻带三个范围。
在通带内:
1-AP≤≤1≤
在阻带中:
≤≤≤
其中为通带截止频率,为阻带截止频率,Ap为通带误差,为阻带误差。
图2-6低通滤波器的幅度特性
与模拟滤波器类似,数字滤波器按频率特性划分为低通、高通、带通、带阻、全通等类型,由于数字滤波器的频率响应是周期性的,周期为2π。
由于频率响应的周期性,频率变量以数字频率来表示,所以数字滤波器设计中必须给出抽样频率。
图2-7为各种数字滤波器理想幅度,可以看出:
1、一个高通滤波器相当于一个全通滤波器减去一个低通滤波器。
2、一个带通滤波器相当于两个低通滤波器相减。
3、一个带阻滤波器相当于一个低通滤波器加上一个高通滤波器。
这里的相加相减都是相当于并联结构。
图2-7中所示的各种数字滤波器理想频率响应只表示了正频率部分,这样的理想频率响应是不可能实现的,原因是频带之间幅度响应是突变的,因而其单位抽样响应是非因果的。
因此要给出实际逼近容限。
数字滤波器的系统函数,它在z平面单位圆上的值为滤波器频率响应,表征数字滤波器频率响应特征的三个参量是幅度平方响应、相位响应和群延时响应。
窗函数的设计指标主要为:
过渡带宽和阻带最小衰减。
3.4窗函数设计FIR滤波器的设计原理
FIR滤波器与IIR滤波器特点不同,设计方法也就不同。
由于FIR系统的冲激响应就是其系统函数各次项的系数,所以设计FIR滤波器的方法之一可以从时域出发,截取有限长的一段冲激响应作为H(z)的系数,冲激响应长度N就是系统函数H(z)的阶数。
只要N足够长,截取的方法合理,总能满足频域的要求。
这种时域设计、频域检验的方法一般要反复几个回合,不像IIRDF设计靠解析公式一次计算成功。
窗函数法设计FIR的基本思想是:
首先根据给定的设计指标求出理想滤波器的频响,其对应的单位样值响应是非因果的无限长序列。
设计要用一个有限长序列来逼近它,最有效的办法是用一个有限长的窗函数截取理想滤波器的单位样值响应,因而窗函数的形状及长度的选择就成为了关键。
在Matlab中常用的窗函数有矩形窗、Hanning窗、Hamming窗、Blackman窗、Kaiser窗等。
这些窗函数各有优缺点,所以要根据实际情况合理选择窗函数类型。
3.4.1.窗函数分为:
矩形窗、三角形窗、汉宁窗(Hanning)、哈明窗、布莱克曼窗、凯塞---贝塞尔窗。
3.4.2.窗函数法设计原理
设数字滤波器的传输函数为,是与其对应的单位脉冲响应,为系统函数。
(式3.1.1)
(式3.1.2)
(式3.1.3)
一般说来,是无限长的,需要求对的一个逼近。
采用窗函数设计法时,可通过对理想滤波器的单位采样响应加窗设计滤波器
(式3.1.4)
其中,是一个长度有限的窗,在区间0≤n≤N外值为0,且关于中间点对称
(式3.1.5)
频率响应根据(式3.1.5),由卷积定理得出
(式3.1.6)
理想的频率响应被窗函数的离散时间傅立叶变换“平滑”了。
采用窗函数设计法设计出来的滤波器的频率响应对理想响应的逼近程度,由两个因素决定:
①主瓣的宽度;②旁瓣的幅度大小。
理想的情况是主瓣的宽度窄,旁瓣的幅度小。
但对于一个长度固定的窗函数来说,这些不能独立地达到最小。
窗函数的一些通用性质为:
1、窗函数的长度N增加,主瓣的宽度减小,使得过渡带变小。
关系为:
NB=C其中:
B是过渡带的宽度;C是取决于窗函数的一个参数。
如矩形窗为4π。
调整N可以有效地控制过渡带的宽度,但N的改变不改变主瓣和旁瓣的相对比例。
随着N值增加,过渡带变窄,波动频率也随着增加,虽然总的幅度有所减少,但截止频率附近的肩峰并不减少,而只是随着N值的增加,肩峰被抑制在愈来愈小的范围内,使肩峰宽度变窄。
2、窗函数的旁瓣的幅度大小取决于窗函数的选择。
选择恰当的窗函数使主瓣包含更多的能量,相应旁瓣的幅度就减小。
旁瓣幅度的减小,可以减少通带和阻带的波动,使通带尽可能趋近水平,阻带尽可能达到最大衰减。
但通常此时过渡带会变宽。
3、取不同的窗函数对幅度特性的整形效果比单纯的增加窗口长度要强得多。
3.4.3设计方法
这种方法也叫傅里叶级数法。
一般是先给出所要求的理想的滤波器的频率响应,要求设计一个FIR滤波器频率响应来逼近。
设计是在时域进行的,因而先由的傅里叶反变换导出,即
(式3.2.1)
由于是矩形频率响应特性,故一定是无限长序列,且是非因果的,而FIR滤波器的必然是有限长的,所以要用有限长的来逼近无限长的,最有效的方法是截断或者说用一个有限长度的窗口函数序列来截取,即
(式3.2.2)
因而窗函数序列的形状及长度的选择就是关键。
我们以一个截止频率为的线性相位的理想矩形幅度特性的低通滤波器为例来讨论。
设低通特性的群延时为,即
(式3.2.3)
这表明,在通带≤范围内,的幅度是均匀的,其值为1,相位是。
利用
(1)式可得
(式3.2.4)
是中心点在的偶对称无限长非因果序列,要得到有限长的,一种最简单的方法就是取矩形窗,即
但是按照线形相位滤波器的约束,必须是偶对称的,对称中心应为长度的一半(N-1)/2,因而必须=(N-1)/2,所以有
(式3.2.5)
将(式3.2.4)代入(式3.25),可得
(式3.2.6)
此时,一定满足这一线性相位的条件。
下面求的傅里叶变换,也就是找出待求FIR滤波器的频率特性,以便能看出加窗处理后究竟对频率响应有何影响。
按照复卷积公式,在时域是相乘、频域上是周期性卷积关系,即
(式3.2.7)
因而逼近的好坏,完全取决于窗函数的频率特性。
窗函数的频率特性为
(式3.2.8)
对矩形窗,则有
(式3.2.9)
也可表示成幅度函数与相位函数
(式3.2.10)
其中
(式3.2.11)
就是频域抽样内插函数,其幅度函数在之内为一个主瓣,两侧形成许多衰减振荡的旁瓣,如果将理想频率响应也写成
(式3.2.12)
则其幅度函数为
(式3.2.13)
3.5用窗函数设计滤波器的步骤
1、根据阻带的衰减,选择合适的窗:
不同的窗有不同的性质:
不同的窗函数,产生泄漏的大小不一样,频率分辨能力也不一样。
信号的截断产生了能量泄漏,而用FFT算法计算频谱又产生了栅栏效应,从原理上讲这两种误差都是不能消除的,但是我们可以通过选择不同的窗函数对它们的影响进行抑制。
(矩形窗主瓣窄,旁瓣大,频率识别精度最高,幅值识别精度最低;布莱克曼窗主瓣宽,旁瓣小,频率识别精度最低,但幅值识别精度最高)。
2、根据窗函数得到的序列经过firl或fir2得到一个滤波器传输函数系数的序列。
1)fir1:
用来设计传统的低通,高通,带通,带阻,多频带FIR滤波器;
调用格式:
b=fir1(N,Wn);
b=fir1(N,Wn,‘high’);
b=fir1(N,Wn,‘stop’);
参数说明:
N:
阶次,滤波器长度为N+1;
Wn:
通带截止频率,其值在0~1之间,1对应Fs/2;
b:
滤波器系数。
在上述所有格式中,若不指定窗函数的类型,fir1自动选择Hamming窗。
2)fir2:
用来设计具有任意幅度响应的FIR滤波器。
调用格式:
b=fir2(N,F,M);
参数说明:
F是频率向量,其值在0~1之间;
M是和F相对应的所希望的幅频相应。
如同fir1,缺省时自动选用Hamming窗。
3)为了观测到设计出来的滤波器的特性,用freqz得到频率响应。
其中在画频率响应的时候我们分为幅度和相位画出。
又因为我们要观测的是衰减的大小程度,以dB为单位,所以我们在画幅度的时候纵坐标应该转换成dB。
4)为了观测是否滤除已知频率,用filte(b,1,a)函数来实现,对信号的滤波实验。
3.6实验所用MATLAB函数说数
1[H,w]=freqz(b,a,N)
b和a分别为离散系统的系统函数分子、分母多项式的系数向量,返回量H则包含了离散系统频响在0~pi范围内N个频率等分点的值(其中N为正整数),w则包含了范围内N个频率等分点。
调用默认的N时,其值是512。
可以先调用freqz()函数计算系统的频率响应,然后利用abs()和angle()函数及plot()函数,绘制出系统的频响曲线。
2Wn=kaiser(N,beta)
列向量wn中返回长度为N的凯塞——贝塞尔窗函数w(n)。
3ceil(x)
是取大于等于x的最小整数。
4fir1
使用窗函数法设计线性相位FIR数字滤波器的工具箱函数。
本函数在3.5中有详细介绍。
4.设计内容
4.1用MATLAB设计程序如下
clear;fs=500;t=(1:
250)/fs;
x=10*cos(2*pi*20*t)+cos(2*pi*80*t)+10*cos(2*pi*120*t);
L=length(x);N=2^(nextpow2(L));Hw=fft(x,N);
figure
(2);subplot(2,1,1);plot(t,x);
gridon;title('滤波前信号x');xlabel('时间/s');%原始信号
subplot(2,1,2);plot((0:
N-1)*fs/L,abs(Hw));%查看信号频谱
gridon;title('滤波前信号频谱图');xlabel('频率/Hz');ylabel('振幅|H(e^jw)|');
Ap=0.25;As=20;%定义通带及阻带衰减
dev=[10^(-As/20),(10^(Ap/20)-1)/(10^(Ap/20)+1),10^(-As/20)];%计算偏移量
mags=[0,1,0];%带通
fcuts=[30,50,90,110];%边界频率
[N,Wn,beta,ftype]=kaiserord(fcuts,mags,dev,fs);%估算FIR滤波器阶数
hh2=fir1(N,Wn,ftype,kaiser(N+1,beta));%FIR滤波器设计
[hn,w1]=freqz(hh2,1,512);%求解数字滤波器的频率响应
figure(5)
stem(hh2)
title('hh2波形')
figure
(1)%绘图
subplot(2,1,1)
plot(w1/pi,20*log10(abs(hn)))
grid
xlabel('频率w');ylabel('幅度/db');
subplot(2,1,2)
plot(w1/pi,angle(hn))
grid
xlabel('频率w');ylabel('相位/rad');
y=cos(2*pi*80*t)
y=filter(hh2,1,x);%滤波
y(1:
ceil(N/2))=[];%群延时N/2,删除无用信号部分
L=length(y);N=2^(nextpow2(L));Hw_2=fft(y,N);
figure(3);subplot(2,1,1);plot(t(1:
L),y);
gridon;title('y=cos(2*pi*80*t)');xlabel('时间/s');
subplot(2,1,2);plot((0:
N-1)*fs/L,abs(Hw_2));%查看信号频谱
gridon;title('滤波后信号y频谱图');xlabel('频率/Hz');ylabel('振幅|H(e^jw)|');
4.2结果分析
实验的第一幅图是滤波器的频率、相位图,第二幅是产生的原始信号的图,后一幅图是滤波后的图,从图上的信号变化可以看出,基本以达到滤波的要求,达到实验目的。
结果的不理想可能是在选取参数上有些不同。
5总结
本次课程设计对我的收获是很大的,经过两天的自我学习及设计滤波器,让我对滤波器有了更深刻的了解,编程中的错误也为我敲响了我的失误点,让我及时改正,编程的不断修改和完善,是对知识了解的一种透彻表现。
对不懂得函数查资料,也是一种对知识了解得渴望。
本次课题的顺利完成对我的学习起到了很大的推进作用。
6参考文献
《数字信号处理》(第三版)高西全丁美玉编著西安电子科技大学出版社
17