實驗四IIR和FIR數(shù)字濾波器設(shè)計及軟件實現(xiàn)實驗報告

1、.數(shù)字信號處理實驗報告實驗四iir數(shù)字濾波器設(shè)計及軟件實現(xiàn)(一)fir數(shù)字濾波器設(shè)計及軟件實現(xiàn)（二）2018 年 11 月 28 日.一、實驗?zāi)康?實驗4_1)（1）熟悉用雙線性變換法設(shè)計iir數(shù)字濾波器的原理與方法；（2）學(xué)會調(diào)用matlab信號處理工具箱中濾波器設(shè)計函數(shù)（或濾波器設(shè)計分析工具fdatool）設(shè)計各種iir數(shù)字濾波器，學(xué)會根據(jù)濾波需求確定濾波器指標(biāo)參數(shù)。（3）掌握iir數(shù)字濾波器的matlab實現(xiàn)方法。（4）通過觀察濾波器輸入輸出信號的時域波形及其頻譜，建立數(shù)字濾波的概念。（實驗4_2）（1）掌握用窗函數(shù)法設(shè)計fir數(shù)字濾波器的原理和方法。（2）掌握用等波紋最佳逼近法設(shè)計fi

2、r數(shù)字濾波器的原理和方法。（3）掌握fir濾波器的快速卷積實現(xiàn)原理。（4）學(xué)會調(diào)用matlab函數(shù)設(shè)計與實現(xiàn)fir濾波器。二、實驗原理與方法（實驗4_1）設(shè)計iir數(shù)字濾波器一般采用間接法（脈沖響應(yīng)不變法和雙線性變換法），應(yīng)用最廣泛的是雙線性變換法?；驹O(shè)計過程是：先將給定的數(shù)字濾波器的指標(biāo)轉(zhuǎn)換成過渡模擬濾波器的指標(biāo)； 設(shè)計過渡模擬濾波器；將過渡模擬濾波器系統(tǒng)函數(shù)轉(zhuǎn)換成數(shù)字濾波器的系統(tǒng)函數(shù)。matlab信號處理工具箱中的各種iir數(shù)字濾波器設(shè)計函數(shù)都是采用雙線性變換法。第六章介紹的濾波器設(shè)計函數(shù)butter、cheby1 、cheby2 和ellip可以分別被調(diào)用來直接設(shè)計巴特沃斯、切比雪夫1

