《數(shù)字信號處理》實驗報告

1、物理與信息工程學院數(shù)字信號處理實驗報告年級：2011級班級：信通4班 姓名：朱明貴學號：111100443老師：李娟福州大學2013 年 11 月實驗一 快速傅里葉變換(FFT)及其應(yīng)用一、實驗?zāi)康?. 在理論學習的基礎(chǔ)上，通過本實驗，加深對FFT的理解，熟悉MATLAB中的有關(guān)函數(shù)。2. 熟悉應(yīng)用FFT對典型信號進行頻譜分析的方法。3. 了解應(yīng)用FFT進行信號頻譜分析過程中可能出現(xiàn)的問題，以便在實際中正確應(yīng)用FFT。4. 熟悉應(yīng)用FFT實現(xiàn)兩個序列的線性卷積和相關(guān)的方法。二、實驗類型演示型三、實驗儀器裝有MATLAB語言的計算機四、實驗原理在各種信號序列中，有限長序列信號處理占有很重要地位，

2、對有限長序列，我們可以使用離散Fouier變換(DFT)。這一變換不但可以很好的反映序列的頻譜特性，而且易于用快速算法在計算機上實現(xiàn)，當序列x(n)的長度為N時，它的DFT定義為：反變換為： 有限長序列的DFT是其Z變換在單位圓上的等距采樣，或者說是序列Fourier變換的等距采樣，因此可以用于序列的譜分析。 FFT并不是與DFT不同的另一種變換，而是為了減少DFT運算次數(shù)的一種快速算法。它是對變換式進行一次次分解，使其成為若干小點數(shù)的組合，從而減少運算量。常用的FFT是以2為基數(shù)的，其長度 。它的效率高，程序簡單，使用非常方便，當要變換的序列長度不等于2的整數(shù)次方時，為了使用以2為基數(shù)的FF

3、T，可以用末位補零的方法，使其長度延長至2的整數(shù)次方。 (一)在運用DFT進行頻譜分析的過程中可能的產(chǎn)生三種誤差1 混疊 序列的頻譜是被采樣信號頻譜的周期延拓，當采樣速率不滿足Nyquist定理時，就會發(fā)生頻譜混疊，使得采樣后的信號序列頻譜不能真實的反映原信號的頻譜。避免混疊現(xiàn)象的唯一方法是保證采樣速率足夠高，使頻譜混疊現(xiàn)象不致出現(xiàn)，即在確定采樣頻率之前，必須對頻譜的性質(zhì)有所了解，在一般情況下，為了保證高于折疊頻率的分量不會出現(xiàn)，在采樣前，先用低通模擬濾波器對信號進行濾波。2泄漏 實際中我們往往用截短的序列來近似很長的甚至是無限長的序列，這樣可以使用較短的DFT來對信號進行頻譜分析，這種截短等

4、價于給原信號序列乘以一個矩形窗函數(shù)，也相當于在頻域?qū)⑿盘柕念l譜和矩形窗函數(shù)的頻譜卷積，所得的頻譜是原序列頻譜的擴展。泄漏不能與混疊完全分開，因為泄漏導(dǎo)致頻譜的擴展，從而造成混疊。為了減少泄漏的影響，可以選擇適當?shù)拇昂瘮?shù)使頻譜的擴散減至最小。3柵欄效應(yīng) DFT是對單位圓上Z變換的均勻采樣，所以它不可能將頻譜視為一個連續(xù)函數(shù)，就一定意義上看，用DFT來觀察頻譜就好像通過一個柵欄來觀看一個圖景一樣，只能在離散點上看到真實的頻譜，這樣就有可能發(fā)生一些頻譜的峰點或谷點被“尖樁的柵欄”所攔住，不能別我們觀察到。 減小柵欄效應(yīng)的一個方法就是借助于在原序列的末端填補一些零值，從而變動DFT的點數(shù)，這一方法實際

5、上是人為地改變了對真實頻譜采樣的點數(shù)和位置，相當于搬動了每一根“尖樁柵欄”的位置，從而使得頻譜的峰點或谷點暴露出來。 (二)用FFT計算線性卷積 用FFT可以實現(xiàn)兩個序列的圓周卷積。在一定的條件下，可以使圓周卷積等于線性卷積。一般情況，設(shè)兩個序列的長度分別為N1和N2，要使圓周卷積等于線性卷積的充要條件是FFT的長度NN1N2-1，對于長度不足N的兩個序列，分別將他們補零延長到N。 當兩個序列中有一個序列比較長的時候，我們可以采用分段卷積的方法。有兩種方法：1 重疊相加法將長序列分成與短序列相仿的片段，分別用FFT對它們作線性卷積，再將分段卷積各段重疊的部分相加構(gòu)成總的卷積輸出。2重疊保留法這

