北理工數(shù)字信號(hào)處理實(shí)驗(yàn)報(bào)告

1、北京理工大學(xué)20112012學(xué)年第一學(xué)期數(shù)字信號(hào)處理實(shí)驗(yàn)報(bào)告數(shù)字信號(hào)處理實(shí)驗(yàn)報(bào)告姓名：仲易班級(jí)：05110901學(xué)號(hào)：20091087實(shí)驗(yàn)1 基2-FFT算法實(shí)現(xiàn)一、實(shí)驗(yàn)?zāi)康?.掌握基2-FFT的原理及具體實(shí)現(xiàn)方法。2.編程實(shí)現(xiàn)基2-FFT算法。3.深刻理解FFT算法的特點(diǎn)。二、實(shí)驗(yàn)設(shè)備與環(huán)境計(jì)算機(jī)，matlab軟件環(huán)境。三、 實(shí)驗(yàn)基礎(chǔ)理論FFT是EDF的一種快速算法，能使DFT的計(jì)算大大簡(jiǎn)化，運(yùn)算時(shí)間縮短。FFT利用了WNnk的三個(gè)固有特性。即對(duì)稱性，周期性和可約性，將長(zhǎng)序列的DFT分解為短序列的DFT，合并了DFT運(yùn)算中的某些項(xiàng)，從而減少了DFT的運(yùn)算量。FFT算法基本上可分為兩大類，即按

2、時(shí)間抽取法和按頻率抽取法。在實(shí)現(xiàn)FFT算法時(shí)，要重點(diǎn)考慮兩個(gè)問題，注意數(shù)據(jù)的讀取和存儲(chǔ)：(1)輸入輸出的排序；（2）蝶形運(yùn)算的實(shí)現(xiàn)。按時(shí)間抽取算法中輸入反序輸出順序，按頻率抽取算法中輸入順序輸出反序；運(yùn)算過程中的每一級(jí)都有N/2個(gè)蝶形運(yùn)算構(gòu)成，每一個(gè)蝶形運(yùn)算單元中，兩個(gè)節(jié)點(diǎn)變量運(yùn)算后得到的結(jié)果為下一列相同位置的節(jié)點(diǎn)變量，而和其他節(jié)點(diǎn)變量無關(guān)，可以采用原位運(yùn)算，節(jié)省存儲(chǔ)單元。另外，蝶形運(yùn)算中的復(fù)系數(shù)WNnk可以存儲(chǔ)為能及時(shí)查閱的系數(shù)表，這樣可以借閱運(yùn)算量，但是需要N/2哥復(fù)數(shù)存儲(chǔ)器。MATLAB中提供了用于計(jì)算FFT的函數(shù)fft，可將實(shí)驗(yàn)中所得到的結(jié)果與利用MATLAB中fft函數(shù)計(jì)算的結(jié)果相比

