1、數字信號處理實驗報告姓名：_賀云天_學號：_2010035015_班級：_電信碩23_實驗日期：2014年10月22日提交日期：2014年10月30日實驗二 頻率采樣型濾波器一、 實驗目的1. 學習使用頻率采樣型結構實現FIR濾波器，初步熟悉FIR濾波器的線性相位特點；2. 直觀體會頻率采樣型濾波器所具有的“濾波器組”特性，即在并聯結構的每條支路上可以分別得到輸入信號的各次諧波；3. 學習使用周期沖激串檢測所實現濾波器的頻域響應。二、 實驗內容頻率采樣型濾波器是由一個梳狀濾波器和若干路諧振器構成的，可用公式表述如下：Hz=1-rNz-NNk=0N-1Hk1-rWN-kz-1(1)其中r值理論上

2、為1，實際中取非常接近1的值。為了使系數為實數，可以將諧振器的共軛復根合并，不失一般性，假設N為偶數，于是可以得到如圖1所示的結構。圖1 N為偶數的實系數頻率采樣型結構濾波器其中0k=2ReH(k)，1k=-2rReH(k)WNk以下實驗中假設頻率采樣型濾波器階數N=16。1. 構造濾波器輸入信號st=k=03sk(t)，其中skt=Akcos(2kf0t+k)，基波頻率f0=50Hz，A0=0.5，A1=1，A2=0.5，A3=2，0=0，1=2，0=，3=-2。設時域信號s(t)的采樣頻率fs=Nf0，繪制出采樣時刻從0到L-1的采樣信號波形，其中采樣點數為L=2N，確認時域信號采樣正確。

3、2. 對采樣信號的第二個周期（n=N,N+1,L-1）進行離散傅里葉變換，畫出幅頻特性和相頻特性圖，觀察并分析其特點。3. 設H(0)=1，H(1)=e-j(N-1)N，H(2)=e-j2(N-1)N，H3=H4=H(13)=0，H14=-e-j14(N-1)N，H(15)=-e-j15(N-1)N，計算濾波器抽頭系數h(n)，n=0,1,N-1，畫出該濾波器的頻譜圖，觀察并分析其幅頻特性和相頻特性。4. 編程實現圖 1 所示的頻率采樣型濾波器結構，其中r=0.999，H(k)取第3步中的值。為了簡化編程，梳妝濾波器可以調用CombFilter.m，諧振器可以調用 Resonator2.m，使

4、用 help CombFilter和help Resonator2查看如何配置參數。將第1步生成的采樣信號通過該濾波器，畫出輸出信號第二個周期（n=N,N+1,L-1）的時域波形和頻譜，并與第2步的頻譜進行對比，觀察并分析二者的區別。5. 分別畫出圖1中前4路諧振器的輸出信號第二個周期（n=N,N+1,L-1）的時域波形，觀察并分析輸出信號的特點。6. 將輸入信號換成周期為N的沖激串，畫出輸出信號第二個周期（n=N,N+1,L-1）的幅頻特性，并與第3步的濾波器幅頻特性進行對比，觀察并分析二者的關系。7. 思考題（1）在第2步的幅頻特性中，各次諧波的幅度與相應的時域信號幅度有什么關系？（2）實

5、驗中為什么要觀察第二個周期，如果直接觀察第一個周期會怎么樣？（3）如果取r=0.95，觀察會出現什么情況。三、 實驗要求1. 按照實驗內容編寫Matlab程序，給出運行結果（圖），并逐項進行分析討論。2. 提交完整的Matlab源程序。3. 回答思考題，撰寫實驗報告。四、 實驗結果與分析1. 構造濾波器輸入信號st=k=03sk(t)，其中skt=Akcos(2kf0t+k)，基波頻率f0=50Hz，A0=0.5，A1=1，A2=0.5，A3=2，0=0，1=2，0=，3=-2。設時域信號s(t)的采樣頻率fs=Nf0，繪制出采樣時刻從0到L-1的采樣信號波形，其中采樣點數為L=2N，確認時域

6、信號采樣正確。輸入信號是直流、一次波、二次波、三次波的疊加，采樣頻率為16f0，高于輸入信號的2倍。輸入信號的時域圖像如上。2. 對采樣信號的第二個周期（n=N,N+1,L-1）進行離散傅里葉變換，畫出幅頻特性和相頻特性圖，觀察并分析其特點。由幅頻特性圖可知，在第0、1、2、3、13、14、15個點上有幅值，數值對稱、相位反對稱，其余點的幅度為0。第0點表示基波分量；第1點和第15點為直流分量，其幅值和是基波分量的2倍；第2點和第14點為二次分量，其幅值和與基波分量相等；第3點和第13點為三次分量，幅值和是基波分量的4倍。同時，從FFT結果可發現，幅值為0的分量相位并不相同，這是由計算機計算精

