課程設(shè)計---iir數(shù)字低通濾波器的matlab實現(xiàn)

1、<p>　　IIR數(shù)字低通濾波器的MATLAB實現(xiàn)</p><p><b>　　目錄</b></p><p><b>　　一、引言4</b></p><p><b>　?。ㄒ唬┰O(shè)計目的4</b></p><p><b>　?。ǘ┰O(shè)計思想4</b

2、></p><p>　　二、IIR數(shù)字濾波器的MATLAB實現(xiàn)</p><p><b>　　三、實驗步驟13</b></p><p>　　四、實驗結(jié)果及截圖14</p><p><b>　　結(jié)果14</b></p><p><b>　　參考文獻15&l

3、t;/b></p><p><b>　　附錄16</b></p><p><b>　　一、引言</b></p><p><b>　　設(shè)計目的</b></p><p>　　（1）熟悉Matlab在信號處理方面的應用</p><p>　?。?）深刻了

4、解頻域上的頻譜特性。（3）了解數(shù)字濾波器的設(shè)計。（4）學會比較濾波前后語音信號的差異。（5）利用GUI工具做一個簡潔大方的圖形用戶界面。</p><p>　　二、IIR數(shù)字低通濾波器的MATLAB實現(xiàn)</p><p>　　（一）IIR數(shù)字濾波器的傳遞函數(shù)及特點</p><p>　　設(shè)IIR濾波器的輸入序列為x(n)，則IIR濾波器的輸入序列x(n)與輸出序列y

5、(n)之間的關(guān)系可以用下面的方程式表示：</p><p>　　其中，和是濾波器的系數(shù)，其中中至少有一個非零。與之相對應的差分方程為：</p><p>　　由傳遞函數(shù)可以發(fā)現(xiàn)無限常單位沖激響應濾波器有如下特點：</p><p>　　單位沖激響應h(n)是無限長的。 </p><p>　　系統(tǒng)傳遞函數(shù)H(z)在有限z平面上有極點存在。</p

6、><p>　　結(jié)構(gòu)上存在著輸出到輸入的反饋，也就是結(jié)構(gòu)上是遞歸型的。</p><p>　　（二）IIR數(shù)字濾波器的設(shè)計與實現(xiàn)</p><p>　　IIR數(shù)字濾波器的設(shè)計有多種方法，如頻率變換法、數(shù)字域直接設(shè)計以及計算輔助設(shè)計等。下面只介紹頻率變換設(shè)計法。首先考慮由模擬低通濾波器到數(shù)字低通濾波器的轉(zhuǎn)換，其基本的設(shè)計過程如下:</p><p>　　將

7、數(shù)字濾波器的技術(shù)指標轉(zhuǎn)換為模擬濾波器的技術(shù)指標；</p><p>　　設(shè)計模擬濾波器G(S)；</p><p>　　將G(S)轉(zhuǎn)換成數(shù)字濾波器H(Z)；</p><p>　　在低通濾波器的設(shè)計基礎(chǔ)上，可以得到數(shù)字高通、帶通、帶阻濾波器的設(shè)計流程如下：</p><p>　　給定數(shù)字濾波器的設(shè)計要求（高通、帶阻、帶通）；</p>&

8、lt;p>　　轉(zhuǎn)換為模擬（高通、帶阻、帶通）濾波器的技術(shù)指標；</p><p>　　轉(zhuǎn)換為模擬低通濾波器的指標；</p><p>　　設(shè)計得到滿足第三步要求的低通濾波器傳遞函數(shù)；</p><p>　　通過頻率轉(zhuǎn)換得到模擬（高通、帶阻、帶通）濾波器；</p><p>　　變換為數(shù)字（高通、帶阻、帶通）濾波器。</p><

