基于matlab的低通滤波器.docx
《基于matlab的低通滤波器.docx》由会员分享，可在线阅读，更多相关《基于matlab的低通滤波器.docx（17页珍藏版）》请在冰点文库上搜索。
基于matlab的低通滤波器
基于matlab的低通滤波器
摘要:
调用MATLAB信号处理工具箱中滤波通过观察滤波器输入输出信号的时域波形及其频谱，建立数字滤波的概念。
应用最广泛的是双线性变换法。
基本设计过程是：
①先将给定的数字滤波器的指标转换成过渡模拟滤波器的指标；②设计过渡模拟滤波器；③将过渡模拟滤波器系统函数转换成数字滤波器的系统函数。
MATLAB信号处理工具箱中的各种IIR数字滤波器设计函数都是采用双线性变换法。
关键词:
滤波器，matlab，c语言，声音
Abstract:
callMATLABsignalprocessingtoolboxfilteringthroughtheobservationfilterofinputandoutputsignalstimedomainwaveformandspectrum,establishtheconceptofdigitalfilter.Oneofthemostwidelyappliedisdoublelineartransformationmethod.Thebasicdesignprocessis:
（1）toagivendigitalfiltersindexconvertedintoanalogfiltertransitionindex;
（2）thedesigntransitionsimulationfilter;（3）transitionsimulationfiltersystemfunctionwillbeconvertedintodigitalfiltersystemfunction.MATLABsignalprocessingtoolboxdigitalfilterfunctiondesignofIIRisbilineartransformationmethod.
Keywords:
filter,matlab,theclanguage,thevoice
一.任务：
用matlab软件设计IIR模拟、数字以及各种窗函数的FIR低通滤波器
二.设计目的：
（1）了解matlab软件的用途以及用法；
（2）了解用冲激响应不变法设计模拟低通滤波器；
（3）了解用脉冲响应不变法设计的巴特沃思数字低通滤波器；
（4）了解基于汉明窗函数的FIR低通滤波器的设计。
三.设计内容：
3.1用冲激响应不变法设计模拟低通滤波器
3.1.1设计内容：
要求按照设计指标设计无限冲激响应IIR巴特沃什模拟低通滤波器。
3.1.2设计原理：
低通滤波器的技术要求用图形表示如下：
通带过渡带阻带
如上图表示了一个频域设计、一维低通滤波器的技术要求图。
ws和wp分别称为通带截止频率和阻带截止频率。
通带频率范围为0≤w≤wp，阻带频率范围ws≤w≤π。
从wp到ws称为过渡带，用Δw表示，在过渡带里，幅频特性单调下降。
在通带和阻带内的衰减一般用数dB表示，通带内允许最大衰减是αs，阻带内允许最小衰减是αp。
定义分别为：
Matlab信号处理工具箱中提供了设计巴特沃思模拟滤波器的函数buttord、buttap和butter，格式如下：
用于计算巴特沃思模拟低通滤波器的阶N和3dB截止频率Wc（即本书中的符号
）。
其中，Wp和Ws分别是滤波器的通带截止频率
和阻止截止频率
，单位为rad/s；Rp和Rs分别是通带最大衰减系数
和阻带最小衰减系数
，单位为dB。
用于计算N阶巴特沃思归一化（
=1）模拟低通滤波器系统函数的零、极点和增益因子，返回长度为N的向量z和p分别给出N个零点和极点，G是滤波器增益。
得到的滤波器系统函数形式如下：
其中，
和
分别是向量z和p的第k个元素。
如果要从零、极点得到系统函数的分子和分母多项式系数向量B和A，可以调用结构转换函数：
结构转换后系统函数的形式为：
其中，M是向量B的长度，N是向量A的长度，
分别是向量B和A的第k个元素。
用于计算巴特沃什模拟滤波器系统函数中分子和分母多项式系数向量B和A，其中N和
分别是滤波器的阶和3dB截止频率
，返回向量B和A中的元素
和
分别是上面的
表示式中的分母和分子系数。
ftype缺省时，设计低通滤波器，S缺省时，设计数字滤波器。
3.1.3设计程序如下：
（1）巴特沃思模拟低通滤波器设计指标如下：
wp=0.4*pi,ws=0.45*pi,Rp=5,Rs=10;
实验所用到的程序如下
wp=0.4*pi;ws=0.45*pi;Rp=5;Rs=10;
[N,wc]=buttord（wp,ws,Rp,Rs,'s'）;
[B,A]=butter（N,wc,'s'）;
w=0:
0.001:
6;
[H,W]=freqs（B,A,w）;
H=20*log10（abs（H））;
plot（w,H）,gridon;
xlabel（‘频率/（rad/s）'）;ylabel（'幅度/dB'）
经matlab仿真以后得到的滤波器如下图所示：
图1.1巴特沃思模拟低通滤波器
运行后在matlab控制窗口上显示如下结果：
N=7
（2）为了方便观察故设计了不同阶次N的滤波器，即Butterworth模拟原型低通滤波器，其程序代码如下：
n=0:
0.01:
2;
forii=1:
4
switchii
case1,N=2;
case2,N=5;
case3,N=10;
case4,N=20;
end[z,p,k]=buttap（N）;
[b,a]=zp2tf（z,p,k）;
[H,w]=freqs（b,a,n）;
magH2=（abs（H））.^2;
holdon;
plot（w,magH2）;
end
xlabel（'w/wc'）;
ylabel（'|H（jw）|^2'）;
title（'Butterworth模拟原型低通滤波器'）;
text（1.5,0.18,'n=2'）
text（1.3,0.08,'n=5'）
text（1.16,0.08,'n=10'）
text（0.93,0.98,'n=20'）
gridon;
仿真后结果图形如图：
图1.2Butterworth模拟原型低通滤波器
（3）利用冲激响应不变法，把系统函数为
的模拟滤波器变换成等价的数字滤波器，采样间隔T=1S。
先经过计算得到：
，则有如下式子
Matlab实验参考程序如下:
B=[1,1];
A=[1,5,6];
T=1;
Fs=1/T;
[Bz,Az]=impinvar（B,A,Fs）;%用冲激响应不变法将模拟滤波器变换成数字滤波器
运行结果：
Bz=1.0000-0.2209
Az=1.0000-0.18510.0067
得到的数据和先前计算得到的数据相同，故利用冲激响应不变化法把模拟滤波器设计成数字滤波器正确。
3.1.4实验结果及分析
（1）步骤1中产生的低通模拟滤波器的系统函数，经计算为：
经过仿真后的结果，与上式对比，所得结果进分析相同。
此处为N=7的低通模拟滤波器
（2）在Butterworth模拟原型低通滤波器中，可以看到，滤波器的幅频特性随着滤波器阶次N的增加而变得越来越好。
（3）步骤3中，利用冲激响应不变法，把系统函数的模拟滤波器变换成等价的数字滤波器，在matlab仿真软件中验证了冲激响应不变法的正确性。
3.2用脉冲响应不变法设计的巴特沃思数字低通滤波器
3.2.1设计内容：
设计巴特沃思数字低通滤波器和椭圆数字低通滤波器，要求通带边界频率fp=2.1kHZ，通带最大衰减Rp=0.5dB；阻带边界频率fs=8kHZ，阻带最小衰减Rs=30dB,采样频率为Fs=20kHZ。
3.2.2设计步骤：
用MATLAB进行典型的数字滤波器的设计，一般步骤如下:
（1）将设计指标归一化处理。
如果采用双线性变换法，还需进行预畸变。
（2）根据归一化频率，确定最小阶数N和频率参数Wn。
可供选用的阶数择函数有:
buttord，cheblord，cheb2ord，ellipord等。
（3）运用最小阶数N设计模拟低通滤波器原型。
模拟低通滤波器的创建函数有:
buttap,cheblap,cheb2ap,ellipap和besselap，这些函数输出的是零极点式形式，还要用zp2tf函数转换成分子分母多项式形式。
如果想根据最小阶数直接设计模拟低通滤波器原型，可用butter,chebyl,cheby2,ellip,bessel等函数，只是注意要将函数中的Wn设为1。
（4）根据第2步的频率参数Wn，模拟低通滤波原型转换模拟低通、高通、带通、带阻滤波器，可用函数分别是：
lp21p,lp2hp,lp2bp,lp2bs。
（5）运用脉冲响应不变法或双线性变法把模滤波器转数字滤波器，调用的函数是impinvar和bilinear。
脉冲响应不变法适用于采样频率大于4倍截止频率的锐截止低通带通滤波器，而双线性变换法适合于相位特性要求不高的各型滤波器。
（6）根据输出的分子分母系数，用tf函数生成H（z）的表达式，再用freqz函数验证设计结果。
3.2.3设计程序如下：
fp=2100;
fs=8000;
Fs=20000;
Rp=0.5;
Rs=30;
T=1/Fs;%设计指标
W1p=fp/Fs*2;W1s=fs/Fs*2;%求归一化频率
[N,Wn]=buttord（W1p,W1s,Rp,Rs,'s'）;%确定butterworth的最小介数N和频率参数Wn
[z,p,k]=buttap（N）;%设计模拟低通原型的零极点增益参数
[bp,ap]=zp2tf（z,p,k）;%将零极点增益转换成分子分母参数
[bs,as]=lp2lp（bp,ap,Wn*pi*Fs）;%将低通原型转换为模拟低通
[bz,az]=impinvar（bs,as,Fs）;%用脉冲响应不变法进行模数变换
sys=tf（bz,az,T）;%给出传输函数H（Z）
[H,W]=freqz（bz,az,512,Fs）;%生成频率响应参数
subplot（2,1,1）;
plot（W,20*log10（abs（H）））;%绘制幅频响应
gridon;%加坐标网格
xlabel（'频率/Hz'）;
ylabel（'振幅/dB'）;
subplot（2,1,2）;
plot（W,abs（H））;gridon;
xlabel（'频率/Hz'）;
ylabel（'振幅/H'）;
运行后的波形如下：
图2.1巴特沃思和椭圆数字低通滤波器
N=4
bz=0.00000.09990.19140.0252
az=1.0000-1.43361.0984-0.41150.0627
3.3基于汉明窗函数的FIR低通滤波器的设计
3.3.1设计具体低通滤波器要求：
设低通滤波器通带截止频率为
，阻带截止频率为
，通带纹波
=0.5dB，最小阻带衰减
=50dB。
3.3.2设计内容及原理：
对于线性相位FIR滤波器通常采用窗函数法设计。
窗函数设计FIR滤波器基本思想，就是从时域出发,把理想的无限长的
用一定形状的窗函数截取成有限长的
，以此
来逼近
，从而使所得到的频率响应
与所要求的理想频率响应
相接近。
理想的数字滤波器频率特性
是无法实现的，FIR的设计就是要寻找一个可以得到的频率特性
来逼近
，这相当于用一个可实现的单位脉冲响应
去逼近一个理想单位脉冲响
。
要想用一个有限长的因果序列去逼近它，最简单的方法是截取n从0~N-1的一段来表示它，即：
（0≦n≦N-1）；当n为其他的数时，
=0。
同时，为了保证线性相位，还要满足偶对称：
，这就好像通过一个窗口观看到的一段
，因此
就表示成
和一个“窗口函数”的乘积。
这样，对
的求解就变为：
这里的
就称为窗口函数。
窗口函数相当于对
在时域上截取一段，势必造成频域的矩形窗口的失真，结果就是截取出的信号也相应失真。
为了补偿这种失真，只有改变原来窗口的形状，修正经过时域截取后的窗口失真。
改善窗口函数的标准：
（1）尽量减少窗口频谱的旁瓣，使能量集中在主瓣中，以减少波动的幅度和次数，提高阻带衰减。
这样，主瓣就不能太窄；
（2）主瓣的宽度应尽量窄，以得到尽量窄的过渡带。
工程中常用的窗函数有以下六种：
矩形窗，三角形窗，汉宁窗，汉明窗，布莱克曼窗，凯泽窗（可变窗函数）。
表1为六种窗函数的性能指标。
窗函数窗普性能指标加窗后滤波器性能指标
----------------------------------------------------------------------------------
旁瓣峰值dB主瓣宽度2π/n过滤宽度2π/n阻带最小衰减
矩形窗-1320.9-21
三角窗-2542.1-25
汉宁窗-3143.1-44
汉明窗-4143.3-53
布莱克曼窗-5765.5-74
凯泽窗-5755-80
表3.1常用传函数的性能指标
3.3.3根据条件采用汉明窗设计程序如下：
wp=0.2*pi;wr=0.4*pi;
trwidth=wr-wp; %过渡带宽度
N=ceil（6.64*pi/trwidth）+1;%滤波器的长度
wc=（wr+wp）/2;
hd=ideal_lp（wc,N）;
w_ham=（hamming（N））';
h=hd.*w_ham;
[db,mag,pha,grd,w]=freqz_m（h,[1]）;
delta_w=2*pi/1000;
Ap=-（min（db（1:
1:
wp/delta_w+1）））;
Ar=-round（max（db（wr/delta_w+1:
1:
501）））;
subplot（221）;stem（n,hd）;title（'理想单位脉冲响应hd（n）'）
subplot（222）;plot（n,w_ham）;title（'海明窗w（n）'）
subplot（223）;stem（n,h）;title（'实际单位脉冲响应h（n）'）
subplot（224）;plot（w/pi,db）;title（'幅度响应（dB）'）
输出结果如下：
图3.1
由上图可以看到：
通带纹波为0.0301dB，最小阻带衰减为55dB。
满足设计指标的要求。
滤波器系数N=35具体系数如下：
-0.0004628854790980980.001027078059784260.002356825273471590.00198992091909062-0.00151370239796531-0.00662950110914604-0.007841668819955144.84287457256860e-180.0142366469196944.022*********0.00935693617842729.023*********-0.0520253274593482.0411********0.03052864357449910.1466635181306000.2555011292338670.300000000000000
下面我把滤波器的阶数人为的设置成15，我们得到如下的输出结果：
图3.2
由上图可以看到：
通带纹波为1.6322dB，最小阻带衰减为16dB。
显然不能满足指标要求。
滤波器系数N=15，具体系数如下：
0.00112415044923825-0.00391515738696378-0.0161188746641399.020*********.021*********0.1251496714246940.2457870449329650.300000000000000
综上所述：
当滤波器的阶数增加时，可以减小通带纹波，增加最小阻带衰减。
但同时计算量也随之增大，成本也要同样增加。
下面，我用MATLAB中的FilterDesignandAnalysis工具进行同样的设计，设置阶数为35，低通，窗口选择hamming，归一化频率
。
结果如下：
图3.3幅频特性曲线
图3.4相频特性曲线
图3.5冲击响应
图3.6零极点图
3.3.4实验结论及分析
FIR滤波器以它优越的性能，在数字信号处理领域中占有很重要的地位，是现代电子技术中必须掌握的设计技术。
本文通过讨论FIR数字滤波器的基本原理和设计步骤，阐述了窗函数法设计FIR数字滤波器的思路。
在MATLAB中,通过窗函数法设计了FIR数字滤波器，并很好的完成了设计要求的指标。
参考文献
[1]张辉，曹丽娜，王勇，通信原理辅导.，西安：
西安电子科技大学出版社，1992.
[2]王兴亮，等.数字通信原理与技术.西安：
西安电子科技大学出版社，2001.
[3]宋祖顺，等.现代通信原理.北京：
电子工业出版社，2001.
[4]赵树杰等，数字信号处理，西安电子科技大学出版社，1997.10.
[5]陈怀琛等，MATLAB及在电子信息课程中的应用，电子工业出版社出版，2002.
[6]李丽王振领，MATLAB工程计算及应用，人民邮电出版社，2001.
[7]陈怀琛等译，《数字信号处理及其MATLAB实现》，电子工业出版社；2002.
升级会员 