7、度造成的。3. 設H(0)=1，H(1)=e-j(N-1)N，H(2)=e-j2(N-1)N，H3=H4=H(13)=0，H14=-e-j14(N-1)N，H(15)=-e-j15(N-1)N，計算濾波器抽頭系數h(n)，n=0,1,N-1，畫出該濾波器的頻譜圖，觀察并分析其幅頻特性和相頻特性。濾波器抽頭系數h(n)可由H(k)作IFFT運算得出，第一行兩幅圖表示出H(k)的幅頻特性和相頻特性。將H(k)內插后觀察，發現H2為低通濾波器。4. 編程實現圖 1 所示的頻率采樣型濾波器結構，其中r=0.999，H(k)取第3步中的值。為了簡化編程，梳妝濾波器可以調用CombFilter.m，諧振器

8、可以調用 Resonator2.m，使用 help CombFilter和help Resonator2查看如何配置參數。將第1步生成的采樣信號通過該濾波器，畫出輸出信號第二個周期（n=N,N+1,L-1）的時域波形和頻譜，并與第2步的頻譜進行對比，觀察并分析二者的區別。第四題輸出圖像：信號通過濾波器后第二個周期的時域圖像、FFT后的頻譜信號直接作FFT后的頻譜幅值：與第2題對比，第0、1、2、14、15點不變，第3點和第13點的數值為0，這是因為三次分量在通過濾波器時被濾掉了，而濾波器在第0、1、2、14、15點的幅值為1，所以通過濾波器后的幅值沒有發生變化。相位：直流分量的初始相位與濾波器

9、H(0)相乘，得到新的相位；一次分量的相位由基波初始相位0與H(1)相位相乘而得；其他分量的相位無法確定。5. 分別畫出圖1中前4路諧振器的輸出信號第二個周期（n=N,N+1,L-1）的時域波形，觀察并分析輸出信號的特點。H(0)輸出波形為直流分量，所以波形的幅度不變；因為H(3)=0，所以H(3)輸出波形也為0，幅度不變；其他兩路分別為基波和二次分量的波形。6. 將輸入信號換成周期為N的沖激串，畫出輸出信號第二個周期（n=N,N+1,L-1）的幅頻特性，并與第3步的濾波器幅頻特性進行對比，觀察并分析二者的關系。輸入信號是周期為N的沖激串時，輸出信號的時域圖像和頻譜濾波器幅頻特性沖激函數在頻域

10、上為1，通過濾波器后被濾掉k=313區域，所以其幅值和濾波器本身的幅值相等，保留了直流、基波、二次分量。五、 思考題1. 在第2 步的幅頻特性中，各次諧波的幅度與相應的時域信號幅度有什么關系？頻率為0時的幅度等于時域波形中直流的幅度，直流、二次、三次諧波的譜線的幅度為時域波形相應諧波對應幅度的一半。因為正弦信號的頻譜中有 和兩條譜線，幅度都為信號的一半，cos(t)=ejw+e-jw2。2. 實驗中為什么要觀察第二個周期，如果直接觀察第一個周期會怎么樣？第二個周期信號與系統的單位脈沖響應卷積后，輸出信號開始變得周期，不會發生失真；直接觀察第一個周期，信號沒有完全進入系統，卷積結果不完整，輸出會

11、發生失真。3. 如果取r=0.95，觀察會出現什么情況。所有諧振器的極點從單位圓向內收縮一點，同時梳狀濾波器的零點也移到r圓上，濾波器的幅頻特性的幅度減小，通帶內的仍然保持線性相位特性。4. 如何理解第3 步與第6 步在工程使用中的區別？兩者均為低通濾波器，但是第三步是進行頻率采樣后由內插函數得到傳遞函數H(z)，從而得到系統頻響H()；而題6中使用了沖激串，這是因為在理想條件下，沖激串的頻響就是系統的頻響，但是在實際生活中很難實現。六、 MATLAB程序代碼1. 以下是MATLAB程序代碼：clc;clear;clf;f0=50;N=16;fs=N*f0; %時域信號s(t)的采樣頻率T=1