9、;p>　　在matlab中設(shè)計IIR濾波器的方法及其它們所用到的函數(shù)如表6－1所示。</p><p>　　表6－1 matlab中設(shè)計IIR濾波器的方法列表</p><p>　　（三） 雙線性變換法</p><p>　　與沖激響應不變法比較，雙線性變換的主要優(yōu)點是靠頻率的非線性關(guān)系得到S平面與Z平面的單值一一對應關(guān)系，整個值對應于單位圓一周。所以從模擬傳遞函數(shù)

10、可直接通過代數(shù)置換得到數(shù)字濾波器的傳遞函數(shù)。</p><p>　　【實例6-4】設(shè)計一個截止頻率為200Hz的數(shù)字低通濾波器，采用頻率為1000Hz。</p><p><b>　　解：源代碼如下：</b></p><p>　　[z,p,k]=buttap(3);</p><p>　　[b,a]=zp2tf(z,p,k);

11、</p><p>　　[bt,at]=lp2lp(b,a,200*2*pi);</p><p>　　[bz,az]=bilinear(bt,at,1000);</p><p>　　freqz(bz,az,512,1000)</p><p>　　結(jié)果如圖6－4所示：</p><p><b>　　bz =<

12、/b></p><p>　　0.0753 0.2259 0.2259 0.0753</p><p><b>　　az =</b></p><p>　　1.0000 -0.8266 0.5154 -0.0865</p><p>　　【實例6-5】基于Butterworth模擬濾波器原型

13、，使用雙線性轉(zhuǎn)換設(shè)計數(shù)字濾波器，其中參數(shù)指標為：</p><p>　　通帶截止頻率： 通帶波動值：</p><p>　　阻帶截止頻率： 阻帶波動值：</p><p>　　解：首先確定濾波器的階數(shù)N，同時根據(jù)確定＝0.5。</p><p>　　接著使用bilinear進行雙線性轉(zhuǎn)換，最后繪制在頻域上的各種圖像，其源代碼如下：</p>

14、<p><b>　　％數(shù)字濾波器指標</b></p><p>　　wp=0.2*pi;</p><p>　　ws=0.3*pi;</p><p><b>　　Rp=1;</b></p><p><b>　　As=15;</b></p><p&g

15、t;　?。?shù)字濾波器指標反轉(zhuǎn)變化為模擬濾波器的參數(shù)</p><p><b>　　T=1;</b></p><p><b>　　fs=1/T;</b></p><p>　　omegap=(2/T)*tan(wp/T);</p><p>　　omegas=(2/T)*tan(ws/T);</p&

16、gt;<p>　　ep=sqrt(10^(Rp/10)-1);</p><p>　　Ripple=sqrt(1/(1+ep*ep));</p><p>　　Attn=1/(10^(As/20));</p><p>　?。utterworth原型模擬濾波器的設(shè)計</p><p>　　[cs,ds]=afd_butt(omegap

17、,omegas,Rp,As);</p><p><b>　?。ルp線性變換</b></p><p>　　[b,a]= bilinear(cs,ds,T,fs);</p><p><b>　?。ヮl域圖像的繪制</b></p><p>　　freqz(b,a);</p><p>　

18、　程序運行后，產(chǎn)生4階的butterworth數(shù)字濾波器，頻率響應如圖6－5所示的波形。</p><p>　　圖6－5 4階Butterworth數(shù)字濾波器</p><p><b>　　實驗步驟</b></p><p>　　錄入一段語音信號，設(shè)置采樣率和采樣點數(shù)。</p><p><b>　　輸入信號的頻譜分析

19、</b></p><p>　　要求首先畫出輸入信號的時域波形；然后對該信號進行快速傅立葉變換，得到信號的頻譜特性。</p><p>　　設(shè)計數(shù)字低通濾波器并畫出數(shù)字低通濾波器的頻率響應。</p><p>　　用濾波器對信號進行濾波，比較濾波前后語音信號的波形及頻譜。</p><p>　　利用Matlab GUI工具做一個圖形用戶界

20、面，窗口簡潔大方。</p><p><b>　　實驗結(jié)果及截圖</b></p><p><b>　　結(jié)果</b></p><p><b>　　總結(jié)</b></p><p>　　通過這次課程設(shè)計，使自己對通信原理中的信號的變換，特別是傅立葉變換及其實現(xiàn)過程有了更深的了解。我對ma

