1、華 北 電 力 大 學 實 驗 報 告華 北 電 力 大 學 實 驗 報 告實驗環境MATLAB 6.5實驗名稱實驗一： FFT的應用實 驗 目 的1、熟悉MATLAB在數字信號處理中的應用。2、掌握利用FFT計算序列線性卷積的基本原理及編程實現。3、掌握對連續信號進行采樣的基本原理和方法，并利用FFT對信號進行頻譜分析。實 驗 原 理1、 離散時間信號的表示：在數字信號處理中，所有信號都是離散時間信號序列，所有的序列都可以表示為如下形式：x(n)=,x(-1),x(0),x(1), 用matlab的stem（n,x）即可實現簡單的表示。 若要表示具有特定采樣頻率的信號，需要定義時間軸向量。2

2、、 快速傅里葉變換：采用時間抽取基2FFT算法。MATLAB中提供了fft和ifft函數來分別計算DFT和IDFT。fft和ifft函數是用機器語言，而不是用MATLAB指令寫成的，因此它的執行速度很快。 fft函數的用法： y=fft(x);%計算x的快速離散傅里葉變換y y=fft(x,N);%計算x的N點FFT。當x的長度大于N時，截斷x；否則補零 ifft函數的用法： y=ifft(x);%計算x的快速離散傅里葉反變換y y=ifft(x,N);%計算x的N點IFFT 3、 利用FFT計算線性卷積：步驟如下：1) 將x1 (n)和x2(n)都延長到N點， N=N1+N2-12) 計算x

3、1(n)的N點DFT，即：X1(k)=FFTx1(n) 3) 計算x2(n)的N點DFT,即：X2(k)=FFTx2(n)4) 計算 Y（k）=X1(k)X2(k)5) 計算的反變換,即y(n)=IFFTX1(k)X2(k)4、 利用FFT對信號進行譜分析：包括振幅譜，相位譜和功率譜。參數：fs=2f0在matlab中，提供了計算模值的函數abs和計算相角的函數angle，其用法如下：Abs(x);Angle(x);實 驗 內 容1.對于兩個序列：x(n)=nR16(n),h(n)=R8(n)（1）在同一圖形窗口中繪出兩序列的時域圖形。（2）利用FFT編程計算兩序列的線性卷積，繪出的時域圖形。

4、設計方案： 首先定義x序列，為離散序列，再通過stem函數將x依次與y序列對應，即可表示兩序列。然后利用FFT計算倆序列的線性卷積，方法與步驟如實驗原理所述，注意為使計算快速，可延長時把兩序列都延長到最接近的2的整數冪的值，通過函數N=2nextpow2(N1+N2-1)實現，然后兩序列補零，求傅里葉變換，相乘，求反變換，依次即可實現。 源程序：%兩序列時域波形N1=16;x1=0:N1-1;y1=x1;N2=8;x2=0:N2-1;y2=ones(1,N2);subplot(2,1,1);stem(x1,y1,b);hold on;stem(x2,y2,r);title(兩序列時域波形);%

5、線性卷積時域波形N=2nextpow2(N1+N2-1);y3=y1(1:N1) zeros(1,N-N1);y4=y2(1:N2) zeros(1,N-N2);Y1=fft(y3,N);Y2=fft(y4,N);Y=Y1.*Y2;y=ifft(Y,N);n=0:N-1;subplot(2,1,2);stem(n,y);title(線性卷積時域波形);2利用FFT對信號進行譜分析 對于連續信號xa(t)=cos(2f1t) +5cos(2f2t) +cos(2f3t) ，其中f1=6.5kHz, f2=7kHz, f3=9kHz, 以采樣頻率fs=32 kHz對其進行采樣，（1）對xa(t)

6、信號采集16點樣本，分別作16點和補零到256點的FFT，并分別繪出對應的幅頻特性曲線。（2）對xa(t)信號采集256點樣本，分別作256點和512點的FFT，并分別繪出對應的幅頻特性曲線。 （3）比較（1）和（2）中的結果，分析采樣點數和傅里葉變換點數對FFT的影響，說明高密度頻譜和高分辨率頻譜的特點與區別。設計方案：若要表示具有特定采樣頻率的信號，先需要根據采樣頻率定義時間軸向量，然后定義函數，依次表示出16個采樣點，補零到256個采樣點，256個采樣點和補零至512個采樣點時的時域波形和頻譜，圖形表示時直接使用stem（）函數即可，補零用x2=x(1:N1) zeros(1,N2-N1

7、)這樣的形式，快速傅里葉變換為fft()函數，幅值用abs()函數直接求解。源程序：f1=6500;f2=7000;f3=9000;fs=32000;t=0:1/fs:1;x=cos(2*f1*pi*t)+5*cos(2*f2*pi*t)+cos(2*f3*pi*t);N1=16;N2=256;N3=512;n=N3-1;%16個采樣點的時域波形，頻譜n1=0:N1-1;x1=x(1:N1);subplot(2,4,1);stem(n1,x1);title(16個采樣點的時域信號x1);Y1=fft(x1,N1);Y1=abs(Y1);k1=0:N1-1;f1=fs/N1*k1;subplot