12、/fs;%采樣間隔L=N*2; %采樣點數t=0:T:(L-1)*T;%構造輸入信號s0=0.5*cos(0);s1=1*cos(2*pi*f0*t+pi/2);s2=0.5*cos(2*pi*2*f0*t+pi);s3=2*cos(2*pi*3*f0*t-pi/2);s=s0+s1+s2+s3;i=0:1:L-1;stem(i,s(i+1);grid on;2. 以下是MATLAB程序代碼：clc;clear;clf;f0=50;N=16;fs=N*f0; %時域信號s(t)的采樣頻率T=1/fs;%采樣間隔L=N*2; %采樣點數t=N*T:T:(L-1)*T;s0=0.5*cos(2*p

13、i*0*f0*t+0);s1=1*cos(2*pi*f0*t+pi/2);s2=0.5*cos(2*pi*2*f0*t+pi);s3=2*cos(2*pi*3*f0*t-pi/2);s=s0+s1+s2+s3;a=fft(s);%FFT運算a%輸出a的各項值y=abs(a);z=angle(a);i=0:1:N-1;subplot(2,1,1);stem(i,y(i+1);subplot(2,1,2);stem(i,z(i+1);3. 以下是MATLAB程序代碼；clear;clc;clf;N=16;H=1,exp(-1i*pi*(N-1)/N),exp(-j*2*pi*(N-1)/N),0,

14、0,0,0,0,0,0,0,0,0,0,-exp(-j*14*pi*(N-1)/N),-exp(-j*15*pi*(N-1)/N);h=ifft(H,N);hH_abs=abs(H);H_angle=angle(H);i=0:1:N-1;subplot(2,2,1);stem(i,H_abs(i+1);title(H 幅頻特性);grid on;subplot(2,2,2);stem(i,H_angle(i+1);title(H 相頻特性);grid on;w=0:pi/200:2*pi;H2=zeros(1,length(w);for i=1:length(w) for n=1:N e=ex

15、p(-j*w(i)*(n-1); H2(i)=H2(i)+h(n)*e; endendw=0:pi/200:2*pi;subplot(2,2,3);plot(abs(H2);title(H2 幅頻特性);grid on;subplot(2,2,4);plot(angle(H2);title(H2 相頻特性);grid on;4. 以下是MATLAB程序代碼；clc;%初始化數據f0=50;N=16;fs=N*f0;L=2*N;T=1/fs;t=0:T:(L-1)*T;s0=0.5*cos(0);s1=1*cos(2*pi*1*f0*t+pi/2);s2=0.5*cos(2*pi*2*f0*t+

16、pi);s3=2*cos(2*pi*3*f0*t-pi/2);s=s0+s1+s2+s3;x=s;r=0.999;y=CombFilter(x,N,r);x=y;z=zeros(1,48);for i=0:1:(N/2)Order=i;H=1 exp(-j*pi*(N-1)/N) exp(-j*2*pi*(N-1)/N) 0 0 0 0 0 0 0 0 0 0 0 -exp(-j*14*pi*(N-1)/N) -exp(-j*15*pi*(N-1)/N);y=Resonator2(x,N,r,Order,H(i+1); z=z+y;endy=z/N; for i=N:1:L-1; y1(i-N

17、+1)=y(i+1); %取出N到L-1的元素end%作時域圖i=0:1:N-1;subplot(2,2,1),stem(i,y1(i+1);xlabel(n);ylabel(振幅);title(時域波形);grid on; %作頻譜y=fft(y1);yy1=abs(y);ang=angle(y);i=0:1:N-1;subplot(2,2,3),stem(i,y1(i+1);xlabel(n);ylabel(幅值);title(幅頻特性);grid on;subplot(2,2,4),stem(i,ang(i+1);xlabel(n);ylabel(幅值);title(相頻特性);grid

18、 on;5. 以下是MATLAB程序代碼；clc;f0=50;N=16;fs=N*f0; T=1/fs;L=N*2; t=0:T:(L-1)*Ts0=0.5*cos(0);s1=1*cos(2*pi*f0*t+pi/2);s2=0.5*cos(2*pi*2*f0*t+pi);s3=2*cos(2*pi*3*f0*t-pi/2);s=s0+s1+s2+s3;%生成信號%前4路諧振器的輸出r=0.999;y1=CombFilter(s,N,r);H=1,exp(-1i*pi*(N-1)/N),exp(-1i*2*pi*(N-1)/N),0,0,0,0,0,0,0,0,0,0,0,exp(-1i*1

19、4*pi*(N-1)/N),exp(-1i*15*pi*(N-1)/N);y_h0=Resonator2(y1,N,r,0,H(1);y_h0=y_h0/N;y_h0=y_h0(N:2*N-1);y_h1=Resonator2(y1,N,r,1,H(2);y_h1=y_h1/N;y_h1=y_h1(N:2*N-1);y_h2=Resonator2(y1,N,r,2,H(3);y_h2=y_h2/N;y_h2=y_h2(N:2*N-1);y_h3=Resonator2(y1,N,r,3,H(4);y_h3=y_h3/N;y_h3=y_h3(N:2*N-1);%輸出圖像k=0:N-1;subplot(2,2,1);stem(k,y_h0);title(H(0)路輸出信號第二個周期);subplot(2,2,2);stem(k,y_h1);title(H(1)路輸出信號第二個周期);subplot(2,2,3);stem(k,y_h