21、tlab 軟件也有了熟練的了解，對工程項目設(shè)計積累的一定經(jīng)驗。起初在面對一個陌生時會有忌憚和不知所措，而通過這次訓練基本可以通過help 指令和doc 指令來自學函數(shù)用法。在程序調(diào)試的時候經(jīng)常會出現(xiàn)錯誤，耐心調(diào)試找到錯誤并搞清楚為什么會出錯則是作為一個程序員的一項基本素質(zhì)。</p><p>　　我們團隊分工明確，成員們勤勞能干，奠定我們高效地完成課程設(shè)計的基礎(chǔ)。而這期間，我們也確實遇到了很多困難，但我們都攻克了。

22、總之，我們在這次課程設(shè)計中學到了很多。</p><p><b>　　參考文獻：</b></p><p>　　【1】通信原理 第六版 樊昌信 曹麗娜編著 國防工業(yè)出版社</p><p>　　【2】信號與系統(tǒng)引論 鄭君里 應啟 楊為理編著 高等教育出版社</p><p>　　【3】 MATLAB及在電子信息類課程

23、中的應用 第二版 唐向宏 岳恒立 鄭雪峰 編著 電子工業(yè)出版社</p><p><b>　　附錄：</b></p><p><b>　　源程序</b></p><p>　　function varargout = idea2(varargin)</p><p>　　% IDEA2 M-file f

24、or idea2.fig</p><p>　　% IDEA2, by itself, creates a new IDEA2 or raises the existing</p><p>　　% singleton*.</p><p><b>　　%</b></p><p>　　% H =