6、種方法在長序列分段時，段與段之間保留有互相重疊的部分，在構(gòu)成總的卷積輸出時只需將各段線性卷積部分直接連接起來，省掉了輸出段的直接相加。 五、實驗內(nèi)容和要求 1、一個連續(xù)信號含兩個頻率分量，經(jīng)采樣得已知，分別為1/16和1/64，觀察其頻譜；當時，不變，其結(jié)果有何不同，為什么？代碼：N=16;n=0:N-1;Df=1/16;Xn=sin(2*pi*0.125*n)+cos(2*pi*(0.125+Df)*n);Xk=fft(Xn,N);subplot(245);stem(n,abs(Xk);xlabel(n);ylabel(X(k);title(N=16,Df=1/16,頻譜圖);subplot

7、(241);stem(n,abs(Xn);xlabel(n);ylabel(X(n);title(N=16,Df=1/16,時序圖);Df=1/64;Xn=sin(2*pi*0.125*n)+cos(2*pi*(0.125+Df)*n);Xk=fft(Xn,N);subplot(246);stem(n,abs(Xk);xlabel(n);ylabel(X(k);title(N=16,Df=1/64,頻譜圖);subplot(242);stem(n,abs(Xn);xlabel(n);ylabel(X(n);title(N=16,Df=1/64,時序圖);N=128;n=0:N-1;Df=1/1

8、6;Xn=sin(2*pi*0.125*n)+cos(2*pi*(0.125+Df)*n);Xk=fft(Xn,N);subplot(247);stem(n,abs(Xk);xlabel(n);ylabel(X(k);title(N=128,Df=1/16,頻譜圖);subplot(243);stem(n,abs(Xn);xlabel(n);ylabel(X(n);title(N=128,Df=1/16,時序圖);Df=1/64;Xn=sin(2*pi*0.125*n)+cos(2*pi*(0.125+Df)*n);Xk=fft(Xn,N);subplot(248);stem(n,abs(Xk

9、);xlabel(n);ylabel(X(k);title(N=128,Df=1/64,頻譜圖);subplot(244);stem(n,abs(Xn);xlabel(n);ylabel(X(n);title(N=128,Df=1/64,時序圖);頻譜圖： 2、用FFT分別實現(xiàn)的和的線性卷積和10點、20點圓周卷積，記錄其波形，并說明他們之間的關(guān)系。代碼：function y = circonv( x1,x2,N )%UNTITLED3 Summary of this function goes here% Detailed explanation goes here if length(x1

10、)N error(N should higher than or equal to the length of x1) end if length(x2)N error() end x1=x1,zeros(1,N-length(x1); x2=x2,zeros(1,N-length(x2); m=0:1:N-1; x2=x2(mod(-m,N)+1); H=zeros(N,N); for n=1:1:N H(n,:)=cirshift(x2,n-1,N); endy=x1*H;function y = cirshift( x,m,N )%UNTITLED4 Summary of this fu

11、nction goes here% Detailed explanation goes here if length(x)N error() end x=x zeros(1,N-length(x); n=0:1:N-1; n=mod(n-m,N); y=x(n+1);n=0:11;x=0.8.n;h=1,1,1,1,1,1;N=17;x=x,zeros(1,5);h=h,zeros(1,11);X=fft(x);H=fft(h);Y=X.*H;y=ifft(Y);subplot(321);stem(x);xlabel(n);ylabel(x);title(序列x(n);subplot(323)

12、;stem(h);xlabel(n);ylabel(h);title(序列h(n);subplot(325);stem(y);xlabel(n);ylabel(y);title(序列y(n),線性卷積);n=0:11;N=20;x=0.8.n;h=1,1,1,1,1,1;y=circonv(x,h,N);x=x,zeros(1,5);h=h,zeros(1,11);subplot(322);stem(x);xlabel(n);ylabel(x);title(序列x(n);subplot(324);stem(h);xlabel(n);ylabel(h);title(序列h(n);subplot(

13、326);stem(y);xlabel(n);ylabel(y);title(序列y(n)，圓周卷積);頻譜圖：說明他們之間的關(guān)系？周期卷積是線性卷積的周期延拓。=+六、思考題 1.實驗中的信號序列和，在單位圓上的Z變換頻譜和會相同嗎？如果不同，說出哪一個低頻分量更多一些，為什么？答：不同；其中，Xd(n)的低頻分量較多，由圖形可以看出，在低頻處，Xd(n)的取值較多，呈遞減趨勢。2. 對一個有限長序列進行DFT等價于將該序列周期延拓后進行DFS展開，因為DFS也只是取其中一個周期來計算，所以FFT在一定條件下也可以用以分析周期信號序列。如果實正弦信號用16點FFT來做DFS運算，得到的頻譜是

14、信號本身的真實譜嗎？為什么？答：可以把有限長非周期序列假設(shè)為一無限長周期序列的一個主直周期，即對有限長非周期序列進行周期延拓，延拓后的序列完全可以采用DFS進行處理，即采用復(fù)指數(shù)。實驗二 IIR數(shù)字濾波器的設(shè)計一、實驗?zāi)康?. 掌握雙線性變換法及脈沖響應(yīng)不變法設(shè)計IIR數(shù)字濾波器的具體設(shè)計方法及其原理，熟悉用雙線性變換法及脈沖響應(yīng)不變法設(shè)計低通、高通和帶通IIR數(shù)字濾波器的MATLAB編程。 2. 觀察雙線性變換及脈沖響應(yīng)不變法設(shè)計的濾波器的頻域特性，了解雙線性變換法及脈沖響應(yīng)不變法的特點。 3. 熟悉Butterworth濾波器、Chebyshev濾波器和橢圓濾波器的頻率特性。 二、實驗類型

15、設(shè)計型三、實驗儀器裝有MATLAB語言的計算機四、實驗原理1 脈沖響應(yīng)不變法用數(shù)字濾波器的單位脈沖響應(yīng)序列模仿模擬濾波器的沖激響應(yīng),讓正好等于的采樣值，即，其中為采樣間隔，如果以及分別表示的拉式變換及的Z變換，則2雙線性變換法 S平面與z平面之間滿足以下映射關(guān)系：s平面的虛軸單值地映射于z平面的單位圓上，s平面的左半平面完全映射到z平面的單位圓內(nèi)。雙線性變換不存在混疊問題。 雙線性變換是一種非線性變換 ，這種非線性引起的幅頻特性畸變可通過預(yù)畸而得到校正。以低通數(shù)字濾波器為例，將設(shè)計步驟歸納如下：1. 確定數(shù)字濾波器的性能指標：通帶臨界頻率、阻帶臨界頻率、通帶波動、阻帶內(nèi)的最小衰減、采樣周期、采

16、樣頻率 ；2. 確定相應(yīng)的數(shù)字角頻率，； 3. 計算經(jīng)過預(yù)畸的相應(yīng)模擬低通原型的頻率，；4. 根據(jù)c和r計算模擬低通原型濾波器的階數(shù)N，并求得低通原型的傳遞函數(shù)； 5. 用上面的雙線性變換公式代入，求出所設(shè)計的傳遞函數(shù)； 6. 分析濾波器特性，檢查其指標是否滿足要求。 五、實驗內(nèi)容和要求1用脈沖響應(yīng)不變法設(shè)計一個巴特沃斯數(shù)字低通濾波器,要求在0-0.2內(nèi)衰耗不大于3dB,在0.6 內(nèi)衰耗不小于60dB，采樣頻率Fs=500 Hz. 代碼：Omegap=0.2*pi*500;Omegas=0.6*pi*500;rp=3;as=60;fs=500;wp=Omegap/fs;ws=Omegas/fs

17、;n,Omegac=buttord(Omegap,Omegas,rp,as,s);b,a=butter(n,Omegac,s);bz,az=impinvar(b,a,fs);w0=wp,ws;hx=freqz(bz,az,w0);h,w=freqz(bz,az);plot(w*fs/(2*pi),abs(h);gridxlabel(頻率hz)ylabel(頻率響應(yīng)幅度)頻譜圖：2 分別用脈沖響應(yīng)不變法和雙線性變換法分別設(shè)計一個巴特沃斯濾波器數(shù)字低通濾波器,使其特性逼近一個巴特沃斯模擬濾波器的性能指標如下：通帶截止頻率為22000（rad/s）,阻帶截止頻率為23000（rad/s），通帶衰耗不

18、大于3dB,阻帶衰耗不小于15dB，采樣頻率Fs=100000 Hz，觀察記錄所設(shè)計數(shù)字濾波器的幅頻特性曲線，記錄帶寬和衰減量，檢查是否滿足要求。比較這兩種方法的優(yōu)缺點。代碼（1）：omegap=2*pi*2000;omegas=2*pi*3000;rp=3;as=15;fs=10000;wp=omegap/fs;ws=omegas/fs;n,omegac=buttord(omegap,omegas,rp,as,s);b,a=butter(n,omegac,s);bz,az=impinvar(b,a,fs);w0=wp,ws;hx=freqz(bz,az,w0);h,w=freqz(bz,az

19、);plot(w*fs/(2*pi),abs(h);gridxlabel(頻率hz)ylabel(頻率響應(yīng)幅度)屏幕圖（1）：代碼（2）：omegap=2*pi*2000;omegas=2*pi*3000;rp=3;as=15;fs=10000;wp=omegap/fs;ws=omegas/fs;omegap1=2*fs*tan(wp/2);omegas1=2*fs*tan(ws/2);n,omegac=buttord(omegap1,omegas1,rp,as,s);b,a=butter(n,omegac,s);bz,az=bilinear(b,a,fs);w0=wp,ws;hx=freqz

20、(bz,az,w0);h,w=freqz(bz,az);plot(w*fs/(2*pi),abs(h);gridxlabel(頻率hz)ylabel(頻率響應(yīng)幅度)頻譜圖（2）：比較優(yōu)缺點：脈沖響應(yīng)不變法的優(yōu)點：1，模擬頻率到數(shù)字頻率的轉(zhuǎn)換時線性的。2，數(shù)字濾波器單位脈沖響應(yīng)的數(shù)字表示近似原型的模擬濾波器單位脈沖響應(yīng)，因此時域特性逼近好 缺點： 會產(chǎn)生頻譜混疊現(xiàn)象，只適合帶限濾波器雙線性變換法優(yōu)點： 克服多值映射得關(guān)系，可以消除頻率的混疊 缺點： 是非線性的，在高頻處有較大的失真。3 編寫濾波器仿真程序，完成對實際采集的心電圖信號序列（具體數(shù)據(jù)見下面）的總響應(yīng)序列，可直接調(diào)用MATLAB fi

21、lter函數(shù)實現(xiàn)仿真。附：人體心電圖采樣信號在測量過程中往往受到工業(yè)高頻干擾，所以，必須經(jīng)過低通濾波處理后，才能作為判斷心臟功能的有用信息。下面的序列就是一個實際心電圖信號采樣序列樣本，其中存在高頻干擾，實驗時，將其作為輸入信號，濾除其中的干擾成分。代碼：clear all;wp=0.2*pi;ws=0.3*pi;rp=1;rs=15;n,wn=buttord(wp/pi,ws/pi,rp,rs);b,a=butter(n,wn);N=0.5/(0.02);figure(1);freqz(b,a,N);grid;xlabel(頻率hz);ylabel(頻率響應(yīng)幅度);xn=-4,-2,0,-4

22、,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0;figure(2);subplot(2,1,1);stem(xn,.);title(心電圖信號采樣序列);yn=filter(b,a,xn);subplot(2,1,2)stem(yn,.);title(濾波后的心電圖信號);頻譜圖：六、思考題1. 雙線性變換法中和之間的關(guān)系是非線性的，在實驗中你注意到這種非

23、線性關(guān)系了嗎？從那幾種數(shù)字濾波器的幅頻特性曲線中可以觀察到這種非線性關(guān)系？ 答：在雙線性變換法中，模擬頻率與數(shù)字頻率不再是線性關(guān)系，所以一個線性相位模擬器經(jīng)過雙線性變換后得到的數(shù)字濾波器不再保持原有的線性相位了。如以上實驗過程中，采用雙線性變化法設(shè)計的butter和cheby1數(shù)字濾波器，從圖中可以看到這種非線性關(guān)系。2. 能否利用公式 完成脈沖響應(yīng)不變法的數(shù)字濾波器設(shè)計？為什么？答：IIR數(shù)字濾波器的設(shè)計實際上是求解濾波器的系數(shù)和，它是數(shù)學上的一種逼近問題，即在規(guī)定意義上（通常采用最小均方誤差準則）去逼近系統(tǒng)的特性。如果在S平面上去逼近，就得到模擬濾波器；如果在z平面上去逼近，就得到數(shù)字濾波

24、器。但是它的缺點是，存在頻率混迭效應(yīng)，故只適用于阻帶的模擬濾波器。實驗三 FIR數(shù)字濾波器的設(shè)計一、實驗?zāi)康?. 掌握用窗函數(shù)法，頻率采樣法及優(yōu)化設(shè)計法設(shè)計FIR濾波器的原理及方法，熟悉相應(yīng)的MATLAB編程。 2. 熟悉線性相位FIR濾波器的幅頻特性和相頻特性。 3. 了解各種不同窗函數(shù)對濾波器性能的影響。 二、實驗類型設(shè)計型三、實驗儀器裝有MATLAB語言的計算機四、實驗原理 線性相位實系數(shù)FIR濾波器按其N值奇偶和h(n)的奇偶對稱性分為四種： 1h(n)為偶對稱，N為奇數(shù)H(ej)的幅值關(guān)于=0，2成偶對稱。 2h(n)為偶對稱，N為偶數(shù)H(ej)的幅值關(guān)于=成奇對稱，不適合作高通。

25、3h(n)為奇對稱，N為奇數(shù)H(ej)的幅值關(guān)于=0，2成奇對稱，不適合作高通和低通。 4h(n)為奇對稱，N為偶數(shù)H(ej) =0、20，不適合作低通。FIR濾波器的常用設(shè)計方法：(一) 窗口法 窗函數(shù)法設(shè)計線性相位FIR濾波器步驟 1. 確定數(shù)字濾波器的性能要求：臨界頻率、，濾波器單位脈沖響應(yīng)長度N。 2. 根據(jù)性能要求，合理選擇單位脈沖響應(yīng)的奇偶對稱性，從而確定理想頻率響應(yīng)的幅頻特性和相頻特性。 3. 求理想單位脈沖響應(yīng)，在實際計算中，可對按M(M遠大于N)點等距離采樣，并對其求IDFT得，用代替。 4. 選擇適當?shù)拇昂瘮?shù)，根據(jù)求所需設(shè)計的FIR濾波器單位脈沖響應(yīng)。 5. 求，分析其幅頻

26、特性，若不滿足要求，可適當改變窗函數(shù)形式或長度N，重復(fù)上述設(shè)計過程，以得到滿意的結(jié)果。 窗函數(shù)的傅氏變換的主瓣決定了的過渡帶寬。的旁瓣大小和多少決定了在通帶和阻帶范圍內(nèi)波動的幅度。常用的幾種窗函數(shù)有：(1) 矩形窗 w(n)=RN(n)； (2) Hanning窗 ； (3) Hamming窗 ； (4) Blackmen窗 ； (5) Kaiser窗 。 式中Io(x)為零階貝塞爾函數(shù)。（二）頻率采樣法頻率采樣法是從頻域出發(fā)，將給定的理想頻率響應(yīng)加以等間隔采樣然后以此作為實際FIR數(shù)字濾波器的頻率特性的采樣值，即令由通過IDFT可得有限長序列將上式代入到Z變換中去可得其中()是內(nèi)插函數(shù)五、實

27、驗內(nèi)容及步驟1. N=50，編程并畫出矩形窗、漢寧窗、海明窗和布萊克曼窗的時域波形和歸一化的幅度譜，并比較各自的主要特點。代碼：wvtool(boxcar(50),hanning(50),hamming(50),blackman(50);grid;頻譜圖：并比較各自的主要特點?矩形窗優(yōu)點是主瓣比較集中，缺點是旁瓣較高，并有負旁瓣，導(dǎo)致變換中帶進了高頻干擾和泄漏，甚至出現(xiàn)負譜現(xiàn)象；漢寧窗主瓣加寬并降低，旁瓣則顯著減小，從減小泄漏觀點出發(fā)，漢寧窗優(yōu)于矩形窗但漢寧窗主瓣加寬，相當于分析帶寬加寬，頻率分辨力下降；海明窗的第一旁瓣衰減為一42dB海明窗的頻譜也是由3個矩形時窗的頻譜合成，但其旁瓣衰減速度

28、為20dB/（10oct），這比漢寧窗衰減速度慢；布萊克曼窗的幅度函數(shù)主要由五部分組成，他們的位移都不同，其幅度也是不同的WRg（w）使旁瓣再進一步抵消。旁瓣峰值幅度進一步增加，其幅度譜主瓣寬度是矩形窗的3倍。設(shè)計程序時用backman函數(shù)調(diào)用。2. N=15，帶通濾波器的兩個通帶邊界分別是，。用漢寧（Hanning）窗設(shè)計此線性相位帶通濾波器，觀察它的實際3dB和20dB帶寬。N=45，重復(fù)這一設(shè)計，觀察幅頻和相位特性的變化，注意長度N變化的影響。N=15漢寧代碼：clc;clear alln=15;w1=0.3;w2=0.5;wn=w1,w2;b2=fir1(n,wn,hanning(n+1);freqz(b2,1);title(漢寧窗，N=15);N=15漢寧頻譜：N=45漢寧代碼：clc;clear alln=45;w1=0.3;w2=0.5;wn=w1,w2;b2=fir1(n,wn,hanning(n+1);freqz(b2,1);title(漢寧窗，N=45);N=45漢寧頻譜：3. 分別改用矩形窗和Blackman窗，設(shè)計2中的帶通濾波器，