3、、切比雪夫2和橢圓模擬和數(shù)字濾波器。本實驗要求讀者調(diào)用如上函數(shù)直接設(shè)計iir數(shù)字濾波器。本實驗的數(shù)字濾波器的matlab實現(xiàn)是指調(diào)用matlab信號處理工具箱函數(shù)filter對給定的輸入信號x(n)進(jìn)行濾波，得到濾波后的輸出信號y(n）。3、 實驗內(nèi)容及步驟（實驗4_1）（1）調(diào)用信號產(chǎn)生函數(shù)mstg產(chǎn)生由三路抑制載波調(diào)幅信號相加構(gòu)成的復(fù)合信號st，該函數(shù)還會自動繪圖顯示st的時域波形和幅頻特性曲線，如圖1所示。由圖可見，三路信號時域混疊無法在時域分離。但頻域是分離的，所以可以通過濾波的方法在頻域分離，這就是本實驗的目的。圖1 三路調(diào)幅信號st的時域波形和幅頻特性曲線（2）要求將st中三路調(diào)幅

4、信號分離，通過觀察st的幅頻特性曲線，分別確定可以分離st中三路抑制載波單頻調(diào)幅信號的三個濾波器（低通濾波器、帶通濾波器、高通濾波器）的通帶截止頻率和阻帶截止頻率。要求濾波器的通帶最大衰減為0.1db,阻帶最小衰減為60db。提示：抑制載波單頻調(diào)幅信號的數(shù)學(xué)表示式為其中，稱為載波，fc為載波頻率，稱為單頻調(diào)制信號，f0為調(diào)制正弦波信號頻率，且滿足。由上式可見，所謂抑制載波單頻調(diào)幅信號，就是2個正弦信號相乘，它有2個頻率成分：和頻和差頻，這2個頻率成分關(guān)于載波頻率fc對稱。所以，1路抑制載波單頻調(diào)幅信號的頻譜圖是關(guān)于載波頻率fc對稱的2根譜線，其中沒有載頻成分，故取名為抑制載波單頻調(diào)幅信號。容易

5、看出，圖1中三路調(diào)幅信號的載波頻率分別為250hz、500hz、1000hz。如果調(diào)制信號m(t)具有帶限連續(xù)頻譜，無直流成分，則就是一般的抑制載波調(diào)幅信號。其頻譜圖是關(guān)于載波頻率fc對稱的2個邊帶（上下邊帶），在專業(yè)課通信原理中稱為雙邊帶抑制載波 (dsb-sc) 調(diào)幅信號,簡稱雙邊帶 (dsb) 信號。如果調(diào)制信號m(t)有直流成分，則就是一般的雙邊帶調(diào)幅信號。其頻譜圖是關(guān)于載波頻率fc對稱的2個邊帶（上下邊帶），并包含載頻成分。（3）編程序調(diào)用matlab濾波器設(shè)計函數(shù)ellipord和ellip分別設(shè)計這三個橢圓濾波器，并繪圖顯示其幅頻響應(yīng)特性曲線。 （4）調(diào)用濾波器實現(xiàn)函數(shù)filte

6、r，用三個濾波器分別對信號產(chǎn)生函數(shù)mstg產(chǎn)生的信號st進(jìn)行濾波，分離出st中的三路不同載波頻率的調(diào)幅信號y1(n)、y2(n)和y3(n)， 并繪圖顯示y1(n)、y2(n)和y3(n)的時域波形，觀察分離效果。(實驗4_2)（1）認(rèn)真復(fù)習(xí)第七章中用窗函數(shù)法和等波紋最佳逼近法設(shè)計fir數(shù)字濾波器的原理；（2）調(diào)用信號產(chǎn)生函數(shù)xtg產(chǎn)生具有加性噪聲的信號xt，并自動顯示xt及其頻譜，如圖所示；圖 具有加性噪聲的信號x(t)及其頻譜如圖（3）請設(shè)計低通濾波器，從高頻噪聲中提取xt中的單頻調(diào)幅信號，要求信號幅頻失真小于0.1db，將噪聲頻譜衰減60db。先觀察xt的頻譜，確定濾波器指標(biāo)參數(shù)。（4）

7、根據(jù)濾波器指標(biāo)選擇合適的窗函數(shù)，計算窗函數(shù)的長度n，調(diào)用matlab函數(shù)fir1設(shè)計一個fir低通濾波器。并編寫程序，調(diào)用matlab快速卷積函數(shù)fftfilt實現(xiàn)對xt的濾波。繪圖顯示濾波器的頻響特性曲線、濾波器輸出信號的幅頻特性圖和時域波形圖。（4）重復(fù)（3），濾波器指標(biāo)不變，但改用等波紋最佳逼近法，調(diào)用matlab函數(shù)remezord和remez設(shè)計fir數(shù)字濾波器。并比較兩種設(shè)計方法設(shè)計的濾波器階數(shù)。提示：matlab函數(shù)fir1和fftfilt的功能及其調(diào)用格式請查閱本書第7章和第8章；采樣頻率fs=1000hz，采樣周期t=1/fs；根據(jù)圖10.6.1(b)和實驗要求，可選擇濾波器

8、指標(biāo)參數(shù)：通帶截止頻率fp=120hz，阻帶截至頻率fs=150hz，換算成數(shù)字頻率，通帶截止頻率，通帶最大衰為0.1db，阻帶截至頻率，阻帶最小衰為60db。實驗程序框圖如圖10.5.2所示，供讀者參考。fs=1000，t=1/fsxt=xtg產(chǎn)生信號xt, 并顯示xt及其頻譜用窗函數(shù)法或等波紋最佳逼近法設(shè)計fir濾波器hn對信號xt濾波：yt=fftfilt(hn,xt)1、計算并繪圖顯示濾波器損耗函數(shù)2、繪圖顯示濾波器輸出信號ytend4、 matlab源代碼、實驗結(jié)果圖像和結(jié)果分析（實驗4_1）實驗代碼（函數(shù)部分單獨分出來）：（mstg函數(shù)）：function st=mstg%產(chǎn)生信號

9、序列變量st,并顯示st的時域波形和頻譜%st=mstg返回三路調(diào)幅信號相加形成的混合信號，長度n=1600n=1600;fs=10000;t=1/fs;tp=n*t;t=0:t:(n-1)*t;k=0:n-1;f=k/tp;fc1=fs/10;fm1=fc1/10;fc2=fs/20;fm2=fc2/10;fc3=fs/40;fm3=fc3/10;xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t); %xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t);xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t);st=xt1+xt2+