8、(2,4,5);stem(f1,Y1);title(信號x1的頻譜);%補零到256個采樣點的時域波形，頻譜n2=0:N1-1;x2=x(1:N1) zeros(1,N2-N1);m1=0:255;subplot(2,4,2);stem(m1,x2);title(補零到256點的時域信號x2);Y2=fft(x2,N2);Y2=abs(Y2);k2=0:N2-1;f2=fs/N2*k2;subplot(2,4,6);stem(f2,Y2);title(信號x2的頻譜);% 256個采樣點的時域波形，頻譜n3=0:N2-1;x3=x(1:N2);subplot(2,4,3);stem(n3,x3

9、);title(256個采樣點的時域信號x3);Y3=fft(x,N2);Y3=abs(Y3);k3=0:N2-1;f3=fs/N2*k3;subplot(2,4,7);stem(f3,Y3);title(信號x3的頻譜);%補零到512個采樣點的時域波形，頻譜n4=0:N2-1;x4=x(1:N2) zeros(1,N3-N2);m2=0:511;subplot(2,4,4);stem(m2,x4);title(補零到512點的時域信號x4);Y4=fft(x4,N3);Y4=abs(Y4);k4=0:N3-1;f4=fs/N3*k4;subplot(2,4,8);stem(f4,Y4);t

10、itle(信號x4的頻譜);實 驗 結 果 及 分 析1.對于兩個序列：x(n)=nR16(n),h(n)=R8(n)（1）在同一圖形窗口中繪出兩序列的時域圖形。（2）利用FFT編程計算兩序列的線性卷積，繪出的時域圖形。實驗結果：實驗分析：兩信號分別為x(n)=nR16(n),h(n)=R8(n)，為時間離散信號，兩個信號同時在第一個圖形中表示出來，結果與信號表達式相同，結果正確。兩信號卷積x(n)*h(n)為所有h(n-k)與x(k)的乘積的累加和。所以在n=15時會出現卷積的最大值。X(n)和h(n)的線性卷積的結果y(n)是一個有限序列，其非零值長度為N1+N2-1,最終matlab實現

11、的結果為22個非零值，實驗結果正確。2利用FFT對信號進行譜分析實驗結果：實驗分析：程序運行結果如上圖所示，可以得到如下結論：當取樣點N=16時，從頻譜圖上幾乎無法看出信號的任何頻率信息。將16點中的信號補零到N=256點時，頻譜的譜線相當密，但從中仍然很難看出信號的頻率成分，稱之為高密度頻譜。它只是在信號X1的頻譜基礎上增加采樣密度，但不增加分辨率，無法提取有用的頻譜成分。當對序列取足夠的采樣點時，如256個采樣點，可以從幅頻特性中清晰地看出信號的頻率成分，稱之為高分辨率頻譜。頻譜應為6.5kHz，7kHz，9kHz。另一組為23kHz,25kHz,25.5kHz。（如6.5k對稱的為-6.

12、5k采樣頻率為32k,32-6.5為25.5，故第二組為含25.5kHz的頻譜，其余同理即可得，如上面第二三個圖所示）問題：高密度頻譜與高分辨率頻譜有哪些特點和區別：答：高密度頻譜與高分辨率頻譜的譜線都非常密集，采樣頻率都非常高。但高密度頻譜是當時間域長度不變時，頻域內對它的頻譜進行提高采樣頻率，細化當前分辨率下的頻譜，故對分辨率沒有改善。而高分辨率頻譜才是我們所希望得到的頻譜結果。實驗名稱設計性實驗一：IIR數字濾波器的設計實 驗 目 的1、本實驗為設計性實驗。2、掌握用雙線性變換法設計IIR數字濾波器的基本原理和設計方法。3、掌握用雙線性變換法設計IIR數字Butterworth濾波器的原