3、較，以此驗(yàn)證結(jié)果的正確性。4. 實(shí)驗(yàn)內(nèi)容及實(shí)驗(yàn)過程1. 編程實(shí)現(xiàn)序列長(zhǎng)度為N=8的按時(shí)間抽取的基2-FFT算法。給定一個(gè)8點(diǎn)序列，采用編寫的程序計(jì)算其DFT，并與MATLAB中fft函數(shù)計(jì)算的結(jié)果相比較，以驗(yàn)證結(jié)果的正確性。代碼如下：x=1 2 3 4 5 6 7 8;m=3; %X序列長(zhǎng)度為2的3次冪N=8;nxd=bin2dec(fliplr(dec2bin(1:N-1,m)+1; %求倒序列y=x(nxd); for mm=1:m %對(duì)x做m次基2分解Nz=2mm;u=1; WN=exp(-i*2*pi/Nz); for j=1:Nz/2for k=j:Nz:Nkp=k+Nz/2; t=

4、y(kp)*u; %蝶形運(yùn)算乘項(xiàng)y(kp)=y(k)-t; %加項(xiàng)y(k)=y(k)+t; endu=u*WN;endendyy1=fft(x) %對(duì)比matlab的ffty = Columns 1 through 3 36.0000 -4.0000 + 9.6569i -4.0000 + 4.0000i Columns 4 through 6 -4.0000 + 1.6569i -4.0000 -4.0000 - 1.6569i Columns 7 through 8 -4.0000 - 4.0000i -4.0000 - 9.6569iy1 = Columns 1 through 3 36

5、.0000 -4.0000 + 9.6569i -4.0000 + 4.0000i Columns 4 through 6 -4.0000 + 1.6569i -4.0000 -4.0000 - 1.6569i Columns 7 through 8 -4.0000 - 4.0000i -4.0000 - 9.6569i2、 編程實(shí)現(xiàn)序列長(zhǎng)度為N=8的按頻率抽取的基2-FFT算法。給定一個(gè)8點(diǎn)序列，采用編寫的程序計(jì)算其DFT，并與MATLAB中fft函數(shù)計(jì)算的結(jié)果相比較，以驗(yàn)證結(jié)果的正確性。代碼如下：L=2048;M=256; x=ones(1,L); n=0:M-1; h=cos(0.2*p

6、i*n); ticy1=conv(x,h);tocticX=fft(x,L+M -1);H=fft(h, L+M -1);Y=X.*H;y2=ifft(Y);toc運(yùn)行結(jié)果如下：（分別為線形卷積和快速卷積）Elapsed time is 0.003352 seconds.Elapsed time is 0.007167 seconds.L=4096;M=256; x=ones(1,L); n=0:M-1; h=cos(0.2*pi*n); ticy1=conv(x,h);tocticX=fft(x,L+M -1);H=fft(h, L+M -1);Y=X.*H;y2=ifft(Y);toc運(yùn)行

7、結(jié)果如下：（分別為線形卷積和快速卷積）Elapsed time is 0.003666 seconds.Elapsed time is 0.030134 seconds.3.將上述FFT程序推廣大序列長(zhǎng)度為N=2v的情況，要求利用原位運(yùn)算。代碼如下：L=2048;M=256; n=0:M-1; h=cos(0.2*pi*n);l=128;x=ones(1,l);q=0;ticH=fft(h,l+M-1);X=fft(x,l+M-1);Y=X.*H;y=ifft(Y);yp=y;for i=2:16q=q+l;yq=zeros(1,q),y;yp=yp,zeros(1,l);yp=yq+yp;e

8、ndtoc運(yùn)行結(jié)果如下：Elapsed time is 0.018875 seconds.4.L=2048;M=256;n=0:M-1;h=cos(0.2*pi*n);x=ones(1,L);Lx=2048;M=256;M1=M-1;N=512;L=N-M1;tich=h,zeros(1,N-M); x=zeros(1,M1),x,zeros(1,N-1); a=floor (Lx+M1-1)/(L)+1; Y=zeros(1,N); for k=0:a-1 xk=x(k*L+1:k*L+N); b=fft(xk,N); C=fft(h,N); Z=b.*C; Y(k+1,:)=ifft(Z,

9、N); endY=Y(:,M:N);Y=(Y(:)toc運(yùn)行結(jié)果如下：Elapsed time is 0.035329 seconds.實(shí)驗(yàn)2 利用DFT分析信號(hào)頻譜一、 實(shí)驗(yàn)?zāi)康?、 加深對(duì)DFT原理的理解。2、 應(yīng)用DFT分析信號(hào)的頻譜。3、 深刻理解利用DFT分析信號(hào)頻譜的原理，分析實(shí)現(xiàn)過程中出現(xiàn)的現(xiàn)象及解決方法。二、 實(shí)驗(yàn)設(shè)備與環(huán)境計(jì)算機(jī)、MATLAB軟件環(huán)境三、 實(shí)驗(yàn)基礎(chǔ)理論1、 DFT與DTFT的關(guān)系序列的N點(diǎn)DFT實(shí)際上就是序列的DTFT在N個(gè)等間隔頻率點(diǎn)上的采樣樣本。2、 利用DFT求DTFT方法一：利用差值運(yùn)算 其中為內(nèi)插函數(shù)方法二：在實(shí)際MATLAB計(jì)算中，采用增加數(shù)據(jù)的長(zhǎng)

10、度N，使得到的DFT譜線更加精細(xì)，用其包絡(luò)近似計(jì)算DTFT。如果數(shù)據(jù)量不足，可采用補(bǔ)零來增加數(shù)據(jù)長(zhǎng)度。3、 利用DFT分析連續(xù)時(shí)間信號(hào)的頻譜1） 確定時(shí)域采樣間隔T，得到離散序列2） 確定截取長(zhǎng)度M，得到M點(diǎn)離散序列，這里為窗函數(shù)3） 確定頻域采樣點(diǎn)數(shù)N，要求4） 利用FFT計(jì)算離散序列的N點(diǎn)DFT，得到5） 由計(jì)算采樣點(diǎn)的近似值。4、 可能用到的MATLAB函數(shù)與代碼實(shí)驗(yàn)中DFT運(yùn)算采用MATLAB中提供的函數(shù)fft來實(shí)現(xiàn)。DTFT可以利用MATLAB矩陣運(yùn)算的方法進(jìn)行計(jì)算。四、 實(shí)驗(yàn)內(nèi)容1、（1）n=0:3;x=2 -1 1 1;w=-pi:0.01*pi:pi;X=x*exp(-j*n*

11、w);subplot(211)plot(w,abs(X);axis tightsubplot(212)plot(w,angle(X)/pi);axis tight(2)在原基礎(chǔ)上增加代碼：n=0:3;x= 2 -1 1 1;N=4;k=0:3;Y=fft(x);w=0:0.01*pi:2*pi;X=x*exp(-j*n*w);m=0:2*pi/4:3/2*pisubplot(211);plot(w,abs(X);xlabel(Omega/pi);title(DTFT-Magnitude);axis tighthold onstem(m,abs(Y);subplot(212);plot(w,an

12、gle(X)/pi);xlabel(Omega/pi);title(DTFT-Phase);hold onstem(m,angle(Y)（3）在原本基礎(chǔ)上設(shè)置代碼subplot(211);Z=fft(x,64); % n=64DFTN=0:63;stem(abs(Z);subplot(212);stem(angle(Z)2、考察序列（1），用DFT估計(jì)的頻譜；將補(bǔ)零加長(zhǎng)到長(zhǎng)度為100點(diǎn)序列用DFT估計(jì)的頻譜。要求畫出相應(yīng)波形。程序代碼：程序輸出結(jié)果：補(bǔ)零后：程序代碼：程序輸出結(jié)果：（2）時(shí)，用DFT估計(jì)的頻譜，并畫出波形。程序代碼：程序輸出結(jié)果：（3）根據(jù)實(shí)驗(yàn)結(jié)果，分析怎樣提高頻譜分辨率。 答

13、：減小截取長(zhǎng)度或減小采樣間隔的方法可以提高頻譜分辨率。3、已知信號(hào)其中=1Hz，=2Hz，=3Hz。從的表達(dá)式可以看出，它包含三個(gè)頻率的正弦波，但是，從其時(shí)域波形來看，似乎是一個(gè)正弦信號(hào)，利用DFT做頻譜分析，確定適合的參數(shù)，使得到的頻譜的頻率分辨率符合需要。程序代碼： 程序輸出結(jié)果：利用DFT近似分析連續(xù)時(shí)間信號(hào)的頻譜。分析采用不同的采樣間隔和截取長(zhǎng)度進(jìn)行計(jì)算的結(jié)果，并最終確定適合的參數(shù)。程序代碼：程序輸出結(jié)果：通過改變步長(zhǎng)可以達(dá)到改變采樣間隔的結(jié)果，輸出不同的頻譜。實(shí)驗(yàn)5、脈沖響應(yīng)不變法設(shè)計(jì)IIR數(shù)字濾波器一、 實(shí)驗(yàn)?zāi)康?、 掌握利用脈沖響應(yīng)不變法設(shè)計(jì)IIR數(shù)字濾波器的原理及具體方法2、

14、加深理解數(shù)字濾波器和模擬濾波器之間的技術(shù)指標(biāo)轉(zhuǎn)化。3、 掌握脈沖相應(yīng)不變法設(shè)計(jì)IIR數(shù)字濾波器的優(yōu)缺點(diǎn)及適用范圍。二、 實(shí)驗(yàn)設(shè)備與環(huán)境計(jì)算機(jī)、MATLAB軟件環(huán)境三、 實(shí)驗(yàn)基礎(chǔ)理論1、 基本原理：從時(shí)域響應(yīng)出發(fā)，使數(shù)字濾波器的單位脈沖響應(yīng)h（n）模仿模擬濾波器的單位沖激響應(yīng)。2、 設(shè)計(jì)步驟：a. 確定數(shù)字濾波器性能指標(biāo)。b. 將數(shù)字濾波器頻率指標(biāo)轉(zhuǎn)換成相應(yīng)的模擬濾波器的頻率指標(biāo)c. 設(shè)計(jì)模擬濾波器實(shí)現(xiàn)方法：脈沖響應(yīng)不變法 r,p,k=residue(b,a) b,a=residue(r,p,k)d. 模擬濾波器數(shù)字化實(shí)現(xiàn)方法：bz,az=impinvar(b,a,fs)四、 實(shí)驗(yàn)內(nèi)容1、設(shè)采樣

15、頻率fs=4khz,采用脈沖響應(yīng)不變法設(shè)計(jì)一個(gè)三階巴特沃斯數(shù)字低通濾波器，其3db截止頻率為fc=1khz.首先設(shè)計(jì)巴特沃斯濾波器：N=3;Wc=0.25*pi;b,a=butter(N,Wc,s)b = 0 0 0 0.4845a = 1.0000 1.5708 1.2337 0.4845因此，然后脈沖響應(yīng)不變法計(jì)算：bz,az=impinvar(b,a)bz = 1.1702 -1.4212 0.6078 -0.0945az = 1.0000 -1.6286 0.9286 -0.1879得到：畫出模擬濾波器和數(shù)字濾波器的圖像進(jìn)行對(duì)比：H,w=freqs(b,a);subplot(221);

16、plot(w/pi,abs(H);grid on;subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H);w=0:500*pi/500;H,w=freqz(bz,az);grid on;subplot(223);plot(w/pi,abs(H);grid on;subplot(224);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;2.設(shè)采樣頻率為fs=10khz，設(shè)計(jì)數(shù)字濾波器滿足如下指標(biāo) 通帶截止頻率：1khz，通帶波動(dòng)：1db， 阻帶截止頻率：1.5khz，阻帶衰減：15db 分別采用巴特沃斯、切比雪

17、夫I型、切比雪夫ii型和橢圓模擬原型濾波器及脈沖響應(yīng)不變法進(jìn)行設(shè)計(jì)。結(jié)合實(shí)驗(yàn)結(jié)果，分別討論采用上述方法設(shè)計(jì)的數(shù)字濾波器是否都能滿足給定指標(biāo)要求，分析脈沖響應(yīng)不變法設(shè)計(jì)IIR數(shù)字濾波器的優(yōu)缺點(diǎn)及適用范圍。（1）巴特沃斯濾波器為：Wp=0.1*pi;Ws=0.15*pi;Rp=1;As=15;N=ceil(log10(10(Rp/10)-1)/(10(As/10)-1)/(2*log10(Wp/Ws)OC1=Wp/(10(Rp/10)-1)(1/(2*N)b,a=butter(N,OC1,s)N = 6OC1 = 0.3516b = 0 0 0 0 0 0 0.0019a = 1.0000 1.3

18、585 0.9227 0.3974 0.1141 0.0208 0.0019得到：bz,az=impinvar(b,a)bz = 1.0e-003 * 0.0000 0.0125 0.2573 0.5194 0.1637 0.0051 0az = 1.0000 -4.6517 9.1374 -9.6829 5.8300 -1.8889 0.2570得到bz，az后，畫出模擬濾波器和數(shù)字濾波器的圖像進(jìn)行對(duì)比：H,w=freqs(b,a);subplot(221);plot(w/pi,abs(H);grid on;subplot(222);plot(w/pi,20*log10(abs(H)/max

19、(abs(H);w=0:500*pi/500;H,w=freqz(bz,az);grid on;subplot(223);plot(w/pi,abs(H);grid on;subplot(224);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;H,w=freqs(b,a);subplot(221);plot(w/pi,abs(H);grid on;subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H);w=0:500*pi/500;H,w=freqz(bz,az);grid on;subplot(223);

20、plot(w/pi,abs(H);grid on;subplot(224);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;（2）切比雪夫型 Wp=0.2*pi; Ws=0.3*pi; Rp=1; As=15; e= e=(10(Rp/10)-1)(1/2)e =0.5088 A=10(As/20)A = 5.6234 N=ceil(acosh(A*A-1)(1/2)/e)/acosh(Ws/Wp)N = 4 b,a=cheby1(N,Rp,Wp,s)b = 0 0 0 0 0.0383a =1.0000 0.5987 0.5740 0.1842 0.

21、0430（3）切比雪夫型： Wp=0.2*pi; Ws=0.3*pi; Rp=1; As=15; N=4; b,a=cheby2(N,As,Ws,s)b =0.1778 -0.0000 1.2637 -0.0000 1.1225a =1.0000 2.3696 3.0322 1.9925 1.1225（4）橢圓模擬原型濾波器： Wp=0.2*pi; Ws=0.3*pi; Rp=1; As=15; e=0.5088; A=5.6234; k1=e/(A*A-1)(1/2)k1 =0.0919 k=Wp/Wsk = 0.6667 K,E=ellipke(k)K =2.0290E =1.2612 K

22、1,E1=ellipke(1-k1*k1)(1/2)K1 =4.1217E1 =1.0077 K2,E2=ellipke(k1)K2 = 1.6089E2 =1.5340 K3,E3=ellipke(1-k*k)(1/2)K3 =2.1483E3 = 1.2140 N=(K*K1)/(K2*K3)N = 2.4195 N2=2; b,a=ellip(N2,Rp,As,Ws,s)b = 0.1778 0 0.9379a =1.0000 0.9120 1.0523結(jié)果分析：（1） 當(dāng)給定的設(shè)計(jì)指標(biāo)相同時(shí)，選用橢圓濾波器所要求的結(jié)束N最低，切比雪夫?yàn)V波器 次之，巴特沃斯濾波器最高；（2） 從通帶的相

23、位相應(yīng)來看，橢圓濾波器雖然提供了最優(yōu)秀的幅度平方相應(yīng)，但通帶上 的相位相應(yīng)非線性較大，而巴特沃斯濾波器在通帶上具有相當(dāng)?shù)木€性相位，切比雪夫?yàn)V波器的相位特征介于兩者之間。實(shí)驗(yàn)七、窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器一、實(shí)驗(yàn)?zāi)康恼莆沾昂瘮?shù)法設(shè)計(jì)FIR數(shù)字濾波器的原理及具體方法。二、實(shí)驗(yàn)設(shè)備及環(huán)境計(jì)算機(jī)、matlab軟件及環(huán)境三、實(shí)驗(yàn)基礎(chǔ)理論1.基本原理 窗函數(shù)設(shè)計(jì)法的基本思想為，首先選擇一個(gè)適當(dāng)?shù)睦硐氲臑V波器Hd（ejw），然后用窗函數(shù)截取它的單位脈沖響應(yīng)hd（n），得到線性相位和因果的FIR濾波器。這種方法的重點(diǎn)是選擇一個(gè)合適的窗函數(shù)和理想濾波器，使設(shè)計(jì)的濾波器的單位響應(yīng)盡力逼向?yàn)V波器的單位響應(yīng)。2.設(shè)

24、計(jì)步驟（1）給定理想濾波器的頻率響應(yīng)Hd（ejw），再通帶上具有單位增益和線性相位，在阻帶上具有零響應(yīng)。（2）確定這個(gè)濾波器的單位脈沖響應(yīng)hd=sincn-n-為了得到一個(gè)h(n)長(zhǎng)度為N的因果相位FIR濾波器，我們令=N-12（3）用窗函數(shù)截取hd（n）得到所設(shè)計(jì)FIR數(shù)字濾波器h(n)h(n)= hd（n）w（n）四、實(shí)驗(yàn)內(nèi)容1.設(shè)計(jì)一個(gè)數(shù)字低通濾波器，其技術(shù)指標(biāo)如下 p=0.2, Rp=0.25db st=0.3, As=50db分別采用矩形窗、漢寧窗、海明窗、布萊克曼窗、凱瑟窗設(shè)計(jì)該濾波器。結(jié)合實(shí)驗(yàn)結(jié)果，分別討論采用上述方法設(shè)計(jì)的數(shù)字濾波器是否都能滿足給定指標(biāo)要求。2. 設(shè)計(jì)一個(gè)數(shù)字低

25、通濾波器，其技術(shù)指標(biāo)如下 下阻帶邊緣：st1=0.2, As=60db 下通帶邊緣：p1=0.35, Rp=1db 上通帶邊緣：p2=0.65, Rp=1db 上阻帶邊緣：st2=0.8, As=60db五、代碼確定長(zhǎng)度：wp=0.2*pi; wst=0.3*pi; width=wst-wp; N=ceil(1.8*pi/width)+1N = 19確定hd(n) n=0:(N-1); wc=(wp+wst)/2; alpha=(N-1)/2; hd=(wc/pi)*sinc(wc/pi)*(n-alpha) wboxcar=boxcar(N); h=hd.*wboxcarsubplot(22

26、1); stem(n,hd,filled);axis tight;xlabel(n);ylabel(hd(n); Hr,w1=zerophase(h); subplot(222);plot(w1/pi,Hr); axis tight;xlabel(omegapi);ylabel(omega); subplot(223);stem(n,h,filled);axis tight;xlabel(n);ylabel(h(n); H,w=freqz(h,1); subplot(224); plot(w/pi,20*log10(abs(H)/max(abs(H); xlabel(omegapi);ylab

27、el(dB) grid on由上圖可清楚看出：加矩形窗不滿足指標(biāo)用海明窗：wp=0.2*pi;wst=0.3*pi;width=wst-wp;N=ceil(6.6*pi/width)+1N = 67n=0:(N-1);wc=(wp+wst)/2;alpha=(N-1)/2;hd=(wc/pi)*sinc(wc/pi)*(n-alpha);whamming=hamming(N);h=hd.*whamming; subplot(221);stem(n,hd,filled);axis tight;xlabel(n);ylabel(hd(n);Hr,w1=zerophase(h);subplot(22

28、2);plot(w1/pi,Hr);axis tight;xlabel(omegapi);ylabel(omega); subplot(223);stem(n,h,filled);axis tight;xlabel(n);ylabel(h(n);H,w=freqz(h,1);subplot(224);plot(w/pi,20*log10(abs(H)/max(abs(H);xlabel(omegapi);ylabel(dB)grid on由上圖可清楚看出：加海明窗滿足指標(biāo) 漢寧窗：代碼：wp = 0.2*pi;wst = 0.3*pi; tr_width = wst-wp; N = ceil(

29、6.2*pi/tr_width)+1; n = 0:(N-1);wc = (wp+wst)/2; alpha = (N-1)/2; hd = (wc/pi)*sinc(wc/pi)*(n-alpha); w_hanning = hanning(N); h = hd.*w_hanning; subplot(221);stem(n,hd,filled); axis tight;xlabel(n);ylabel(hd(n); Hr,w1 = zerophase(h); subplot(222);plot(w1/pi,Hr); axis tight;xlabel(omega/pi);ylabel(H(omega); subplot(223);stem(n,h,filled); axis t