基于MATLAB的数字滤波器的设计与仿真.docx
《基于MATLAB的数字滤波器的设计与仿真.docx》由会员分享，可在线阅读，更多相关《基于MATLAB的数字滤波器的设计与仿真.docx（11页珍藏版）》请在冰豆网上搜索。
基于MATLAB的数字滤波器的设计与仿真
一、课题简介
本课题是基于MATLAB的数字滤波器的设计与仿真，采用MATLAB软件设计与仿真。
有限冲击响应数字滤波器（FIR）具有突出的优点:
系统总是稳定的、易于实现线性相位、允许设计多通带（或多阻带）滤波器。
首先在了解有限冲击响应数字滤波器的基本概念和数学模型的前提下，给出有限冲击响应数字滤波器具有线性相位的条件，以及有限冲击响应数字滤波器的各种结构及其特点。
其次，由于在实际工程设计限冲击响应数字滤波的时候，窗函数设计法和频率采样法都存在设计精度不高，运算量大，边缘频率不容易确定的缺点。
而优化设计法恰能弥补上述方法的不足，能很好的逼近理想数字滤波器。
最后，在Simulink环境下建立一个数字滤波器系统仿真模型，用优化设计法和频率采样法分别设计相同指标的滤波器。
把原始信号和干扰信号同时输入，两种方法设计的滤波器分别在仿真模型中滤除干扰。
以仿真图的形式直观的给出滤波器的性能。
二、设计过程
⒈有限长单位冲激响应（FIR）滤波器的基本结构
⑴直接型：
如图1-1可以看出直接型结构共需要N个乘法器，若系数不对称则不能设计线性相位。
图1-1FIR滤波器的直接型结构
⑵级联型：
将H（z）分解成实系数二阶因子的乘积形式
（1.1）
这种结构的每一节控制一对共轭极点，因此调整传输零点方便，但是这种结构所需的系数和所需的乘法运算比直接型多，所以这种结构使用的比较少。
图1-2FIR滤波器的级联型结构
⑶频率抽样型：
把一个有限长序列（长度为N点）的z变换H（z）在单位圆上作N等分抽样，就得到H（k），其主值序列就等于h（n）的离散傅里叶变换H（k）。
用H（k）表示的H（z）的内插公式为
（1.2）
（1.3）
其中
为梳状滤波器，
为谐振器。
谐振器的极点正好与梳状滤波器的零点相抵消，保证了网络的稳定性。
N个并联谐振器与梳状滤波器级联后，得到图1-3的频率抽样结构。
图1-3FIR滤波器的频率抽样型结构
2.FIR数字滤波器的设计方法
2.1窗函数设计法流程图如2-1所示：
图2-1窗函数设计流程
常用的窗函数有：
矩形窗、汉宁窗、海明窗、布莱克曼窗、凯塞窗、三角窗等。
在窗函数中，凯塞窗是比较灵活的一种窗函数，调整凯塞窗中的参数β的大小，我们可以得到不同性能的滤波器。
窗函数设计法是一种常用的设计数字滤波器的方法，不过在设计精度和性能方面并不是十分理想。
几种窗函数的性能比较如下表1所示：
表1各种窗函数的参数表
窗函数
过渡带宽度（P/N）
阻带最小衰减（dB）
旁瓣峰
矩形窗
4
21
13
三角窗
8
25
25
汉宁窗
8
44
31
海明窗
8
53
41
凯塞窗（β=5.6）
7.442
60
51
布莱克曼窗
12
74
71
2.2频率采样设计法流程图如2-2所示：
频率采样法的缺点是在信号的幅值出不容易控制，而且其设计过程不好理解。
图2-2频率采样设计流程
2.3优化设计法流程图如2-3所示：
交替定理为优化法提供了一组很有用的充分必要条件。
这里的优化算法也叫做Remez算法，是一种可以用来设计线性相位的FIR滤波器的高效算法。
图2-3优化设计流程
3.不同设计方法的比较
本文将举例要求用MATLAB实现一个线性相位FIR数字低通滤波器，参数要求：
通带临界频率fp=2000Hz，通带允许波动Rp=ldB，阻带临界频率fs=4000Hz，阻带衰减Rs=50db，截止频率Fc=（fp+fs）／2=3000Hz，采样频率Fs=22050Hz。
3.1用窗口设计法没计的m程序代码如下：
clear;closeall;clc;
fp=2000;fs=4000;Fs=22050;Rp=1;Rs=5;%输入设计指标
W1p=fp/（Fs/2）;Wls=fs/（Fs/2）;W1c=（W1p+Wls）/2;%求归一化频率
%选择窗函数的类型（Rs=50选哈明窗）并估计窗口长度M
dW=Wls-W1p;M=ceil（6.6/dW）;
ifmod（M,2）==0;N=M+I,elseN=M;end;%选用第一类滤波器
n=0:
N-1;r=（N-1）/2;%用于计算理想低通单位脉冲响应
hdn=sin（W1c*pi*（（n-r）+eps））/（pi*（（n-r）+eps））;
win=hamming（N）;hn=hdn.*win';%用哈明窗加窗
%以上三行可用hn=firl（N-1，W1c，hamming（N））直接求得
freqz（hn,1,512,Fs）;gridon;%绘制结果并加网络
程序运行结果如下：
N=37（hn只取了前20个数据，图如3-1所示）
hn=0.00040.00150.00190.0008-0.0024-0.0059
-0.00610.00020.0l130.0l890.0130-0.0093-0.0374
-0.0482-0.01950.05430.15320.23850.27210.2385
滤波器的阶数为37阶（即满足要求的滤波器的最低阶数为37阶），如果把窗函数的阶数降低，则最后设计得到的滤波器性能指标不符合设计要求。
除了矩形窗之外，一般的窗函数能满足实际工程的需求，但是在精度要求比较高的场合，窗函数设计往往达不到要求。
图3-1幅频及相频响应曲线
3.2采用频率采样法设计的m程序代码如下：
clear;closeall;clc;
fc=3000;Fs=22050;W1c=fc/（Fs/2）;%指标归一化
N=37;%选用阶数N为奇的第一类滤波器
%根据约束务件确定H（k）的值
w=[0:
N-1]*2*pi/N;m=fix（W1c/（2/N）+1）;
abs_H=[ones（1,m）,zeros（1,N-2*m+1）,ones（1,m-1）];
angles_H=-w*（N-1）/2;
H=abs_H.*exp（j*angles_H）;
hn=real（ifft（H））;%求单位脉冲响应
freqz（hn,1,512,Fs）;gridon;%绘制结果并加网络
程序运行结果如下（列出hn前20个数据，图如3-2所示）：
hn=-0.0242-0.00460.0l910.02800.0142-0.0125
-0.0310-0.02520.003l0.03330.03990.0122-0.0348
-0.0656-0.04540.03570.15290.25620.29730.2562
频率采样法采用了插值逼近理论，当理想频响曲线越平缓，则内插值越接近理想值，逼近误差越小，故设计时常在不连续的边缘加一些过渡点，以减小通带和阻带的波动。
可见，此方法设计滤波器的阶数比用窗口函数法设计的滤波器阶数少。
图3-2幅频及相频响应曲线
3.3采用优化设计法设计的m程序代码如下：
clear;closeall;clc;
fp=2000;fs=4000;Fs=22050;Rp=1;Rs=50;%设计指标
W1p=fp/（Fs/2）;W1s=fs/（Fs/2）;%求归一化频率
F=[W1p,W1s];A=[1,0];%设置边界频率和幅度要求
DEVp=（10^（Rp/20）-1）/（10^（Rp/20）+1）;DEVs=10^（-Rs/20）;
DEV=[DEVp,DEVs];%各频带内的波纹要求
[N,Fo,Ao,W]=remezord（F,A,DEV）;%确定remez参数
hn=remez（N,Fo,Ao,W）;%调用remez函数进行设计
freqz（hn,1,512,Fs）;gridon;%绘制结果并加网络
程序运行结果如下：
N=18（列出hn前11个数据，图如3-3所示）
hn=-0.0005-0.0122-0.0272-0.0400-0.0347
0.00160.06970.15290.22180.24850.2218
如果设计的滤波器指标不合要求，可将remezord输出的参数适当修改。
一般说来，改变remez的阶数，通带阻带指标都会更好；减小通带加权，可能增加阻带衰耗，但通带波动增加；加大阻带加权，将增加阻带衰减。
从理论和数值上分析，相比窗函数法和频率采样法，采用优化设计法设计的滤波器系数较少。
图3-3幅频及相频响应曲线
4.数字滤波器的设计与分析工具（FDAToo1）
工具设计法是利用MATLAB提供的滤波器设计与分析工具（FDAToo1）进行设计的一种方法。
该窗口分为上下两部分：
上面是设计结果显示；下面用来设定所需的技术参数。
FDATool需设置的参数主要有响应类型、设计方法、滤波器阶数及选项、频率参数和幅度参数等项目。
不同类型和不同方法的滤波器设计参数不尽相同，图4-1给出的是前述实例的设计指标。
滤波器的设计工作完成后，可以借助于Export操作导出所设计滤波器的系统函数H（z）。
图4-1利用FDAtool的滤波器设计窗口
三、结果分析
利用MATLAB设计FIR数字滤波器有多种方法。
窗口设计法概念清楚，但临界频率难以控制，而且当理想频响是任意曲线时，求hd（n）也比较困难；频率采样法在取样点处精确保证频响要求，克服了窗口设计法转折频率不易控制缺点，但不能确保截止频率自由取值；优化设计法在相同阶数下，可以获得更好的频率特性和衰耗特性，但通带不平滑；FDATool工具设计法操作简单、
设计方便，但不好理解。
无论哪种方法，都要注意选择类型，以便确定阶数，如果不满足指标，应逐渐增加或减小阶数直到符合要求。
四、数字滤波器的应用-软件实现（滤噪）
程序流程框图如图5-1所示：
图5-1程序框图
functionxt=xtg
%xt=xtg产生一个长度为N，有加性高频噪声的单频调幅信号xt，N=1000
N=1000;Fs=1000;T=1/Fs;Tp=N*T;%采样频率Fs=1000Hz
t=0:
T:
（N-1）*T;
fc=Fs/10;f0=fc/10;%载波频率fc=Fs/10=100Hz，调制正弦波f0=fc/10=10Hz
mt=cos（2*pi*f0*t）;%产生单频调制信号mt，频率为f0
ct=cos（2*pi*fc*t）;%产生载波正弦波信号ct，频率为fc
xt=mt.*ct;%相乘产生单频调幅信号xt
nt=2*rand（1,N）-1;%产生随机噪声nt
%设计高通滤波器hn，用于滤除噪声nt中的低频成分，生成高通噪声
fp=150;fs=200;Rp=0.1;As=70;%滤波器指标
fb=[fp,fs];m=[0,1];%计算remezord函数所需参数f，m，dev
dev=[10^（-As/20）,（10^（Rp/20）-1）/（10^（Rp/20）+1）];
[n,fo,mo,W]=remezord（fb,m,dev,Fs）;%确定remez函数所需参数
hn=remez（n,fo,mo,W）;%调用remez函数进行设计，用于滤除噪声nt中的低频成分
yt=filter（hn,1,10*nt）;%滤除随机噪声中的低频成分，产生高通噪声yt
xt=xt+yt;%噪声加信号
fst=fft（xt,N）;k=0:
N-1;f=k/Tp;
subplot（2,2,1）;plot（t,nt）;grid;xlabel（'t/s'）;ylabel（'n（t）'）;
axis（[0,Tp/5,min（nt）,max（nt）]）;title（'（a）滤波前的波形'）;
subplot（2,2,2）;plot（t,yt）;grid;xlabel（'t/s'）;ylabel（'y（t）'）;
axis（[0,Tp/5,min（yt）,max（yt）]）;title（'（b）滤波后的波形'）;
subplot（2,2,3）;plot（t,xt）;grid;xlabel（'t/s'）;ylabel（'x（t）'）;
axis（[0,Tp/5,min（xt）,max（xt）]）;title（'（c）信号加噪声波形'）;
subplot（2,2,4）;plot（f,abs（fst）/max（abs（fst）））;grid;title（'（d）信号加噪声频谱'）;
axis（[0,Fs/2,0,1.2]）;xlabel（'f/Hz'）;ylabel（'幅度'）;
运行结果如图6-1所示：
图6-1滤波前后波形显示
五、进一步需要解决的问题
多阅读相关文献资料，进一步完善软件部分的设计，通过观察运行出的图形及数据，能更加直观的认识滤波器的性能。
下一步工作的重点是在Simulink环境下如何修改滤波器的仿真参数，从而很好的实现对滤波器的可视化建摸与仿真。
六、后期工作安排
继续研究课题，使设计更加完善。
对该课题进行分析整理，完成毕业论文，准备答辩。
升级会员 