13、理和設計方法。實 驗 原 理1、IIR數字濾波器設計的原理：IIR數字濾波器的系統函數H（z）是z-1的有理分式，設計IIR數字濾波器，就是要確定系統函數H(z)的階數N以及分子分母多項式的系數ak和bk，使其頻率特性滿足指定的要求。用雙線性變化法設計IIR數字濾波器的步驟如下：（1）確定數字濾波器的通帶截止頻率wp、阻帶截止頻率ws、通帶最大衰減Rp和阻帶最小衰減Rs；（2）對截止頻率wp和ws進行頻率預畸變；（3）計算模擬濾波器的階數N和頻率Wc，進而求得s平面上的極點，從而構成模擬濾波器的傳遞函數Ha(s)，使得滿足給定的技術指標；（4）利用雙線性變換公式將模擬濾波器的傳遞函數Ha(s)

14、變換成數字濾波器的系統函數H(z)；（5） 若需設計其他形式的濾波器，則根據頻率變換公式將數字低通濾波器原型變換到相應形式的濾波器即可。模擬域頻率與數字域頻率的關系為Wc=（2/T）tan(w/2)這是一種非線性的關系。這種非線性關系使得模擬濾波器和數字濾波器的頻率響應與對應頻率的關系上發生了畸變，也造成了相位的非線性變化。為保證各邊界頻率點為預先指定的頻率，在確定模擬低通濾波器的系統函數之前，必須按上式進行頻率預畸變。2、Matlab中應用到的函數：在MATLAB的數字信號處理工具箱中，提供了用雙線性變換法設計數字Butterworth低通濾波器的相關函數。（1）buttord：Butter

15、worth濾波器階數選擇函數N,Wn=buttord(Wp,Ws,Rp,Rs, s);其中，Wp為通帶截止角頻率，Ws為阻帶截止角頻率，Rp為通帶最大衰減，Rs為阻帶最小衰減；N是符合要求的濾波器最小階數，Wn是Butterworth濾波器固有角頻率（3dB）。Wp、Ws、Wn均歸一化在之間，單位為 rad/s；s表示用于模擬濾波器，去掉則用于數字濾波器。 （2）buttap：Butterworth模擬低通濾波器的建立函數Z,P,K=buttap(N);給出N階Butterworth模擬濾波器的零點向量Z、極點向量P和增益K。產生的濾波器在左半平面的單位圓附近有N個極點，沒有零點。（3）zp2

16、tf：零極點增益模型到傳遞函數模型的轉換函數Bap,Aap=zp2tf(Z,P,K);其中，Z、P、K分別為零極點增益模型的零點向量、極點向量和增益；Bap、Aap分別為傳遞函數分子和分母的多項式系數向量。（4）lp2lp：從低通原型向低通的轉換函數b,a=lp2lp(Bap,Aap,Wn);把截止角頻率為1 rad/s的模擬低通原型濾波器轉換成截止角頻率為Wn的模擬低通濾波器。（5）bilinear：雙線性變換函數bz,az=bilinear(b,a,Fs);用雙線性變換法把模擬低通濾波器轉換為數字低通濾波器。其中，Fs是采樣頻率，bz和az分別是傳遞函數的分子和分母的多項式系數向量。（6）

17、freqz：數字濾波器的頻率響應函數H,W=freqz(B,A,N);返回數字濾波器均勻分布在上的N點頻率向量W和與之對應的N點頻率響應向量H。A和B分別是濾波器系統函數分子和分母的多項式系數向量；N最好選用2的整數冪，以便使用FFT快速運算，N的缺省值為512。H,F=freqz(B,A,N,Fs)對在0,Fs/2上等間隔采樣N點，采樣點頻率及相應的頻率響應值分別記錄在F和H中。實 驗 內 容IIR數字濾波器的設計 用雙線性變換法設計一個IIR數字Butterworth低通濾波器。技術指標為：通帶截止頻率fp=1kHz ，阻帶截止頻率fs=1.5kHz ，通帶衰減Rp1dB，阻帶衰減Rs 4

18、0dB ，采樣頻率Fs=10kHz。繪出濾波器的幅頻特性曲線和相頻特性曲線，判斷設計是否符合要求。設計方案：按照雙線性變換法設計IIR濾波器的步驟，首先指標轉換，w=2pi*f，因為后面用到buttord()函數，所以指標轉換時歸一化處理。對wp,ws進行預畸變，將數字域頻率轉化成模擬域頻率，利用buttord函數求濾波器階數，通過函數由模擬濾波器的零點向量，極點向量，增益求模擬低通的傳遞函數。再通過雙線性變換函數bilinear()將模擬低通轉換為數字低通。通過freqz函數求數字濾波器的頻率響應函數即可。源程序：clear all;fp=1000;fs=1500;Fs=10000;Rp=1

