窗函数法设计FIR数字低通滤波器设计毕业设计.docx
《窗函数法设计FIR数字低通滤波器设计毕业设计.docx》由会员分享,可在线阅读,更多相关《窗函数法设计FIR数字低通滤波器设计毕业设计.docx(16页珍藏版)》请在冰豆网上搜索。
窗函数法设计FIR数字低通滤波器设计毕业设计
基于窗函数法的FIR数字低通滤波器设计
摘要
数字滤波器是一种用来过滤时间离散信号的数字系统,通过对抽样数据进行数学处理来达到频域滤波的目的。
根据其单位冲激响应函数的时域特性可分为两类:
无限冲激响应(IIR)滤波器和有限冲激响应(FIR)滤波器。
与IIR滤波器相比,FIR的实现是非递归的,总是稳定的;更重要的是,FIR滤波器在满足幅频响应要求的同时,可以获得严格的线性相位特性。
因此,它在高保真的信号处理,如数字音频、图像处理、数据传输、生物医学等领域得到广泛应用。
滤波器的设计是信号处理的核心问题之一。
根据FIR滤波器的原理,提出了FIR滤波器的窗函数设计法,给出了在MATLAB环境下,用窗函数法设计FIR滤波器的过程和设计实例。
通过利用不同的窗函数方法设计FIR滤波器,对所设计的滤波器进行分析比较,得出各种方法设计的滤波器的优缺点及其不同的使用场合,从而可以在设计滤波器时能够正确的选择FIR数字滤波器的窗函数的选取及设计方法。
关键词:
FIR滤波器,MATLAB,窗函数
1概述
1.1FIR滤波器简介
FIR数字滤波器设计最简单的方法是窗函数法,通常也称为傅立叶级数法。
它是在时域进行的,因而必须由理想滤波器的频率响应
推导出其单位冲激响应
,在设计一个FIR数字滤波器的单位冲激响应
去逼近
。
根据冲激响应的时域特性,数字滤波器可分为无限长冲激响应(IIR)和有限长冲激响应滤波器(FIR),FIR的突出优点是:
系统总是稳定的、易于实现线性相位、允许设计多通带(或多阻带)滤波器,但与IIR相比,在满足同样阻带衰减的情况下需要的阶数较高,滤波器的阶数越高,占用的运算时间越多,因此在满足指标要求的情况下应尽量减少滤波器的阶数。
1.2窗函数设计法
窗函数设计法是一种通过截短和计权的方法使无限长非因果序列成为有限长脉冲应响应序列的设计方法,通常在设计滤波器之前,应该先根据具体的工程应用确定滤波器的技术指标,在大多数实际应用中,数字滤波器常常被用来实现选频操作,所以指标的形式一般为在频域中以分贝值给出的相对幅度响应和相位响应。
用窗函数法设计FIR滤波器的步骤如下:
(1)根据过渡带宽及阻带衰减要求,选择窗函数的类型并估计窗口长度N(或阶数M=N-1),窗函数类型可根据最小阻带衰减As独立选择,因为窗口长度N对最小阻带衰减As没有影响,在确定窗函数类型以后,可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N,设待求滤波器的过渡带宽为Δω,它与窗口长度N近似成反比,窗函数类型确定后,其计算公式也确定了,不过这些公式是近似的,得出的窗口长度还要在计算中逐步修正,原则是在保证阻带衰减满足要求的情况下,尽量选择较小的N,在N和窗函数类型确定后,即可调用MATLAB中的窗函数求出窗函数W(n)。
(2)根据待求滤波器的理想频率响应求出理想单位脉冲响应hd(n),如果给出待求滤波器频率应为Hd(ejw),则理想的单位脉冲响应可以用下面的傅里叶反变换式求出:
(1)
在一般情况下,hd(n)是不能用封闭公式表示的,需要采用数值方法表示;从w=0到w=2π采样N点,采用离散傅里叶反变换(IDFT)即可求出。
(3)计算滤波器的单位脉冲响应h(n),它是理想单位脉冲响应和窗函数的乘积。
(4)算技术指标是否满足要求,为了计算数字滤波器在频域中的特性,可调用freqz子程序,如果不满足要求,可根据具体情况,调整窗函数类型或长度,直到满足要求为止。
使用窗函数法设计时要满足以下两个条件:
窗谱主瓣尽可能地窄,以获得较陡的过渡带;
尽量减少窗谱的最大旁瓣的相对幅度,也就是使能量尽量集中于主瓣,减小峰肩和纹波,进行增加阻带的衰减。
窗函数的选择原则是:
[1]具有较低的旁瓣幅度,尤其是第一旁瓣的幅度;
[2]旁瓣的幅度下降的速率要快,以利于增加阻带的衰减;
[3]主瓣的宽度要窄,这样可以得到比较窄的过渡带。
通常上述的几点难以同时满足。
实际中设计FIR数字滤波器往往要求是线性相位的,因此要求
满足线性相位的条件,即要求w(n)满足:
(2)
所以,窗函数不仅有截短的作用,而且能够起到平滑的作用,在很多领域得到了应用。
2设计原理
2.1基本原理
设计低通FIR数字滤波器,寻求一系统函数
,使其频率响应
逼近滤波器要求的理想频率响应
,其对应的单位脉冲响应
(3)
其中
如果所希望的滤波器的理想的频率响应函数为
,则其对应的单位脉冲响应为
(4)
窗函数设计法的基本原理是用有限长单位脉冲响应序列
逼近
。
由于
往往是无限长序列,而且是非因果的,所以用窗函数
将
截断,并进行加权处理,得到:
(5)
就作为实际设计的FIR数字滤波器的单位脉冲响应序列,其频率响应函数
为
(6)
式中,N为所选窗函数
的长度。
2.2典型的窗函数
(1)矩形窗(RectangleWindow)
(7)
其频率响应和幅度响应分别为:
,
(8)
(2)三角形窗(BartlettWindow)
(9)
其频率响应为:
(3)汉宁(Hanning)窗,又称升余弦窗
(10)
其频率响应和幅度响应分别为:
(4)汉明(Hamming)窗,又称改进的升余弦窗
(11)
其幅度响应为:
(5)布莱克曼(Blankman)窗,又称二阶升余弦窗
(12)
其幅度响应:
(6)凯塞(Kaiser)窗
(13)
其中:
β是一个可选参数,用来选择主瓣宽度和旁瓣衰减之间的交换关系,一般说来,β越大,过渡带越宽,阻带越小衰减也越大。
I0(·)是第一类修正零阶贝塞尔函数。
若阻带最小衰减表示为
,β的确定可采用下述经验公式:
若滤波器通带和阻带波纹相等即δp=δs时,滤波器节数可通过下式确定:
式中:
我们知道,用窗函数法设计的滤波器性能取决于窗函数
的类型及窗口长度N的取值。
设计过程中,要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N。
各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见下表1。
表2-1各种窗函数的基本参数
窗函数
旁瓣峰值幅度/dB
过渡带宽
阻带最小衰减/dB
矩形窗
-13
4π/N
-12
汉明窗
-41
8π/N
-53
布莱克曼窗
-57
12π/N
-74
这样选定窗函数类型和长度N之后,求出单位脉冲响应
,并按照式(6)求出
。
是否满足要求,要进行演算。
一般在
尾部加零使长度满足2的整数次幂,以便用FFT计算
。
如果要观察细节,补零点数增多即可。
如果
不满足要求,则要重新选择窗函数类型和长度N,再次验算,直至满足要求。
如果要求线性相位特性,则
还必须满足
(14)
根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。
要根据所设计的滤波特性正确选择其中一类,例如,要设计线性相位低通特性,可以选择
这一类,而不能选择
这一类。
我们在设计滤波器时,希望窗谱主瓣尽可能窄,以获得较陡的过渡带,同时尽可能减小最大旁瓣的相对幅度,我们发现这是不可能的,因为这本身就是一个矛盾体,所以在设计滤波器时只是根据实际情况来选择合适的窗函数。
3几种数字低通滤波器的窗函数设计
3.1采用矩形窗设计FIR数字低通滤波器
图3-1采用矩形窗设计FIR数字低通滤波器的仿真图
3.2采用汉明窗设计FIR数字低通滤波器
图3-2采用汉明窗设计FIR数字低通滤波器的仿真图
3.3采用布莱克曼窗设计FIR数字低通滤波器
图3-3采用布莱克曼窗设计FIR数字低通滤波器的仿真图
图3-4采用布莱克曼窗冲激响应仿真图
图3-5采用布莱克曼窗对含噪语音信号及滤波后的语音信号仿真图
参考文献
[1]从玉良.数字信号处理原理及其MATLAB实现[M].北京:
电子工业出版社.2009.7
[2]胡广书.数字信号处理理论、算法与实现[M].北京:
清华大学出版社.2003,8
[3]万永革.数字信号处理的MATLAB实现[M].北京:
科学出版社.2007.
[4]薛山.MATLAB基础教程[M].北京:
清华大学出版社.2011.
[5]陈怀琛.数字信号处理教程—MATLAB释义与实现[M].北京:
电子工业出版社.2002
[6]缪家鼎,徐文娟.光电技术[M].杭州:
浙江大学出版社,1994:
112-156.
[7]李晓东,张庆红,叶瑾琳.气候学研究的若干理论问题[J].北京大学学报:
自然科学版,1999,35
(1):
101-106.
[8]郑开青.通讯系统模拟及软件[D].北京:
清华大学无线电系,1987.
[9]OnlineComputerLibraryCenter,Inc.HistoryofOCLC[EB/OL].[2000-01-08].http:
//www.oclc.org/about/history/default.htm.
附录
矩形窗设计程序
passrad=0.4*pi;
w2=boxcar(16);
n=1:
1:
16;
hd=sin(passrad*(n-8))./(pi*(n-8));
hd(8)=passrad/pi;
h2=hd.*rot90(w2);
title('designedbyHanningwindow');
[mag2,rad]=freqz(h2);
subplot(2,1,1);
plot(rad,20*log10(abs(mag2)));
gridon;
[h2,w2]=freqz(h2,1,100,2);
subplot(2,1,2);
plot(w2,unwrap(angle(h2)));
gridon;
汉明窗设计程序
passrad=0.4*pi;
w2=hamming(32);
n=1:
1:
32;
hd=sin(passrad*(n-16))./(pi*(n-16));
hd(16)=passrad/pi;
h2=hd.*rot90(w2);
title('designedbyHanningwindow');
[mag2,rad]=freqz(h2);
subplot(2,1,1);
plot(rad,20*log10(abs(mag2)));
gridon;
[h2,w2]=freqz(h2,1,100,2);
subplot(2,1,2);
plot(w2,unwrap(angle(h2)));
gridon;
布莱克曼窗函数设计程序
passrad=0.4*pi;
w2=blackman(64);
n=1:
1:
64;
hd=sin(passrad*(n-32))./(pi*(n-32));
hd(32)=passrad/pi;
h2=hd.*rot90(w2);
title('designedbyHanningwindow');
[mag2,rad]=freqz(h2);
subplot(2,1,1);
plot(rad,20*log10(abs(mag2)));
gridon;
[h2,w2]=freqz(h2,1,100,2);
subplot(2,1,2);
plot(w2,unwrap(angle(h2)));
gridon;
用布莱克曼语音信号进行滤波,并得滤波前后信号的时域波形及频谱程序
clear;clc;
[x,fs]=wavread('D:
\SHE.wav');
sound(x,fs);
ifx/2==0
x=[x'0];
end
fnoise=10000;T=length(x)/fs;
t=T/length(x):
T/length(x):
T;
forj=1:
length(x)
noise(j)=sin(2*pi*fnoise*t(j));
end
x=x+noise';
fp=1000;
fst=1200;
delta_w=2*pi*(fst-fp)/fs;
wc=pi*(fst+fp)/fs;
N=31;
ifN/2==0
N=N+1;
end
tau=(N-1)/2;
forn=1:
N
h(n)=sin(wc*(n-tau))/(pi*(n-tau));
end
h((N-1)/2)=(fst+fp)/fs;
wn=(blackman(N));
h1=h(n).*wn';
subplot(2,1,1);stem(h1);
F=abs(fft(h1));
subplot(2,1,2);
plot(F(1:
(N-1)/2));
X=abs(fft(x));
figure;subplot(2,1,1);
plot(X(1:
(length(x)+1)/2));
y=conv(x,h1);
Y=abs(fft(y));
subplot(2,1,2);
plot(Y(1:
(length(x)+1)/2));
sound(x,fs);
sound(y,fs);