课程设计IIR数字低通滤波器的MATLAB实现Word下载.docx
《课程设计IIR数字低通滤波器的MATLAB实现Word下载.docx》由会员分享，可在线阅读，更多相关《课程设计IIR数字低通滤波器的MATLAB实现Word下载.docx（21页珍藏版）》请在冰豆网上搜索。
设IIR滤波器的输入序列为x（n），则IIR滤波器的输入序列x（n）与输出序列y（n）之间的关系可以用下面的方程式表示：
其中，
和
是滤波器的系数，其中
中至少有一个非零。
与之相对应的差分方程为：
由传递函数可以发现无限常单位冲激响应滤波器有如下特点：
（1）单位冲激响应h（n）是无限长的。
（2）系统传递函数H（z）在有限z平面上有极点存在。
（3）结构上存在着输出到输入的反馈，也就是结构上是递归型的。
（二）IIR数字滤波器的设计与实现
IIR数字滤波器的设计有多种方法，如频率变换法、数字域直接设计以及计算辅助设计等。
下面只介绍频率变换设计法。
首先考虑由模拟低通滤波器到数字低通滤波器的转换，其基本的设计过程如下:
（1）将数字滤波器的技术指标转换为模拟滤波器的技术指标；
（2）设计模拟滤波器G（S）；
（3）将G（S）转换成数字滤波器H（Z）；
在低通滤波器的设计基础上，可以得到数字高通、带通、带阻滤波器的设计流程如下：
（1）给定数字滤波器的设计要求（高通、带阻、带通）；
（2）转换为模拟（高通、带阻、带通）滤波器的技术指标；
（3）转换为模拟低通滤波器的指标；
（4）设计得到满足第三步要求的低通滤波器传递函数；
（5）通过频率转换得到模拟（高通、带阻、带通）滤波器；
（6）变换为数字（高通、带阻、带通）滤波器。
在matlab中设计IIR滤波器的方法及其它们所用到的函数如表6－1所示。
表6－1matlab中设计IIR滤波器的方法列表
方法
描述
函数
模拟原型法
采用经典低通滤波器作为连续域上的设计模型，通过频率变换得到IIR数字滤波器，最后进行离散化处理
完整设计函数：
Beself,butter,cheby1,cheby2,ellip
滤波器的阶估计函数：
Buttord,cheb1ord,cheb2ord,ellipord
低通模拟滤波器原型函数：
beselap,buttap,cheb1ap,cheb2ap,ellipap
频域变换函数：
Lp2bp,lp2bs,lp2hp,lp2lp
其他函数：
Bilinear,impinvar
直接设计方法
直接在离散时域上估计线性的幅度响应
yulewalk
通用butterworth
设计方法
使用butterworth设计低通数字滤波器
Maxflat
参数建模方法
寻找接近于所需要设计的滤波器的通用原型
时域上的建模函数：
Lpc,prony,stmcb
频域上的建模函数：
Invfreqs,invfreqz
（三）双线性变换法
与冲激响应不变法比较，双线性变换的主要优点是靠频率的非线性关系得到S平面与Z平面的单值一一对应关系，整个值对应于单位圆一周。
所以从模拟传递函数可直接通过代数置换得到数字滤波器的传递函数。
【实例6-4】设计一个截止频率为200Hz的数字低通滤波器，采用频率为1000Hz。
解：
源代码如下：
[z,p,k]=buttap（3）;
[b,a]=zp2tf（z,p,k）;
[bt,at]=lp2lp（b,a,200*2*pi）;
[bz,az]=bilinear（bt,at,1000）;
freqz（bz,az,512,1000）
结果如图6－4所示：
bz=
0.07530.22590.22590.0753
az=
1.0000-0.82660.5154-0.0865
【实例6-5】基于Butterworth模拟滤波器原型，使用双线性转换设计数字滤波器，其中参数指标为：
通带截止频率：
通带波动值：
阻带截止频率：
阻带波动值：
首先确定滤波器的阶数N，同时根据
确定
＝0.5。
接着使用bilinear进行双线性转换，最后绘制在频域上的各种图像，其源代码如下：
％数字滤波器指标
wp=0.2*pi;
ws=0.3*pi;
Rp=1;
As=15;
％将数字滤波器指标反转变化为模拟滤波器的参数
T=1;
fs=1/T;
omegap=（2/T）*tan（wp/T）;
omegas=（2/T）*tan（ws/T）;
ep=sqrt（10^（Rp/10）-1）;
Ripple=sqrt（1/（1+ep*ep））;
Attn=1/（10^（As/20））;
％butterworth原型模拟滤波器的设计
[cs,ds]=afd_butt（omegap,omegas,Rp,As）;
％双线性变换
[b,a]=bilinear（cs,ds,T,fs）;
％频域图像的绘制
freqz（b,a）;
程序运行后，产生4阶的butterworth数字滤波器，频率响应如图6－5所示的波形。
图6－54阶Butterworth数字滤波器
3、实验步骤
（1）录入一段语音信号，设置采样率和采样点数。
（2）输入信号的频谱分析
要求首先画出输入信号的时域波形；
然后对该信号进行快速傅立叶变换，得到信号的频谱特性。
（3）设计数字低通滤波器并画出数字低通滤波器的频率响应。
（4）用滤波器对信号进行滤波，比较滤波前后语音信号的波形及频谱。
（5）利用MatlabGUI工具做一个图形用户界面，窗口简洁大方。
4、实验结果及截图
结果
5．总结
通过这次课程设计，使自己对通信原理中的信号的变换，特别是傅立叶变换及其实现过程有了更深的了解。
我对matlab软件也有了熟练的了解，对工程项目设计积累的一定经验。
起初在面对一个陌生时会有忌惮和不知所措，而通过这次训练基本可以通过help指令和doc指令来自学函数用法。
在程序调试的时候经常会出现错误，耐心调试找到错误并搞清楚为什么会出错则是作为一个程序员的一项基本素质。
我们团队分工明确，成员们勤劳能干，奠定我们高效地完成课程设计的基础。
而这期间，我们也确实遇到了很多困难，但我们都攻克了。
总之，我们在这次课程设计中学到了很多。
参考文献：
【1】通信原理第六版樊昌信曹丽娜编著国防工业出版社
【2】信号与系统引论郑君里应启杨为理编著高等教育出版社
【3】MATLAB及在电子信息类课程中的应用第二版唐向宏岳恒立郑雪峰编著电子工业出版社
附录：
源程序
functionvarargout=idea2（varargin）
%IDEA2M-fileforidea2.fig
%IDEA2,byitself,createsanewIDEA2orraisestheexisting
%singleton*.
%
%H=IDEA2returnsthehandletoanewIDEA2orthehandleto
%theexistingsingleton*.
%IDEA2（'
CALLBACK'
hObject,eventData,handles,...）callsthelocal
%functionnamedCALLBACKinIDEA2.Mwiththegiveninputarguments.
Property'
'
Value'
...）createsanewIDEA2orraisesthe
%existingsingleton*.Startingfromtheleft,propertyvaluepairsare
%appliedtotheGUIbeforeidea2_OpeningFunctiongetscalled.An
%unrecognizedpropertynameorinvalidvaluemakespropertyapplication
%stop.Allinputsarepassedtoidea2_OpeningFcnviavarargin.
%*SeeGUIOptionsonGUIDE'
sToolsmenu.Choose"
GUIallowsonlyone
%instancetorun（singleton）"
.
%Seealso:
GUIDE,GUIDATA,GUIHANDLES
%Copyright2002-2003TheMathWorks,Inc.
%Edittheabovetexttomodifytheresponsetohelpidea2
%LastModifiedbyGUIDEv2.513-Jun-201216:
27:
25
%Begininitializationcode-DONOTEDIT
gui_Singleton=1;
gui_State=struct（'
gui_Name'
mfilename,...
'
gui_Singleton'
gui_Singleton,...
gui_OpeningFcn'
@idea2_OpeningFcn,...
gui_OutputFcn'
@idea2_OutputFcn,...
gui_LayoutFcn'
[],...
gui_Callback'
[]）;
ifnargin&
&
ischar（varargin{1}）
gui_State.gui_Callback=str2func（varargin{1}）;
end
ifnargout
[varargout{1:
nargout}]=gui_mainfcn（gui_State,varargin{:
}）;
else
gui_mainfcn（gui_State,varargin{:
%Endinitializationcode-DONOTEDIT
%---Executesjustbeforeidea2ismadevisible.
functionidea2_OpeningFcn（hObject,eventdata,handles,varargin）
%Thisfunctionhasnooutputargs,seeOutputFcn.
%hObjecthandletofigure
%eventdatareserved-tobedefinedinafutureversionofMATLAB
%handlesstructurewithhandlesanduserdata（seeGUIDATA）
%varargincommandlineargumentstoidea2（seeVARARGIN）
%Choosedefaultcommandlineoutputforidea2
handles.output=hObject;
%Updatehandlesstructure
guidata（hObject,handles）;
globalSampleRate
SampleRate=44100;
globalSecLength;
SecLength=60;
%因为wav文件是十分钟的，不能全部读入，所以将其分为10段，一段60秒
globalSecNum;
SecNum=600/SecLength;
globalCurrentSection
%UIWAITmakesidea2waitforuserresponse（seeUIRESUME）
%uiwait（handles.figure1）;
%---Outputsfromthisfunctionarereturnedtothecommandline.
functionvarargout=idea2_OutputFcn（hObject,eventdata,handles）
%varargoutcellarrayforreturningoutputargs（seeVARARGOUT）;
%Getdefaultcommandlineoutputfromhandlesstructure
varargout{1}=handles.output;
%---Executesonbuttonpressinpushbutton1.
functionInsert_Callback（hObject,eventdata,handles）
%hObjecthandletopushbutton1（seeGCBO）
globalCurrentSection;
globalSampleRate;
globalMyrate;
globalSoundData;
globalMycode;
SampleRate=Myrate;
%Sam=Myrate;
%disp（Sam）;
[FileName,PathName]=uigetfile（'
*.wav'
）;
data=wavread（strcat（PathName,FileName）,[SecLength*（CurrentSection-1）*（SampleRate）+1,SecLength*CurrentSection*SampleRate]）;
SoundData=data（1:
10:
length（data）,1）;
%N=Mycode;
n=0:
length（SoundData）-1;
t=n/SampleRate;
figure,plot（t,SoundData）;
%CurrentSection=CurrentSection+1
%---Executesonbuttonpressinpushbutton2.
functionSure_Callback（hObject,eventdata,handles）
%hObjecthandletopushbutton2（seeGCBO）
globalMyrate;
Mycode=str2num（get（handles.edit1,'
string'
））;
Myrate=str2num（get（findobj（'
tag'
edit2'
）,'
%disp（Myrate）;
%---Executesonbuttonpressinpushbutton3.
functionDFT_Callback（hObject,eventdata,handles）
%hObjecthandletopushbutton3（seeGCBO）
globaly;
N=Mycode;
%disp（N）;
y=fft（SoundData,N）;
%进行fft变换
mag=abs（y）;
%求幅值
n=1:
N;
f=n*SampleRate/length（y）;
%进行对应的频率转换
figure;
plot（f,mag）;
%做频谱图
xlabel（'
频率（Hz）'
ylabel（'
幅值'
title（'
频谱'
grid;
%Py=2*（y.*conj（y））/N;
%计算功率谱密度Py
%figure;
%plot（f,Py）;
%xlabel（'
%ylabel（'
功率谱密度'
%title（'
x=120.0*cos（2*3.14*sf*t）功率谱密度'
%grid;
%point=get（handles.edit1,'
%fs=SampleRate;
%m=length（SoundData）;
%Windowlength
%n=pow2（nextpow2（m））;
%Transformlength
%y=fft（SoundData,N）;
%DFT%做n点的FFT，n为2的多少次幂，次幂值为大于或等于length（SoundData）的绝对值的2的最小次幂
%n=（0:
N-1）%Frequencyrange%图上显示的最高频率为采样频率
%power=y.*conj（y）/n;
%PoweroftheDFT%计算复共轭值
%f=n*SampleRate/2;
%figure,plot（f,y）;
Frequency（Hz）'
DFT'
{\bf傅里叶变换}'
%pause
（1）;
%---Executesonbuttonpressinpushbutton4.
functionDTpass_Callback（hObject,eventdata,handles）
%hObjecthandletopushbutton4（seeGCBO）
%%根据原始信号与噪声信号频谱图，选择巴特沃斯低通滤波器
%%滤波器性能指标：
取wp=1500*2*pirad/s,ws=2400*2*pirad/s,Rp=2,Rs=15
%选择滤波器的最小阶数
wp=1000*2*pi;
ws=6000*2*pi;
Rp=2;
Rs=15;
[N1,Wn]=buttord（wp,ws,Rp,Rs,'
s'
%创建butterworth模拟滤波器
[Z,P,K]=buttap（N1）;
%把滤波器的零极点模型转化为传递函数模型
[Bap,Aap]=zp2tf（Z,P,K）;
%把模拟滤波器原型转换成截止频率为Wn的低通滤波器
[b,a]=lp2lp（Bap,Aap,Wn）;
%用双线性变换法实现模拟滤波器到数字滤波器的转换
[bz,az]=bilinear（b,a,SampleRate）;
%绘制频率响应曲线
[H,W]=freqz（bz,az）;
figure（3）;
subplot（211）,
plot（W*SampleRate/（2*pi）,abs（H））;
grid
butterworth低通滤波器频响图'
频率/Hz'
ylabel（'
频率响应幅度'
%%%%%%%用所设计的滤波器对语音+噪声信号y滤波%%%%%%%
y1=filter（bz,az,y）;
n=（0:
Mycode-1）*SampleRate/length（y1）;
subplot（212）,plot（n,y1）;
滤波后信号波形图'
%sound（x,fs）;
%播放语音信号
%sound（y,fs）;
%播放语音+噪声信号
%sound（y1,fs）;
%播放滤波后的信号
%---Executesonbuttonpressinpushbutton5.
functionquit_Callback（hObject,eventdata,handles）
%hObjecthandletopushbutton5（seeGCBO）
close;
msgbox（'
谢谢使用，再见！
'
functionedit1_Callback（hObject,eventdata,handles）
%hObjecthandletoedit1（seeGCBO）
%eventdatareserved-tobedefinedin
升级会员 