25、IDEA2 returns the handle to a new IDEA2 or the handle to</p><p>　　% the existing singleton*.</p><p><b>　　%</b></p><p>　　% IDEA2('CALLBACK',hObject,even

26、tData,handles,...) calls the local</p><p>　　% function named CALLBACK in IDEA2.M with the given input arguments.</p><p><b>　　%</b></p><p>　　% IDEA2('Proper

27、ty','Value',...) creates a new IDEA2 or raises the</p><p>　　% existing singleton*. Starting from the left, property value pairs are</p><p>　　% applied to the GUI before id

28、ea2_OpeningFunction gets called. An</p><p>　　% unrecognized property name or invalid value makes property application</p><p>　　% stop. All inputs are passed to idea2_OpeningFcn via v

29、arargin.</p><p><b>　　%</b></p><p>　　% *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one</p><p>　　% instance to run (singleton)"

30、;.</p><p><b>　　%</b></p><p>　　% See also: GUIDE, GUIDATA, GUIHANDLES</p><p>　　% Copyright 2002-2003 The MathWorks, Inc.</p><p>　　% Edit the above text to mo

31、dify the response to help idea2</p><p>　　% Last Modified by GUIDE v2.5 13-Jun-2012 16:27:25</p><p>　　% Begin initialization code - DO NOT EDIT</p><p>　　gui_Singleton = 1;</p>

32、<p>　　gui_State = struct('gui_Name', mfilename, ...</p><p>　　'gui_Singleton', gui_Singleton, ...</p><p>　　'gui_OpeningFcn', @idea2_OpeningFcn, ...</p>

33、<p>　　'gui_OutputFcn', @idea2_OutputFcn, ...</p><p>　　'gui_LayoutFcn', [] , ...</p><p>　　'gui_Callback', []);</p><p>　　if nargin && ischar(

34、varargin{1})</p><p>　　gui_State.gui_Callback = str2func(varargin{1});</p><p><b>　　end</b></p><p>　　if nargout</p><p>　　[varargout{1:nargout}] = gui_mainfcn(

35、gui_State, varargin{:});</p><p><b>　　else</b></p><p>　　gui_mainfcn(gui_State, varargin{:});</p><p><b>　　end</b></p><p>　　% End initialization co

36、de - DO NOT EDIT</p><p>　　% --- Executes just before idea2 is made visible.</p><p>　　function idea2_OpeningFcn(hObject, eventdata, handles, varargin)</p><p>　　% This function has no

37、 output args, see OutputFcn.</p><p>　　% hObject handle to figure</p><p>　　% eventdata reserved - to be defined in a future version of MATLAB</p><p>　　% handles structure with

38、 handles and user data (see GUIDATA)</p><p>　　% varargin command line arguments to idea2 (see VARARGIN)</p><p>　　% Choose default command line output for idea2</p><p>　　handles.ou

39、tput = hObject;</p><p>　　% Update handles structure</p><p>　　guidata(hObject, handles);</p><p>　　global SampleRate</p><p>　　SampleRate=44100;</p><p>　　glo

40、bal SecLength;</p><p>　　SecLength=60; %因為wav文件是十分鐘的，不能全部讀入，所以將其分為10段，一段60秒</p><p>　　global SecNum;</p><p>　　SecNum=600/SecLength;</p><p>　　global CurrentSecti

41、on</p><p>　　% UIWAIT makes idea2 wait for user response (see UIRESUME)</p><p>　　% uiwait(handles.figure1);</p><p>　　% --- Outputs from this function are returned to the command line

42、.</p><p>　　function varargout = idea2_OutputFcn(hObject, eventdata, handles) </p><p>　　% varargout cell array for returning output args (see VARARGOUT);</p><p>　　% hObject handl

43、e to figure</p><p>　　% eventdata reserved - to be defined in a future version of MATLAB</p><p>　　% handles structure with handles and user data (see GUIDATA)</p><p>　　% Get defa

44、ult command line output from handles structure</p><p>　　varargout{1} = handles.output;</p><p>　　% --- Executes on button press in pushbutton1.</p><p>　　function Insert_Callback(hObj

45、ect, eventdata, handles)</p><p>　　% hObject handle to pushbutton1 (see GCBO)</p><p>　　% eventdata reserved - to be defined in a future version of MATLAB</p><p>　　% handles st

46、ructure with handles and user data (see GUIDATA)</p><p>　　global SecLength;</p><p>　　global SecNum;</p><p>　　global CurrentSection;</p><p>　　global SampleRate;</p&g

47、t;<p>　　global Myrate;</p><p>　　global SoundData;</p><p>　　global Mycode;</p><p>　　SampleRate = Myrate;</p><p>　　%Sam =Myrate;</p><p>　　%disp(Sam);&

48、lt;/p><p>　　[FileName,PathName] = uigetfile('*.wav');</p><p>　　data=wavread(strcat(PathName,FileName),[SecLength*(CurrentSection-1)*(SampleRate)+1,SecLength*CurrentSection*SampleRate]);<

49、/p><p>　　SoundData=data(1:10:length(data),1);</p><p>　　%N= Mycode;</p><p>　　n= 0:length(SoundData)-1;</p><p>　　t= n/SampleRate;</p><p>　　figure,plot(t,SoundDa

50、ta);</p><p>　　%CurrentSection = CurrentSection + 1</p><p>　　% --- Executes on button press in pushbutton2.</p><p>　　function Sure_Callback(hObject, eventdata, handles)</p>&l

51、t;p>　　% hObject handle to pushbutton2 (see GCBO)</p><p>　　% eventdata reserved - to be defined in a future version of MATLAB</p><p>　　% handles structure with handles and user data (se

52、e GUIDATA)</p><p>　　global Myrate ;</p><p>　　global Mycode;</p><p>　　Mycode = str2num(get(handles.edit1,'string'));</p><p>　　Myrate = str2num(get(findobj('t

53、ag','edit2'),'string'));</p><p>　　%disp(Myrate);</p><p>　　% --- Executes on button press in pushbutton3.</p><p>　　function DFT_Callback(hObject, eventdata, handl

54、es)</p><p>　　% hObject handle to pushbutton3 (see GCBO)</p><p>　　% eventdata reserved - to be defined in a future version of MATLAB</p><p>　　% handles structure with handles

55、and user data (see GUIDATA)</p><p>　　global SampleRate;</p><p>　　global SoundData;</p><p>　　global Mycode;</p><p><b>　　global y;</b></p><p>　　N

56、 = Mycode;</p><p><b>　　%disp(N);</b></p><p>　　y = fft(SoundData,N);%進行fft變換</p><p>　　mag=abs(y);%求幅值</p><p>　　n= 1:N;f = n*SampleRate/length(y);%進行對應的頻率轉(zhuǎn)換&l

57、t;/p><p><b>　　figure;</b></p><p>　　plot(f,mag);%做頻譜圖</p><p>　　xlabel('頻率(Hz)');</p><p>　　ylabel('幅值');</p><p>　　title('頻譜'

58、);</p><p><b>　　grid;</b></p><p>　　%Py =2*(y.*conj(y))/N; %計算功率譜密度Py</p><p><b>　　%figure;</b></p><p>　　%plot(f,Py);</p><p>　　%xla

59、bel('頻率(Hz)');</p><p>　　%ylabel('功率譜密度');</p><p>　　%title('x=120.0*cos(2*3.14*sf*t)功率譜密度');</p><p><b>　　%grid;</b></p><p>　　%point =

60、 get(handles.edit1,'string');</p><p>　　%fs = SampleRate; </p><p>　　%m = length(SoundData); % Window length</p><p>　　%n = pow2(nextpow2(m)); % Transform length</p><

61、;p>　　%y = fft(SoundData,N); % DFT%做n點的FFT，n為2的多少次冪，次冪值為大于或等于length(SoundData)的絕對值的2的最小次冪</p><p>　　%n = (0:N-1) % Frequency range%圖上顯示的最高頻率為采樣頻率</p><p>　　%power = y.*conj(y)/n; % Power of the

62、DFT%計算復共軛值</p><p>　　%f = n*SampleRate/2;</p><p>　　%figure,plot(f,y);</p><p>　　%xlabel('Frequency (Hz)');</p><p>　　%ylabel('DFT');</p><p>　　

63、%title('{\bf 傅里葉變換}');</p><p>　　%pause(1);</p><p>　　% --- Executes on button press in pushbutton4.</p><p>　　function DTpass_Callback(hObject, eventdata, handles)</p>

64、<p>　　% hObject handle to pushbutton4 (see GCBO)</p><p>　　% eventdata reserved - to be defined in a future version of MATLAB</p><p>　　% handles structure with handles and user data (

65、see GUIDATA)</p><p><b>　　global y;</b></p><p>　　global SampleRate;</p><p>　　%%根據(jù)原始信號與噪聲信號頻譜圖，選擇巴特沃斯低通濾波器</p><p>　　%%濾波器性能指標：取wp=1500*2*pi rad/s,ws=2400*2*pi

66、 rad/s,Rp=2,Rs=15</p><p>　　%選擇濾波器的最小階數(shù)</p><p>　　wp=1000*2*pi;ws=6000*2*pi;Rp=2;Rs=15;</p><p>　　[N1,Wn]=buttord(wp,ws,Rp,Rs,'s');</p><p>　　%創(chuàng)建butterworth模擬濾波器<

67、;/p><p>　　[Z,P,K]=buttap(N1);</p><p>　　%把濾波器的零極點模型轉(zhuǎn)化為傳遞函數(shù)模型</p><p>　　[Bap,Aap]=zp2tf(Z,P,K);</p><p>　　%把模擬濾波器原型轉(zhuǎn)換成截止頻率為Wn的低通濾波器</p><p>　　[b,a]=lp2lp(Bap,Aap,W

68、n);</p><p>　　%用雙線性變換法實現(xiàn)模擬濾波器到數(shù)字濾波器的轉(zhuǎn)換</p><p>　　[bz,az]=bilinear(b,a,SampleRate);</p><p><b>　　%繪制頻率響應曲線</b></p><p>　　[H,W]=freqz(bz,az);</p><p>

69、　　figure(3);</p><p>　　subplot(211),</p><p>　　plot(W*SampleRate/(2*pi),abs(H)); grid </p><p>　　title('butterworth低通濾波器頻響圖');xlabel('頻率/Hz'); ylabel('頻率響應幅度'

70、;);</p><p>　　%%%%%%%用所設(shè)計的濾波器對語音+噪聲信號y濾波%%%%%%%</p><p>　　global Mycode;</p><p>　　y1=filter(bz,az,y);</p><p>　　n=(0:Mycode-1)*SampleRate/length(y1);</p><p>　

71、　subplot(212),plot(n,y1);title('濾波后信號波形圖');</p><p>　　xlabel('頻率/Hz'); ylabel('頻率響應幅度');</p><p>　　%sound(x,fs);%播放語音信號</p><p>　　%sound(y,fs);%播放語音+噪聲信號</p

72、><p>　　%sound(y1,fs);%播放濾波后的信號</p><p>　　% --- Executes on button press in pushbutton5.</p><p>　　function quit_Callback(hObject, eventdata, handles)</p><p>　　% hObject h

73、andle to pushbutton5 (see GCBO)</p><p>　　% eventdata reserved - to be defined in a future version of MATLAB</p><p>　　% handles structure with handles and user data (see GUIDATA)</p>&

74、lt;p><b>　　close;</b></p><p>　　msgbox('謝謝使用，再見！');</p><p>　　function edit1_Callback(hObject, eventdata, handles)</p><p>　　% hObject handle to edit1 (see GCB

75、O)</p><p>　　% eventdata reserved - to be defined in a future version of MATLAB</p><p>　　% handles structure with handles and user data (see GUIDATA)</p><p>　　% Hints: get(hObjec

76、t,'String') returns contents of edit1 as text</p><p>　　% str2double(get(hObject,'String')) returns contents of edit1 as a double</p><p>　　%A=get(hObject,'string');

77、</p><p>　　%guidata(hObject,handles); %必須要加這一句</p><p>　　% --- Executes during object creation, after setting all properties.</p><p>　　function edit1_CreateFcn(hObject, eventdata, han

78、dles)</p><p>　　% hObject handle to edit1 (see GCBO)</p><p>　　% eventdata reserved - to be defined in a future version of MATLAB</p><p>　　% handles empty - handles not created

79、 until after all CreateFcns called</p><p>　　% Hint: edit controls usually have a white background on Windows.</p><p>　　% See ISPC and COMPUTER.</p><p><b>　　if ispc</b

80、></p><p>　　set(hObject,'BackgroundColor','white');</p><p><b>　　else</b></p><p>　　set(hObject,'BackgroundColor',get(0,'defaultUicontrolBack

81、groundColor'));</p><p><b>　　end</b></p><p>　　function edit2_Callback(hObject, eventdata, handles)</p><p>　　% hObject handle to edit2 (see GCBO)</p><p&g

82、t;　　% eventdata reserved - to be defined in a future version of MATLAB</p><p>　　% handles structure with handles and user data (see GUIDATA)</p><p>　　% Hints: get(hObject,'String') r

83、eturns contents of edit2 as text</p><p>　　% str2double(get(hObject,'String')) returns contents of edit2 as a double</p><p>　　% --- Executes during object creation, after setting a

84、ll properties.</p><p>　　function edit2_CreateFcn(hObject, eventdata, handles)</p><p>　　% hObject handle to edit2 (see GCBO)</p><p>　　% eventdata reserved - to be defined in a fu

85、ture version of MATLAB</p><p>　　% handles empty - handles not created until after all CreateFcns called</p><p>　　% Hint: edit controls usually have a white background on Windows.</p>&

86、lt;p>　　% See ISPC and COMPUTER.</p><p><b>　　if ispc</b></p><p>　　set(hObject,'BackgroundColor','white');</p><p><b>　　else</b></p&