10、xt3;fxt=fft(st,n);%subplot(3,1,1);plot(t,st);grid;xlabel('t/s');ylabel('s(t)');axis(0,tp/8,min(st),max(st);title('(a) s(t)的波形');subplot(3,1,2);stem(f,abs(fxt)/max(abs(fxt),'.');grid;title('(b) s(t)的頻譜');axis(0,fs/5,0,1.2);xlabel('f/hz');ylabel('幅度&

11、#39;);（myplot函數(shù)）：function myplot(b,a)%時域離散系統(tǒng)損耗函數(shù)繪圖%b為系統(tǒng)函數(shù)分子多項式系數(shù)向量%a為系統(tǒng)函數(shù)分母多項式系數(shù)向量h,w=freqz(b,a,1000);m=abs(h);plot(w/pi,20*log10(m/max(m);grid on;xlabel('omega/pi');ylabel('幅度(db)')axis(0,1,-150,50);title('損耗函數(shù)曲線');（tplot函數(shù)）：function tplot(xn,t,yn)%時域序列連續(xù)曲線繪圖函數(shù)% xn:信號數(shù)據(jù)序列，yn

12、:繪圖信號的縱坐標(biāo)名稱（字符串）% t為采樣間隔n=0:length(xn)-1;t=n*t;plot(t,xn);xlabel('t/s');ylabel(yn);axis(0,t(end),min(xn),1.2*max(xn)（實驗4_1源代碼）：% iir數(shù)字濾波器設(shè)計及軟件實現(xiàn)clear all;close allfs=10000;t=1/fs; %采樣頻率%調(diào)用信號產(chǎn)生函數(shù)mstg產(chǎn)生由三路抑制載波調(diào)幅信號相加構(gòu)成的復(fù)合信號st st=mstg;%低通濾波器設(shè)計與實現(xiàn)=fp=280;fs=450;wp=2*fp/fs;ws=2*fs/fs;rp=0.1;rs=60;

13、 %df指標(biāo)（低通濾波器的通、阻帶邊界頻）n,wp=ellipord(wp,ws,rp,rs); %調(diào)用ellipord計算橢圓df階數(shù)n和通帶截止頻率wpb,a=ellip(n,rp,rs,wp); %調(diào)用ellip計算橢圓帶通df系統(tǒng)函數(shù)系數(shù)向量b和ay1t=filter(b,a,st); %濾波器軟件實現(xiàn)% 低通濾波器設(shè)計與實現(xiàn)繪圖部分figure(5);subplot(2,1,1);myplot(b,a); %調(diào)用繪圖函數(shù)myplot繪制損耗函數(shù)曲線yt='y_1(t)'subplot(2,1,2);tplot(y1t,t,yt); %調(diào)用繪圖函數(shù)tplot繪制濾波器輸

14、出波形%帶通濾波器設(shè)計與實現(xiàn)=fpl=440;fpu=560;fsl=275;fsu=900;wp=2*fpl/fs,2*fpu/fs;ws=2*fsl/fs,2*fsu/fs;rp=0.1;rs=60; n,wp=ellipord(wp,ws,rp,rs); %調(diào)用ellipord計算橢圓df階數(shù)n和通帶截止頻率wpb,a=ellip(n,rp,rs,wp); %調(diào)用ellip計算橢圓帶通df系統(tǒng)函數(shù)系數(shù)向量b和ay2t=filter(b,a,st); %濾波器軟件實現(xiàn)% 帶通濾波器設(shè)計與實現(xiàn)繪圖部分figure(3);subplot(2,1,1);myplot(b,a); %調(diào)用繪圖函數(shù)m

15、yplot繪制損耗函數(shù)曲線yt='y_2(t)'subplot(2,1,2);tplot(y2t,t,yt); %調(diào)用繪圖函數(shù)tplot繪制濾波器輸出波形%高通濾波器設(shè)計與實現(xiàn)=fp=890;fs=600;wp=2*fp/fs;ws=2*fs/fs;rp=0.1;rs=60; %df指標(biāo)（低通濾波器的通、阻帶邊界頻）n,wp=ellipord(wp,ws,rp,rs); %調(diào)用ellipord計算橢圓df階數(shù)n和通帶截止頻率wpb,a=ellip(n,rp,rs,wp,'high'); %調(diào)用ellip計算橢圓帶通df系統(tǒng)函數(shù)系數(shù)向量b和ay3t=filter(

16、b,a,st); %濾波器軟件實現(xiàn)% 高低通濾波器設(shè)計與實現(xiàn)繪圖部分figure(4);subplot(2,1,1);myplot(b,a); %調(diào)用繪圖函數(shù)myplot繪制損耗函數(shù)曲線yt='y_3(t)'subplot(2,1,2);tplot(y3t,t,yt); %調(diào)用繪圖函數(shù)tplot繪制濾波器輸出波形 function myplot(b,a)%時域離散系統(tǒng)損耗函數(shù)繪圖%b為系統(tǒng)函數(shù)分子多項式系數(shù)向量%a為系統(tǒng)函數(shù)分母多項式系數(shù)向量h,w=freqz(b,a,1000);m=abs(h);plot(w/pi,20*log10(m/max(m);grid on;xlab

17、el('omega/pi');ylabel('幅度(db)')axis(0,1,-80,5);title('損耗函數(shù)曲線'); function tplot(xn,t,yn)%時域序列連續(xù)曲線繪圖函數(shù)% xn:信號數(shù)據(jù)序列，yn:繪圖信號的縱坐標(biāo)名稱（字符串）% t為采樣間隔n=0:length(xn)-1;t=n*t;plot(t,xn);xlabel('t/s');ylabel(yn);axis(0,t(end),min(xn),1.2*max(xn) function st=mstgn=2000;fs=10000;t=1/f

18、s;tp=n*t;t=0:t:(n-1)*t;k=0:n-1;f=k/tp;fc1=fs/10;fm1=fc1/10;fc2=fs/20;fm2=fc2/10;fc3=fs/40;fm3=fc3/10;xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t);xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t);xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t);st=xt1+xt2+xt3;fxt=fft(st,n);subplot(3,1,1)plot(t,st);grid;xlabel('t/s');ylabel(

19、's(t)');axis(0,tp/8,min(st),max(st);title('(a) s(t)的波形')subplot(3,1,2)stem(f,abs(fxt)/max(abs(fxt),'.');grid;title('(b) s(t)的頻譜')axis(0,fs/5,0,1.2);xlabel('f/hz');ylabel('幅度');（實驗4_2）（xtg函數(shù)）：function xt=xtg(n)%實驗五信號x(t)產(chǎn)生,并顯示信號的幅頻特性曲線%xt=xtg(n) 產(chǎn)生一個長度為

20、n,有加性高頻噪聲的單頻調(diào)幅信號xt,采樣頻率fs=1000hz%載波頻率fc=fs/10=100hz,調(diào)制正弦波頻率f0=fc/10=10hz.n=2000;fs=1000;t=1/fs;tp=n*t;t=0:t:(n-1)*t;fc=fs/10;f0=fc/10; %載波頻率fc=fs/10，單頻調(diào)制信號頻率為f0=fc/10;mt=cos(2*pi*f0*t); %產(chǎn)生單頻正弦波調(diào)制信號mt，頻率為f0ct=cos(2*pi*fc*t); %產(chǎn)生載波正弦波信號ct，頻率為fcxt=mt.*ct; %相乘產(chǎn)生單頻調(diào)制信號xtnt=2*rand(1,n)-1; %產(chǎn)生隨機噪聲nt%=設(shè)計高通

21、濾波器hn,用于濾除噪聲nt中的低頻成分,生成高通噪聲=fp=150; fs=200;rp=0.1;as=70; % 濾波器指標(biāo)fb=fp,fs;m=0,1; % 計算remezord函數(shù)所需參數(shù)f,m,devdev=10(-as/20),(10(rp/20)-1)/(10(rp/20)+1);n,fo,mo,w=remezord(fb,m,dev,fs); % 確定remez函數(shù)所需參數(shù)hn=remez(n,fo,mo,w); % 調(diào)用remez函數(shù)進(jìn)行設(shè)計,用于濾除噪聲nt中的低頻成分yt=filter(hn,1,10*nt); %濾除隨機噪聲中低頻成分，生成高通噪聲yt%=xt=xt+y

22、t; %噪聲加信號fst=fft(xt,n);k=0:n-1;f=k/tp;subplot(3,1,1);plot(t,xt);grid;xlabel('t/s');ylabel('x(t)');axis(0,tp/5,min(xt),max(xt);title('(a) 信號加噪聲波形')subplot(3,1,2);plot(f,abs(fst)/max(abs(fst);grid;title('(b) 信號加噪聲的頻譜')axis(0,fs/2,0,1.2);xlabel('f/hz');ylabel(

23、9;幅度')（主程序源代碼）：clear all;close all;%=調(diào)用xtg產(chǎn)生信號xt, xt長度n=1000,并顯示xt及其頻譜,=n=1000;xt=xtg(n);fp=120; fs=150;rp=0.2;as=60;fs=1000; t=1/fs; % 輸入給定指標(biāo)% (1) 用窗函數(shù)法設(shè)計濾波器wc=(fp+fs)/fs; %理想低通濾波器截止頻率(關(guān)于pi歸一化）b=2*pi*(fs-fp)/fs; %過渡帶寬度指標(biāo)nb=ceil(11*pi/b); %blackman窗的長度nhn=fir1(nb-1,wc,blackman(nb);hw=abs(fft(hn,

24、1024); % 求設(shè)計的濾波器頻率特性ywt=fftfilt(hn,xt,n); %調(diào)用函數(shù)fftfilt對xt濾波figure(2);subplot(3,1,1);myplot(hn,xt); %調(diào)用繪圖函數(shù)myplot繪制損耗函數(shù)曲線y1t='y_w(t)'subplot(3,1,2);tplot(ywt,t,y1t);% (2) 用等波紋最佳逼近法設(shè)計濾波器fb=fp,fs;m=1,0; % 確定remezord函數(shù)所需參數(shù)f,m,devdev=(10(rp/20)-1)/(10(rp/20)+1),10(-as/20);ne,fo,mo,w=remezord(fb,m

25、,dev,fs); % 確定remez函數(shù)所需參數(shù)hn=remez(ne,fo,mo,w); % 調(diào)用remez函數(shù)進(jìn)行設(shè)計hw=abs(fft(hn,1024); % 求設(shè)計的濾波器頻率特性yet=fftfilt(hn,xt,n); % 調(diào)用函數(shù)fftfilt對xt濾波figure(3);subplot(3,1,1);myplot(hn,xt); %調(diào)用繪圖函數(shù)myplot繪制損耗函數(shù)曲線y2t='y_e(t)'subplot(3,1,2);tplot(yet,t,y2t)實驗圖像：（如下圖）實驗4_1圖像：低通濾波器損耗函數(shù)及其分離出的調(diào)幅信號y1(t)實驗4_1圖像：帶通

26、濾波器損耗函數(shù)及其分離出的調(diào)幅信號y2(t)實驗4_1圖像：高通濾波器損耗函數(shù)及其分離出的調(diào)幅信號y3(t)實驗4_2圖像實驗4_2圖像5、 思考題及解答（實驗4_1思考題）（1）請閱讀信號產(chǎn)生函數(shù)mstg，確定三路調(diào)幅信號的載波頻率和調(diào)制信號頻率。（2）信號產(chǎn)生函數(shù)mstg中采樣點數(shù)n=800，對st進(jìn)行n點fft可以得到6根理想譜線。如果取n=1000，可否得到6根理想譜線？為什么？n=2000呢？請改變函數(shù)mstg中采樣點數(shù)n的值，觀察頻譜圖驗證您的判斷是否正確。（3）修改信號產(chǎn)生函數(shù)mstg，給每路調(diào)服信號加入載波成分，產(chǎn)生調(diào)幅（am）信號，重復(fù)本實驗，觀察am信號與抑制載波調(diào)幅信號的

27、時域波形及其頻譜的差別。提示：am信號表示式： s(t)=ad+amcos(2f0t)cos(2fct) adam（實驗4_2思考題）(1)如果給定通帶截止頻率和阻帶截止頻率以及阻帶最小衰減,如何用窗函數(shù)法設(shè)計線性相位低通濾波器?請寫出設(shè)計步驟.(2)如果要求用窗函數(shù)法設(shè)計帶通濾波器,且給定通帶上、下截止頻率為和，阻帶上、下截止頻率為和，試求理想帶通濾波器的截止頻率。（3）解釋為什么對同樣的技術(shù)指標(biāo)，用等波紋最佳逼近法設(shè)計的濾波器階數(shù)低？解答：（實驗4_1）思考題（1）第一路調(diào)幅信號的載波頻率fc1=1000hz第一路調(diào)幅信號的調(diào)制信號頻率fm1=100hz第二路調(diào)幅信號的載波信號頻率fc2=500hz第二路調(diào)幅信號的調(diào)制信號頻率fm2=500hz第三路調(diào)幅信號的載波頻率fc3=250hz第三路調(diào)幅信號的調(diào)制信號頻率fm3=25hz思考題（2） 因為信號st是周期序列，譜分析時要求觀察時間為整數(shù)倍周期。所以，本題的一般解答方法是，先確定信號st的周期，在判斷所給采樣點數(shù)n對應(yīng)的觀察時間tp=n