利用MATLAB结合频率取样法设计数字高通FIR滤波器Word文件下载.docx
《利用MATLAB结合频率取样法设计数字高通FIR滤波器Word文件下载.docx》由会员分享,可在线阅读,更多相关《利用MATLAB结合频率取样法设计数字高通FIR滤波器Word文件下载.docx(16页珍藏版)》请在冰豆网上搜索。
1FIR数字滤波器1
1.1FIR滤波器的特点1
1.2FIR数字滤波器设计方法2
1.3线性相位FIR数字滤波器的条件和特点2
1.3.1线性相位条件2
1.3.2线性相位FIR滤波器的幅度特性与相位特性3
2利用频率采样法设计FIR滤波器4
2.1用频率采样法设计滤波器的基本原理4
2.2线性相位的约束条件5
2.3逼近误差及其改进措施6
2.3.1产生误差的原因6
2.3.2减小误差的方法7
2.4频率采样法的特点8
3频率取样法的数字高通滤波器的实现8
3.1MATLAB的介绍8
3.2设计条件8
3.3设计程序9
3.4调试结果11
4心得体会12
附录14
摘要
MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。
它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中,为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案,并在很大程度上摆脱了传统非交互式程序设计语言(如C、Fortran)的编辑模式,代表了当今国际科学计算软件的先进水平。
本文介绍了如何利用MATLAB仿真软件系统及数字信号处理所学知识利用频率采样法设计一个数字高通滤波器。
以此来巩固课堂理论学习,并能用所学理论知识正确分析信号处理的基本问题和解释信号处理的基本现象。
关键字:
MATLAB;
数字信号处理;
数字滤波器;
频率采样法
Abstract
MATLABisreleasedbytheUnitedStatesmathworksmainlyforscientificcomputing,visualizationandinteractiveprogramdesignedhigh-techcomputingenvironment.Itnumericalanalysis,matrixcomputation,scientificdatavisualizationaswellasnon-lineardynamicsystemsmodelingandsimulation,andmanyotherpowerfulintegratedinaneasy-to-useWindowsenvironment,scientificresearch,engineeringdesignandtheneedforeffectivenumericaltheeditmodemanyscientificfieldsprovidesacomprehensivesolution,andinlargeparttogetridofthetraditionalnon-interactiveprogramminglanguage(suchasC,Fortran),onbehalfoftheadvancedleveloftoday'
sinternationalscientificcomputingsoftware.
ThisarticledescribeshowtouseMATLABsimulationsoftwaresystemsanddigitalsignalprocessinglearnedknowledgeusingfrequencysamplingmethodtodesignadigitalhigh-passfilter.Inordertoconsolidatethetheoreticalclassroomlearning,andbasicquestionsandexplainbasicsignalprocessingphenomenoncanbelearnedtheoreticalknowledgetocorrectlyanalyzethesignalprocessing.
Keyword:
digitalsignalprocessing;
digitalfilter;
frequencysamplingmethod
1FIR数字滤波器
1.1FIR滤波器的特点
FIR滤波器的脉冲响应h(n)是有限长的(0≤n≤N-1),其z变换为:
(式1)
它是z-1的(N-1)阶多项式,在有限z平面(0<
n<
∞)上有(N-1)个零点,而极点位于z平面原点z=0处,且有(N-1)阶。
FIR滤波器的基本结构可以理解为一个分节的延时线,把每一节的输出加权累加,可得到滤波器的输出,FIR滤波器的冲激响应h(n)是有限长的,数学上M阶FIR滤波器可以表示为:
y(n)=
(式2)
其系统函数为:
H(z)=
(式3)普通的直接型FIR滤波器结构如图1所示。
图1FIR滤波器的直接型结构
FIR滤波器最突出的优点有2个:
一是只要对h(n)附加一定的条件,很容易获得严格的线性相位特性;
二是由于H(z)的极点位于原点z=0处,始终满足稳定条件,所以FIR滤波器永远稳定。
三是FIR滤波器由于单位脉冲响应是有限长的,因而可以用快速傅里叶变换(FFT)算法来实现过滤信号,从而可大大提高运算效率。
但是,要取得很好的衰减特性,FIR滤波器H(z)的阶次比IIR滤波的要高。
1.2FIR数字滤波器设计方法
IIR滤波器设计中的各种变换法对FIR滤波器设计是不适用的,这是因为那里是利用有理分式的系统函数,而FIR滤波器的系统函数只是z-1的多项式。
FIR的设计任务是选择有限长度的脉冲响应h(n),得到系统函数H(z),使幅频特性满足技术指标要求,同时使相频特性达到线性相位。
常用设计方法:
(1)窗函数法
(2)频率采样法
(3)切比雪夫等波纹逼近法。
人们最感兴趣的是FIR滤波器具有线性相位的相频特性。
对非线性相位的FIR滤波器,一般可以用IIR滤波器来代替,因为同样幅度特性,IIR滤波器所需阶数比FIR滤波器的阶数要少得多。
1.3线性相位FIR数字滤波器的条件和特点
1.3.1线性相位条件
对于长度为N的h(n),传输函数为
(式4)
H(ejω)=Hg(ω)ejθ(ω)(式5)
式中,Hg(ω)称为幅度特性,θ(ω)称为相位特性。
注意,这里Hg(ω)不同于|H(ejω)|,Hg(ω)为ω的实函数,可能取负值,而|H(ejω)|总是正值。
H(ejω)线性相位是指θ(ω)是ω的线性函数,即
θ(ω)=-τω,τ为常数(式6)
如果θ(ω)满足θ(ω)=θ0-τω,θ0是起始相位
严格地说,此时θ(ω)不具有线性相位,但以上两种情况都满足群时延是一个常数,即
(式7)
也称这种情况为线性相位。
1.3.2线性相位FIR滤波器的幅度特性与相位特性
线性相位FIR滤波器的幅度特性与相位特性如下图:
图2线性相位FIR滤波器的幅度特性与相位特性一览表
在设计时,要注意选择合适的h(n)对称形式(奇或偶)和h(n)长度N(奇数或偶数)。
如要设计高通滤波器,只能选情况1和情况4;
要设计低通滤波器,只能选情况1和情况2。
2利用频率采样法设计FIR滤波器
2.1用频率采样法设计滤波器的基本原理
待设计的滤波器的传输函数用Hd(ejω)表示,可按下列思路进行设计:
1它在ω=0到2π之间等间隔采样N点,得到Hd(k)
(式8)
2N点Hd(k)进行IDFT,得到h(n)
(式9)
式中,h(n)作为所设计的滤波器的单位取样响应。
3h(n)求系统函数H(z)
(式10)
将插值公式重写如下
(式11)
此式就是直接利用频率采样值Hd(k)形成滤波器的系统函数。
用频率采样法设计线性相位滤波器的条件:
FIR滤波器具有线性相位的条件是h(n)是实序列,且满足h(n)=h(N–1–n),其传输函数应满足的条件是
(式12)
(式13)
(式14)
(式15)
且Hg(π)=0。
在ω=0~2π之间等间隔采样N点,将ω=ωk代入式(4~7)中,并写成k的函数:
(式16)
(式17)
,N为奇数(式18)
,N为偶数且(式19)
(式20)
说明N等于奇数时Hg(k)对(N–1)/2偶对称,N等于偶数时,Hg(k)对N/2奇对称,且Hg(N/2)=0。
对于高通滤波器,这里N只能取奇数。
截止频率为ωc,采样点数N,Hg(k)和θ(k)用下面公式计算
(式21)
以上是用频率采样法设计滤波器的基本原理。
2.2线性相位的约束条件
以h(n)为偶对称,N为奇数的情况进行分析。
1)FIR的频响具有线性相位的一般表达式
当h(n)为偶对称,N为奇数时,则
(式22)
而且幅度函数H(w)应为偶对称,即
(式23)
2)采样值H(k)具有线性相位的约束
(式24)
其中,
表示采样值的模(纯标量),
表示其相角。
因此,在采样点上具有线性相位的条件应为:
(式25)
而且,
必须满足偶对称,即:
(式26)
实际滤波器的传输函数
,与理想的传输函数Hd(ejω)间存在误差,如图2
图3频率采样的响应
需要讨论逼近误差问题及其改进措施。
2.3逼近误差及其改进措施
2.3.1产生误差的原因
从图3可看出,实际的H(ejω)与理想的Hd(ejω)相比,误差主要体现在一是通带和阻带出现波动,二是过渡带加宽,与窗函数设计法情况类似,产生误差的原因可从时域和频域两方面进行分析。
从时域分析:
如果Hd(ejω)有间断点,那么相应单位取样响应hd(n)应是无限长的。
这样,由于时域混叠,引起所设计的h(n)和hd(n)有偏差。
为此,希望在频域的采样点数N加大。
N愈大,设计出的滤波器愈逼近待设计的滤波器Hd(ejω)。
从频域分析:
在采样点ω=2πk,k=0,1,2,…,N-1,Ф(ω-2πk/N)=1,因此,采样点处H(ejωk)(ωk=2πk/N)与H(k)相等,逼近误差为0。
在采样点之间,H(ejω)由有限项的H(k)Ф(ω-2πk/N)之和形成。
其误差和Hd(ejω)特性的平滑程度有关,特性愈平滑的区域,误差愈小;
特性曲线间断点处,误差最大。
表现形式为间断点用倾斜线取代,且间断点附近形成振荡特性,使阻衰减减小,往往不能满足技术要求。
2.3.2减小误差的方法
最直观的想法是增加采样点数,即加大N值,由于过渡带就等于采样间隔(参看图3),即
(式27)
所以加大N,可使过渡带变窄,但增加要适当,否则会增加滤波器体积与成本。
但是,增加N并不会改善滤波器的阻带衰减特性,因为Hd(ejω)是理想矩形,无论怎样增多频率采样的点数,在通、阻带交界处,幅值总是从1突变到0,会引起较大的起伏振荡。
为使逼近误差更小,和窗口法的平滑截断一样,通过在理想频率响应的不连续点的边缘上加一些过渡的抽样点,减小频带边缘的突变,也就减小了起伏振荡,增大了阻带最小衰减。
一般过渡带取一、二、三点抽样值即可得到满意结果。
如在低通设计中,不加过渡点时,阻带最小衰减为-20dB,加三个过渡点(最优设计)则可达-80dB到-95dB左右。
加过渡点的示意如图4所示。
图4理想低通滤波器增加过渡点
增加过渡点,可使阻带衰减明显提高,但付出的代价是过渡带加宽,可通过下式加大N来调整。
m=0,1,2,3…(式28)
2.4频率采样法的特点
频率采样法设计滤波器最大的优点是直接从频率域进行设计,比较直观,也适合于设计具有任意幅度特性的滤波器。
但边界频率不易控制。
如果增加采样点数N,对确定边界频率有好处,但会增加滤波器的成本。
因此,它适合于窄带滤波器的设计。
3频率取样法的数字高通滤波器的实现
3.1MATLAB的介绍
MATLAB是矩阵实验室(MatrixLaboratory)的简称,是美国MathWorks公司出品的商业数学软件,用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境,主要包括MATLAB和Simulink两大部分。
3.2设计条件
ws=0.6pi,wp=0.8pi,通带波动1dB,阻带衰减50dB,M=33。
3.3设计程序
%频率采样技术:
高通
%ws=0.6pi,wp=0.8pi,Rp=1dB,As=50dB
%M=33,T1=0.1095;
T2=0.598;
M=33;
alpha=(M-1)/2;
l=0:
M-1;
wl=(2*pi/M)*l;
T1=0.1095;
Hrs=[zeros(1,11),T1,T2,ones(1,8),T2,T1,zeros(1,10)];
Hdr=[0,0,1,1];
wdl=[0,0.6,0.8,1];
k1=0:
floor((M-1)/2);
k2=floor((M-1)/2)+1:
angH=[-alpha*(2*pi)/M*k1,alpha*(2*pi)/M*(M-k2)];
H=Hrs.*exp(j*angH);
h=real(ifft(H,M));
[db,mag,pha,grd,w]=freqz_m(h,1);
[Hr,ww,a,L]=Hr_Type1(h);
subplot(1,1,1)
subplot(2,2,1);
plot(wl(1:
17)/pi,Hrs(1:
17),'
o'
wdl,Hdr);
axis([0,1,-0.1,1.1]);
title('
高通:
M=33,T1=0.1095,T2=0.598'
)
xlabel('
'
);
ylabel('
Hr(k)'
set(gca,'
XTickMode'
'
manual'
XTick'
[0;
.6;
.8;
1])
XTickLabelMode'
XTickLabels'
['
0'
;
.6'
.8'
1'
])
YTickMode'
YTick'
[0,0.109,0.59,1]);
grid
subplot(2,2,2);
stem(l,h);
axis([-1,M,-0.4,0.4])
title('
脉冲响应'
h(n)'
text(M+1,-0.4,'
n'
subplot(2,2,3);
plot(ww/pi,Hr,wl(1:
振幅响应'
频率(单位:
pi)'
Hr(w)'
subplot(2,2,4);
plot(w/pi,db);
axis([0,1,-100,10]);
幅度响应'
xlabel('
ylabel('
分贝数'
Manual'
[-50;
0]);
YTickLabelMode'
YTickLabels'
50'
%
3.4调试结果
图5频率采样技术:
高通,最优法
结果分析:
第一幅图为要高通滤波器原型,可以看到它在过渡带添加了两个采样点,以增加阻带衰减;
第二幅图为系统函数单位脉冲响应图形,可以看出,它以中点成偶对称,由于采样点数为奇数,故在对称轴处有取值;
第三幅图(左下)为根据频率取样法设计出的滤波器振幅响应,可以看出它在采样点处的取值与原高通滤波器精确一致,在其他点处与原高通滤波器取值逼近有上下波动;
第四幅图为用分贝数表示的幅度响应,可以看到采用线性最优法设计的高通滤波器的阻带衰减大于50db。
设计取得了良好的效果。
4心得体会
Matlab的课程设计做到现在已经基本接近尾声了,既然学习一门课程,简单的总结是必须要有的。
以前在《信号与系统》和《数字信号处理》的实验中已经接触过matlab,所以上手并不是很难,不过在设计的时候还是遇到了不少问题,首先是对频率取样法掌握的不到位,重新学习了频率取样法后,发现如何利用程序实现频率取样法成了一个问题。
通过自己在网上查找资料,看从图书馆借来的书以及对照着老师的PPT,不断的调试,终于做出了成果。
课程设计虽然做完了,但现在学的这点知识还远远不够,特别是这个软件的函数非常多,要能够熟练运用我们还有很多要学习。
不过我觉得Matlab的函数设计都比较合理,她总是从函数本身的意义出发命名,这使我们记不会很难。
总之这次课程设计完成的还算顺利,虽然也遇到过一些问题,但通过和同学讨论一起学习都能解决。
当然,我们也都明白matlab的确是一个很实用的工具,在今后的学习中我们会不断的边学边运用它,而且我们还可以将它用在我们专业的学习中。
5参考文献
[1]泉,数字信号处理原理与实现,电子工业出版社,2009
[2]郭仕剑,MATLAB7.x数字信号处理,人民邮电出版社,2007
[3]陈怀琛,MATLAB及在电子信息课程中的应用,电子工业出版社,2006
[4]高会生,MATLAB实用教程(第2版),电子工业出版社,2010
[5]陈怀琛,数字信号处理教程——MATLAB释义与实现,电子工业出版社,2004
附录
辅助函数1
function[Hr,w,a,L]=Hr_Type1(h);
%计算第一种低通滤波器设计的振幅响应Hr(w)
%-----------------------------------------------------------
%[Hr,w,a,L]=Hr_Type1(h)
%Hr=振幅响应
%w=在[0pi]区间内计算Hr的500个频率点
%a=第一种低通滤波器的系数
%L=Hr的阶次
%h=第一种低通滤波器的频率响应
M=length(h);
L=(M-1)/2;
a=[h(L+1)2*h(L:
-1:
1)];
%1乘(L+1)行向量rowvector
n=[0:
1:
L];
%(L+1)乘1列向量
w=[0:
500]'
*pi/500;
Hr=cos(w*n)*a'
辅助函数2
function[db,mag,pha,grd,w]=freqz_m(b,a);
%freqz子程序的改进版本
%------------------------------------
%[db,mag,pha,grd,w]=freqz_m(b,a);
%db=[0到pi弧度]区间内的相对振幅(db)
%mag=[0到pi弧度]区间内的绝对振幅
%pha=[0到pi弧度]区间内的相位响应
%grd=[0到pi弧度]区间内的群迟延
%w=[0到pi弧度]区间内的501个频率样本向量
%b=Ha(z)的分子多项式系数(对FIRb=h)
%a=Ha(z)的分母多项式系数(对FIR:
a=[1])
[H,w]=freqz(b,a,1000,'
whole'
H=(H(1:
501))'
w=(w(1:
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
pha=angle(H);
%pha=unwrap(angle(H));
grd=grpdelay(b,a,w);
%grd=diff(pha);
%grd=[grd
(1)grd];
%grd=[0grd(1:
500);
grd;
grd(2:
501)0];
%grd=median(grd)*500/pi;