數(shù)字信號(hào)處理課程設(shè)計(jì)實(shí)驗(yàn)

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

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

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

4、設(shè)計(jì)方案： 首先定義x序列，為離散序列，再通過(guò)stem函數(shù)將x依次與y序列對(duì)應(yīng)，即可表示兩序列。然后利用FFT計(jì)算倆序列的線性卷積，方法與步驟如實(shí)驗(yàn)原理所述，注意為使計(jì)算快速，可延長(zhǎng)時(shí)把兩序列都延長(zhǎng)到最接近的2的整數(shù)冪的值，通過(guò)函數(shù)N=2nextpow2(N1+N2-1)實(shí)現(xiàn)，然后兩序列補(bǔ)零，求傅里葉變換，相乘，求反變換，依次即可實(shí)現(xiàn)。 源程序：%兩序列時(shí)域波形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(兩序列時(shí)域波形);%

5、線性卷積時(shí)域波形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(線性卷積時(shí)域波形);2利用FFT對(duì)信號(hào)進(jìn)行譜分析 對(duì)于連續(xù)信號(hào)xa(t)=cos(2f1t) +5cos(2f2t) +cos(2f3t) ，其中f1=6.5kHz, f2=7kHz, f3=9kHz, 以采樣頻率fs=32 kHz對(duì)其進(jìn)行采樣，（1）對(duì)xa(t)

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

7、)這樣的形式，快速傅里葉變換為fft()函數(shù)，幅值用abs()函數(shù)直接求解。源程序：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個(gè)采樣點(diǎn)的時(shí)域波形，頻譜n1=0:N1-1;x1=x(1:N1);subplot(2,4,1);stem(n1,x1);title(16個(gè)采樣點(diǎn)的時(shí)域信號(hào)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(信號(hào)x1的頻譜);%補(bǔ)零到256個(gè)采樣點(diǎn)的時(shí)域波形，頻譜n2=0:N1-1;x2=x(1:N1) zeros(1,N2-N1);m1=0:255;subplot(2,4,2);stem(m1,x2);title(補(bǔ)零到256點(diǎn)的時(shí)域信號(hào)x2);Y2=fft(x2,N2);Y2=abs(Y2);k2=0:N2-1;f2=fs/N2*k2;subplot(2,4,6);stem(f2,Y2);title(信號(hào)x2的頻譜);% 256個(gè)采樣點(diǎn)的時(shí)域波形，頻譜n3=0:N2-1;x3=x(1:N2);subplot(2,4,3);stem(n3,x3

9、);title(256個(gè)采樣點(diǎn)的時(shí)域信號(hào)x3);Y3=fft(x,N2);Y3=abs(Y3);k3=0:N2-1;f3=fs/N2*k3;subplot(2,4,7);stem(f3,Y3);title(信號(hào)x3的頻譜);%補(bǔ)零到512個(gè)采樣點(diǎn)的時(shí)域波形，頻譜n4=0:N2-1;x4=x(1:N2) zeros(1,N3-N2);m2=0:511;subplot(2,4,4);stem(m2,x4);title(補(bǔ)零到512點(diǎn)的時(shí)域信號(hào)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(信號(hào)x4的頻譜);實(shí) 驗(yàn) 結(jié) 果 及 分 析1.對(duì)于兩個(gè)序列：x(n)=nR16(n),h(n)=R8(n)（1）在同一圖形窗口中繪出兩序列的時(shí)域圖形。（2）利用FFT編程計(jì)算兩序列的線性卷積，繪出的時(shí)域圖形。實(shí)驗(yàn)結(jié)果：實(shí)驗(yàn)分析：兩信號(hào)分別為x(n)=nR16(n),h(n)=R8(n)，為時(shí)間離散信號(hào)，兩個(gè)信號(hào)同時(shí)在第一個(gè)圖形中表示出來(lái)，結(jié)果與信號(hào)表達(dá)式相同，結(jié)果正確。兩信號(hào)卷積x(n)*h(n)為所有h(n-k)與x(k)的乘積的累加和。所以在n=15時(shí)會(huì)出現(xiàn)卷積的最大值。X(n)和h(n)的線性卷積的結(jié)果y(n)是一個(gè)有限序列，其非零值長(zhǎng)度為N1+N2-1,最終matlab實(shí)現(xiàn)

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

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

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

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

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

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

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

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

19、.5;Rs=25;%指標(biāo)轉(zhuǎn)換wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Fs=Fs/Fs;%預(yù)畸變Wp=2*tan(wp/2); Ws=2*tan(ws/2);%由階數(shù)，零點(diǎn)向量，極點(diǎn)向量，增益實(shí)現(xiàn)模擬低通的傳遞函數(shù)N,Wn=buttord(Wp,Ws,Rp,Rs,s); z,p,k=buttap(N); Bap,Aap=zp2tf(z,p,k); %由低通原型轉(zhuǎn)向低通b,a=lp2lp(Bap,Aap,Wn); %雙線性變換實(shí)現(xiàn)模擬低通變?yōu)閿?shù)字低通bz,az=bilinear(b,a,Fs);%頻率響應(yīng)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;實(shí) 驗(yàn) 結(jié) 果 及 分 析實(shí)驗(yàn)結(jié)果：實(shí)驗(yàn)分析：上圖為，程序運(yùn)行后濾波器的幅度特性和相頻特性。從圖中可以看出：幅度響應(yīng)曲線在通帶截止頻率（f=1000hz）即在通帶內(nèi)有更大部分的幅度接近于0，可以從圖中看出，通帶截止頻率為1000Hz，阻帶截止頻率1500Hz為下降3dB是的頻率值，縱坐標(biāo)換成dB時(shí)可清楚看到。幅

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

22、易得到嚴(yán)格的線性相位，所以FIR總是最穩(wěn)定的。FIR數(shù)字濾波器主要采用非遞歸結(jié)構(gòu)，而且可以借助FFT計(jì)算。因?yàn)樗O(shè)計(jì)的是FTR數(shù)字濾波器，其單位沖激響應(yīng)h(n)必然是有限長(zhǎng)的，所以要用有限長(zhǎng)的h(n)來(lái)逼近無(wú)限長(zhǎng)的hd(n)。最有效的方法昰截?cái)鄅d(n)，即用有限長(zhǎng)的窗函數(shù)w(n)來(lái)截取hd(n)，表示h(n)= hd(n) w(n),這種方法稱為窗函數(shù)法。用窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器的步驟如下：（1） 給出所要設(shè)計(jì)的FIR數(shù)字濾波器的技術(shù)指標(biāo)，如通帶截止頻率wp、阻帶截止頻率ws、通帶衰減Rp和阻帶衰減As；（2） 根據(jù)允許的過(guò)渡帶寬度及阻帶衰減，初步選擇窗函數(shù)和N值。（3） 若選用理想

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