信息工程專業(yè)數(shù)字信號處理課程實驗2報告

1、數(shù)字信號處理實驗報告：劉誠征班級：信息 11 班 學號：2010052072頻率采樣型濾波器實驗二頻率采樣型濾波器一、實驗目的1.通過該實驗學會使用頻率采樣型結構實現(xiàn)FIR濾波器，初步熟悉FIR濾波器的線性相位特點。2.通過該實驗直觀體會頻率采樣型濾波器所具有的“濾波器組”特性，即在并聯(lián)結構的每條支可以分別得到輸入信號的各次諧波。3.通過該實驗學會如何使用周期沖激串檢測所實現(xiàn)濾波器的頻域響應。二、實驗內容頻率采樣型濾波器是由一個梳狀濾波器和若干路諧振器的，可用公式表述如下：其中r值理論上為1，實際中取非常接近1的值。為了使系數(shù)為實數(shù)，可以將諧振器的共軛復根合并，不失一般性，假設N為偶數(shù)，于是可

2、以得到如圖1所示的結構。圖1 N為偶數(shù)的實系數(shù)頻率采樣型結構濾波器其中，。實驗中假設頻率采樣型濾波器階數(shù)N =16 。三 實驗結果及分析1. 構造濾波器輸入信號頻率，其中，基波。設時域信號s(t)的采樣頻率fs =N f0，繪制出采樣時刻從0到，L-1的采樣信號波形，其中采樣點數(shù)為L=2N，確認時域信號采樣正確。通過將正弦波形與采樣序列同時繪制出來，可以直觀確認采樣信號正確。2. 對采樣信號的第二個周期（n = N, N +1, , L 1）進行離散出幅頻特性和相頻特性圖，觀察并分析其特點。變換，畫幅頻特性圖顯示4至12點對應幅值為零，第3點和13點最大值，符合共軛對稱關系。相頻特性圖反應了相

3、位對應情況。3.計算濾波器抽頭系數(shù)h(n), n =0,1, N析其幅頻特性和相頻特性。1;畫出該濾波器的頻譜圖，觀察并分抽頭系數(shù)。幅頻特性顯示該濾波器為帶阻濾波器。相頻特性圖如上圖所示。4. 編程實現(xiàn)圖 1 所示的頻率采樣型濾波器結構，其中r=0.999，H(k)取第3步中的值。為了簡化編程，梳妝濾波器可以調用CombFilter.m，諧振器可以調用 Resonator2.m，使用help CombFilter 和help Resonator2 查看如何配置參數(shù)。將第1 步生成的采樣信號通過該濾波器，畫出輸出信號第二個周期（n = N, N +1, L1）的時域波形和頻譜，并與第2 步的頻譜

4、進行對比，觀察并分析二者的區(qū)別。信號通過濾波器后，幅值發(fā)生期望變化，相頻特性圖存在一定“異?！?，但是經(jīng)過分析后發(fā)現(xiàn)，是由于算法中用到幅值作為分母，所以對應幅值為零的地方出現(xiàn)異常。5.（選做）分別畫出圖1 中前4 路諧振器的輸出信號第二個周期（n = N, N +1,L1）的時域波形，觀察并分析輸出信號的特點。第一路為直流分量，第二路第三路諧振系數(shù)決定波形如圖，第四路剛好恒為零。7. 思考并回答下列問題（1） 在第2 步的幅頻特性中，各次諧波的幅度與相應的時域信號幅度系？關答：（2） 實驗中為什么要觀察第二個周期，如果直接觀察第一個周期會怎么樣？答：因為該濾波器工作原理要求有一個周期的延時，所以