19、.5;Rs=25;%指標轉換wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Fs=Fs/Fs;%預畸變Wp=2*tan(wp/2); Ws=2*tan(ws/2);%由階數，零點向量，極點向量，增益實現模擬低通的傳遞函數N,Wn=buttord(Wp,Ws,Rp,Rs,s); z,p,k=buttap(N); Bap,Aap=zp2tf(z,p,k); %由低通原型轉向低通b,a=lp2lp(Bap,Aap,Wn); %雙線性變換實現模擬低通變為數字低通bz,az=bilinear(b,a,Fs);%頻率響應h,w=freqz(bz,az,512,Fs); %輸出幅頻h1=abs(h)

20、;subplot(2,1,1);h2=20*log10(h1+eps)/max(h1);plot(w*10000,h2);title(幅頻特性);grid on; %輸出相頻subplot(2,1,2);plot(w*10000,angle(h);title(相頻特性);grid on;實 驗 結 果 及 分 析實驗結果：實驗分析：上圖為，程序運行后濾波器的幅度特性和相頻特性。從圖中可以看出：幅度響應曲線在通帶截止頻率（f=1000hz）即在通帶內有更大部分的幅度接近于0，可以從圖中看出，通帶截止頻率為1000Hz，阻帶截止頻率1500Hz為下降3dB是的頻率值，縱坐標換成dB時可清楚看到。幅

21、頻特性圖中可以清晰的看到數字濾波器的通帶、過渡帶、阻帶三部分。IIR數字濾波器不能做到嚴格的線性相位。 所以結果符合設計要求。問題：沖激響應不變法和雙線性變換法比較： 答：相同點：首先設計模擬濾波器，再將模擬濾波器轉換為數字濾波器。不同點：沖激響應不變法的幅頻響應有失真，頻率之間呈線性。 雙線性變換法的幅頻響應無混迭 ，頻率之間有失真。實驗名稱設計性實驗二：FIR數字濾波器的設計實 驗 目 的1、本實驗為設計性實驗。2、掌握用窗函數法設計IIR數字濾波器的基本原理和設計方法。3、掌握用窗函數法設計線性相位FIR數字低通濾波器的編程實現。實 驗 原 理1、設計原理：FIR數字濾波器的獨特優點是容

22、易得到嚴格的線性相位，所以FIR總是最穩定的。FIR數字濾波器主要采用非遞歸結構，而且可以借助FFT計算。因為所要設計的是FTR數字濾波器，其單位沖激響應h(n)必然是有限長的，所以要用有限長的h(n)來逼近無限長的hd(n)。最有效的方法昰截斷hd(n)，即用有限長的窗函數w(n)來截取hd(n)，表示h(n)= hd(n) w(n),這種方法稱為窗函數法。用窗函數法設計FIR數字濾波器的步驟如下：（1） 給出所要設計的FIR數字濾波器的技術指標，如通帶截止頻率wp、阻帶截止頻率ws、通帶衰減Rp和阻帶衰減As；（2） 根據允許的過渡帶寬度及阻帶衰減，初步選擇窗函數和N值。（3） 若選用理想

23、低通逼近，求出理想低通的沖激響應hd(n)。理想低通的截止頻率選擇為wc=(wp+ws)/2（4） 將hd(n)與窗函數相乘得FIR數字濾波器的沖激響應h(n)（5） 計算FIR數字濾波器的頻率響應，并驗證是否達到所要求的指標。2、matlab中用到的函數：在MATLAB的數字信號處理工具箱中提供了用窗函數法設計線性相位FIR數字濾波器的函數fir1，它具有標準低通、帶通、高通和帶阻等類型。fir1：濾波器設計函數B=fir1(N,Wn, ftype,window);其中，N為FIR濾波器的階數， Wn為濾波器截止頻率，取值范圍【0,1】rad/s，對于帶通、帶阻濾波器，Wn=W1,W2，且W140dB，與實驗的技術指標要求相同。觀察相頻特性，可以看出FIR濾波器有更嚴格的線性相位。通過比較如上兩個窗函數海明窗和凱澤窗，可