5、只有觀察第二個周期才能看到正常的現(xiàn)象。（3） 如果取r=0.95，觀察會出現(xiàn)什么情況。答：濾波器性能將變差。附錄：實驗代碼%繪制出采樣時刻從 0 到L ?1 的采樣信號波形A=0.5 1 0.5 2;Fi=0 pi/2 pi -pi/2; f0=50;syms t;S=0;N=16;for k=1:4S=S+A(k)*cos(2*pi*(k-1)*f0*t+Fi(k);endNum=0:2*N-1;t=Num/(N*f0);S=eval(S);subplot(3,3,1) hold on stem(Num/(16*f0),S ezplot(S,0,2/f0) title(采樣圖),r)%對采樣

6、信號的第二個周期（n = N, N +1,?, L%出幅頻特性和相頻特性圖?1）進行離散變換，畫FouS=fft(S(N+1):(2*N);subplot(3,3,2)stem(0:N-1,abs(FouS),g) title(幅頻特性) subplot(3,3,3)stem(0:N-1,angle(FouS),r) title(相頻特性)%濾波器抽頭系數(shù) h(n)% syms w w=0:0.01:2*piHk=1 exp(-1i*pi*(N-1)/N) exp(-2i*pi*(N-1)/N)-exp(-14i*pi*(N-1)/N) -exp(-15i*pi*(N-1)/N);0 0 0

7、0 0 0 0 0 0 0 0h=ifft();subplot(3,3,4) stem(0:N-1,h) title(抽頭系數(shù)) for k=1:Nif k=1H_e_jw=Hk(k).*sin(w-2.*pi.*(k-1)./N).*N./2).*exp(-1i.*(w-2.*pi.*(k-1)./N).*(N-1)./2)./(N.*sin(w-2.*pi.*(k-1)./N)./2); elseH_e_jw=H_e_jw+Hk(k).*sin(w-2.*pi.*(k-1)./N).*N./2).*exp(-1i.*(w-2.*pi.*(k-1)./N).*(N-1)./2)./(N.*s

8、in(w-2.*pi.*(k-1)./N)./2); endend subplot(3,3,5)% ezplot(abs(H_e_jw),0,2*pi) plot(w,abs(H_e_jw)title(濾波器幅頻特性) subplot(3,3,6)% w=0:0.02:2*pi;plot(w,angle(H_e_jw),-)% ezplot(angle(H_e_jw),0,2*pi) title(濾波器相頻特性)%實現(xiàn)圖 1 所示的頻率采樣型濾波器結構r=0.999; temp=0;bFilter(S for k=0:N/2-1,N,r);temp=temp+Resonator2(y1,N,r

9、,k,Hk(k+1);endy=temp/N; y=y(N+1:2*N); subplot(3,3,7) stem(0:N-1,y)title(第二周期時域波形)temp=fft(y); subplot(3,3,8) stem(0:N-1,abs(temp) title(幅頻特性) subplot(3,3,9)stem(0:N-1,angle(temp),r) title(相頻特性)%繪制出采樣時刻從 0 到L ?1 的采樣信號波形A=0.5 1 0.5 2;Fi=0 pi/2 pi -pi/2; f0=50;syms t;S=0;N=16;for k=1:4S=S+A(k)*cos(2*pi

10、*(k-1)*f0*t+Fi(k);endNum=0:2*N-1;t=Num/(N*f0);S=eval(S);syms wHk=1 exp(-1i*pi*(N-1)/N) exp(-2i*pi*(N-1)/N)-exp(-14i*pi*(N-1)/N) -exp(-15i*pi*(N-1)/N);0 0 0 0 0 0 0 0 0 0 0for k=1:Nif k=1H_e_jw=Hk(k)*sin(w-2*pi*(k-1)/N)*N/2)*exp(-1i*(w-2*pi*(k-1)/N)*(N-1)/2)/(N*sin(w-2*pi*(k-1)/N)/2); elseH_e_jw=H_e_

11、jw+Hk(k)*sin(w-2*pi*(k-1)/N)*N/2)*exp(-1i*(w-2*pi*(k-1)/N)*(N-1)/2)/(N*sin(w-2*pi*(k-1)/N)/2); endendx=S; r=0.999L=length(x); y=zeros(1,L+N);if L=Ny(1:N)=x(1:N);for n=N+1:Ly(n)=x(n)-x(n-N)*rN;endy(L+1:L+N)=-x(L+1-N:L)*rN;elsey(1:L)=x(1:L);for n=N+1:N+Ly(n)=-x(n-N)*rN;endendx=y; forOrder=0:3 H=Hk(Ord

12、er+1); if OrderN/2display(Order should not be greater Resonator2!);elseL=length(x); f=zeros(1,L);y=zeros(1,L);nN/2infunctionifmod(N,2)=0if Order=0y(1)=x(1)*H;for n=2:Ly(n)=H*x(n)+y(n-1)*r;endelseif Order=N/2 y(1)=x(1)*H;for n=2:Ly(n)=H*x(n)-y(n-1)*r;end elsew=exp(1i*2*pi*Order/N); alfa0=2*real(H); a

13、lfa1=-2*r*real(H.*w);beta1=2*r*cos(2*pi*Order/N); beta2=-r2;f(1)=x(1); y(1)=f(1)*alfa0;f(2)=x(2)+f(1)*beta1; y(2)=f(2)*alfa0+f(1)*alfa1; for n=3:Lf(n)=x(n)+f(n-1)*beta1-f(n-2)*r2; y(n)=f(n)*alfa0+f(n-1)*alfa1;endendelseif Order=0y(1)=x(1)*H;for n=2:Ly(n)=H*x(n)+y(n-1)*w*r;end elsew=exp(1i*2*pi*Order

14、/N); alfa0=2*real(H); alfa1=-2*r*real(H.*w);beta1=2*r*cos(2*pi*Order/N); beta2=-r2;f(1)=x(1); y(1)=f(1)*alfa0;f(2)=x(2)+f(1)*beta1; y(2)=f(2)*alfa0+f(1)*alfa1; for n=3:Lf(n)=x(n)+f(n-1)*beta1-f(n-2)*r2; y(n)=f(n)*alfa0+f(n-1)*alfa1;endendendendtemp=y/N; temp=temp(N+1:2*N); subplot(2,2,Order+1) stem(

15、0:N-1,temp)end% CombFilter.m - by pcmu%Description functionbFilter(x,N)Comb Parametersx:N:r: Return y:filterInput signal Order of the comb RadiusfilterOutput signalfunctionbFilter(x,N,r)L=length(x); y=zeros(1,L+N); if L=Ny(1:N)=x(1:N);for n=N+1:Ly(n)=x(n)-x(n-N)*rN;endy(L+1:L+N)=-x(L+1-N:L)*rN;elsey

16、(1:L)=x(1:L);for n=N+1:N+Ly(n)=-x(n-N)*rN;endend%Resonator2.m - by pcmuDescriptionfunction y=Resonator2(x,N,r,Order,H)Resonator of order 2 Parametersto obtain real-valuedharmonicsx:N:r: Order: H:Return y:Input signal Order of the RadiusOrder of thecomb filterresonatorMagnification coeficient of the

17、resonatorOutput signalfunctiony=Resonator2(x,N,r,Order,H)if OrderN/2display(Order Resonator2!); elseL=length(x); f=zeros(1,L);y=zeros(1,L); if mod(N,2)=0shouldnotbegreaternN/2infunctionif Order=0y(1)=x(1)*H;for n=2:Ly(n)=H*x(n)+y(n-1)*r;endelseif Order=N/2 y(1)=x(1)*H;for n=2:Ly(n)=H*x(n)-y(n-1)*r;e

18、ndelsew=exp(1i*2*pi*Order/N); alfa0=2*real(H);alfa1=-2*r*real(H.*w); beta1=2*r*cos(2*pi*Order/N); beta2=-r2;f(1)=x(1); y(1)=f(1)*alfa0;f(2)=x(2)+f(1)*beta1; y(2)=f(2)*alfa0+f(1)*alfa1; for n=3:Lf(n)=x(n)+f(n-1)*beta1-f(n-2)*r2; y(n)=f(n)*alfa0+f(n-1)*alfa1;endend elseif Order=0y(1)=x(1)*H;for n=2:Ly(n)=H*x(