數(shù)字信號處理實驗報告-7

數(shù)字信號處理實驗報告專業(yè)班級電信1101姓名曾文學(xué)號20111186020指導(dǎo)老師吳莉華中科技大學(xué)武昌分校2014年5月21日實驗一信號、系統(tǒng)及系統(tǒng)響應(yīng)實驗?zāi)康?1)加深對離散線性移不變(LSI)系統(tǒng)基本理論的理解，明確差分方程與系統(tǒng)函數(shù)之間的關(guān)系。(2)初步了解用MATLAB語言進行離散時間系統(tǒng)研究的基本方法。(3)掌握求解離散時間系統(tǒng)沖激響應(yīng)和階躍響應(yīng)程序的編寫方法，了解常用子函數(shù)。(4)通過實驗進一步理解卷積定理，了解卷積的過程。(5)了解MATLAB中有關(guān)卷積的子函數(shù)及其應(yīng)用方法。實驗原理(1)離散LSI系統(tǒng)的響應(yīng)與激勵由離散時間系統(tǒng)的時域和頻域分析方法可知，一個線性移不變離散系統(tǒng)可以用線性常系數(shù)差分方程表示：系統(tǒng)函數(shù)H(z)反映了系統(tǒng)響應(yīng)與激勵的關(guān)系。一旦上式中的bm和ak的數(shù)據(jù)確定了，則系統(tǒng)的性質(zhì)也就確定了。其中特別注意：a0必須進行歸一化處理，即a0＝1。對于復(fù)雜信號激勵下的線性系統(tǒng)，可以將激勵信號在時域中分解為單位脈沖序列或單位階躍序列，把這些單元激勵信號分別加于系統(tǒng)求其響應(yīng)，然后把這些響應(yīng)疊加，即可得到復(fù)雜信號加于系統(tǒng)的零狀態(tài)響應(yīng)。因此，求解系統(tǒng)的沖激響應(yīng)和階躍響應(yīng)尤為重要。由圖1－1可以看出一個離散LSI系統(tǒng)響應(yīng)與激勵的關(guān)系。同時，圖1－1顯示了系統(tǒng)時域分析方法和z變換域分析法的關(guān)系。如果已知系統(tǒng)的沖激響應(yīng)h(n)，則對它進行z變換即可求得系統(tǒng)函數(shù)H(z)；反之，知道了系統(tǒng)函數(shù)H(z)，對其進行z逆變換，即可求得系統(tǒng)的沖激響應(yīng)h(n)。y(n)=x(n)*h(n)x(n)h(n)y(n)=x(n)*h(n)x(n)h(n)H(z)H(z)Y(z)=X(z)H(z)Y(z)=X(z)H(z)X(z)X(z)圖1－1離散LSI系統(tǒng)響應(yīng)與激勵的關(guān)系圖1－1離散LSI系統(tǒng)響應(yīng)與激勵的關(guān)系(2)離散LSI系統(tǒng)的線性卷積由理論學(xué)習(xí)我們已知，對于線性移不變離散系統(tǒng)，任意的輸入信號x(n)可以用及其位移的線性組合來表示，即當(dāng)輸入為 時，系統(tǒng)的輸出y(n)=h(n)，由系統(tǒng)的線性移不變性質(zhì)可以得到系統(tǒng)對x(n)的響應(yīng)y(n)為稱為離散系統(tǒng)的線性卷積，簡記為也就是說，如果已知系統(tǒng)的沖激響應(yīng)，將輸入信號與系統(tǒng)的沖激響應(yīng)進行卷積運算，即可求得系統(tǒng)的響應(yīng)。實驗內(nèi)容(1)已知一個因果系統(tǒng)的差分方程為滿足初始條件y(-1)=0,x(-1)=0,求系統(tǒng)的沖激響應(yīng)和階躍響應(yīng)。編寫仿真程序，并調(diào)試得到結(jié)果，進行分析。(2)已知兩個信號序列：f1=0.8n(0<n<20)f2=u(n)(0<n<10)求兩個序列的卷積和。編寫仿真程序，并調(diào)試得到結(jié)果，進行分析。程序：①a=[1,0,1/3,0];%分母多項式系數(shù)b=[1/6,1/2,1/2,1/6];%分子多項式系數(shù)N=32;n=0:N-1;%一維數(shù)組，含N個分量，間隔為1hn=impz(b,a,n);%單位沖激響應(yīng)gn=dstep(b,a,n);%單位階躍響應(yīng)subplot(1,2,1),stem(n,hn,'k');title('系統(tǒng)的單位沖激響應(yīng)');ylabel('h(n)');xlabel('n');subplot(1,2,2),stem(n,gn,'k');title('系統(tǒng)的單位階躍響應(yīng)');ylabel('g(n)');xlabel('n');②nf1=0:20;f1=0.8.^nf1;subplot(2,2,1);stem(nf1,f1,'filled');title('f1(n)');nf2=0:10;lf2=length(nf2);f2=ones(1,lf2);subplot(2,2,2);stem(nf2,f2,'filled');title('f2(n)');y=conv(f1,f2);subplot(2,1,2);stem(y,'filled');title('y(n)');。實驗二頻域采樣1.實驗?zāi)康模?1)掌握頻率域采樣會引起時域周期化的概念，以及頻率域采樣定理及其對頻域采樣點數(shù)選擇的指導(dǎo)作用。(2)會用MATLAB語言進行頻域抽樣與恢復(fù)時程序的編寫方法。實驗原理：了解頻域采樣定理的要點，掌握采樣理論的結(jié)論：“頻域采樣時域信號周期延拓”。實驗內(nèi)容：（1）頻域采樣理論的驗證。給定信號如下：（2）編寫程序分別對頻譜函數(shù)在區(qū)間上等間隔采樣32和16點，得到，再分別對進行32點和16點IFFT，得到。（3）分別畫出、的幅度譜，并繪圖顯示x(n)、的波形，進行對比和分析，驗證總結(jié)頻域采樣理論。程序：①M=27;N=32;n=0:M-1;%產(chǎn)生M長三角波序列x(n)xa=1:ceil(M/2);%生成序列n+1xb=floor(M/2):-1:1;%生成序列27-n；floor是向下取整，ceil是向上取整xn=[xa,xb];Xk=fft(xn,1024);%1024點FFT[x(n)],用于近似序列x(n)的FTX32k=fft(xn,32);%32點FFT[x(n)]x32n=ifft(X32k);%32點IFFT[X32(k)]得到x32(n)X16k=X32k(1:2:N);%隔點抽取X32k得到X16(K)x16n=ifft(X16k,N/2);%16點IFFT[X16(k)]得到x16(n)subplot(3,2,2);stem(n,xn,'.');%stem畫離散序列圖boxon%給圖形加邊框title('(b)三角波序列x(n)');xlabel('n');ylabel('x(n)');axis([0,32,0,20])k=0:1023;wk=2*k/1024;%產(chǎn)生1024點DFT對應(yīng)的采樣點頻率(關(guān)于π歸一化值)subplot(3,2,1);%畫子圖，3表示行數(shù)，2表示列數(shù)，1表示當(dāng)前子圖的序號數(shù)(以行元素優(yōu)先順序排列)。plot(wk,abs(Xk));%繪制1024點DFT的幅頻特性圖title('(a)FT[x(n)]');xlabel('\omega/\pi');%'\'是轉(zhuǎn)義符號等價于w/piylabel('|X(e^j^\omega)|');%等價于e的jw次方axis([0,1,0,200])%確定軸的范圍，橫軸從0~1，縱軸從0~到200k=0:N/2-1;subplot(3,2,3);stem(k,abs(X16k),'.');boxontitle('(c)16點頻域采樣');xlabel('k');ylabel('|X_1_6(k)|');%下劃線表示下標(biāo)(轉(zhuǎn)義字符‘_’)axis([0,8,0,200])n1=0:N/2-1;subplot(3,2,4);stem(n1,x16n,'.');boxontitle('(d)16點IDFT[X_1_6(k)]');xlabel('n');ylabel('x_1_6(n)');axis([0,32,0,20])k=0:N-1;subplot(3,2,5);stem(k,abs(X32k),'.');boxontitle('(e)32點頻域采樣');xlabel('k');ylabel('|X_3_2(k)|');axis([0,16,0,200])n1=0:N-1;subplot(3,2,6);stem(n1,x32n,'.');boxontitle('(f)32點IDFT[X_3_2(k)]');xlabel('n');ylabel('x_3_2(n)');axis([0,32,0,20])②wp=0.25*pi;%濾波器的通帶截止頻率ws=0.4*pi;%濾波器的阻帶截止頻率Rp=1;As=15;%輸入濾波器的通阻帶衰減指標(biāo)ripple=10^(-Rp/20);Attn=10^(-As/20);Fs=100;T=1/Fs;Omgp=(2/T)*tan(wp/2);%原型通帶頻率預(yù)修正Omgs=(2/T)*tan(ws/2);%原型阻帶頻率預(yù)修正[n,Omgc]=buttord(Omgp,Omgs,Rp,As,'s');%計算階數(shù)n和3dB截止頻率[z0,p0,k0]=buttap(n);%歸一化原型設(shè)計ba=k0*real(poly(z0));%求原型濾波器系數(shù)baa=real(poly(p0));%求原型濾波器系數(shù)a[ba1,aa1]=lp2lp(ba,aa,Omgc);%變換為模擬低通濾波器系數(shù)b,a[bd,ad]=bilinear(ba1,aa1,Fs)%用雙線性變換法求數(shù)字濾波器系數(shù)b,a[sos,g]=tf2sos(bd,ad);%由直接型變換為級聯(lián)型[H,w]=freqz(bd,ad);dbH=20*log10((abs(H)+eps)/max(abs(H)));%化為分貝值subplot(2,2,1),plot(w/pi,abs(H));ylabel('|H|');title('幅度響應(yīng)');axis([0,1,0,1.1]);set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);gridsubplot(2,2,2),plot(w/pi,angle(H)/pi);ylabel('\phi');title('相位響應(yīng)');axis([0,1,-1,1]);set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);set(gca,'YTickMode','manual','YTick',[-1,0,1]);gridsubplot(2,2,3),plot(w/pi,dbH);title('幅度響應(yīng)(dB)');ylabel('dB');xlabel('頻率(\pi)');axis([0,1,-40,5]);set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);set(gca,'YTickMode','manual','YTick',[-50,-15,-1,0]);gridsubplot(2,2,4),zplane(bd,ad);axis([-1.1,1.1,-1.1,1.1]);title('零極圖');思考題：如果序列x(n)的長度為M，希望得到其頻譜在上的N點等間隔采樣，當(dāng)N<M時，如何用一次最少點數(shù)的DFT得到該頻譜采樣？答：先對原序列x(n)以N為周期進行周期延拓后取主值區(qū)序列，再計算N點DFT則得到N點頻域采樣：實驗三信號的快速傅里葉變換實驗1、實驗?zāi)康模菏箤W(xué)生進一步了解快速付里葉變換的理論；及用Matlab軟件設(shè)計快速付里葉變換的方法，使學(xué)生進一步了解數(shù)字信號的頻譜概念。2、實驗內(nèi)容與步驟：1）根據(jù)數(shù)字信號序列的長度決定2的整數(shù)次冪N，用Matlab軟件設(shè)計出快速付里葉變換及反變換的程序，或按照范例程序進行修改，輸入信號序列，運行程序，檢查程序是否有錯。2）添加繪圖語句，畫出數(shù)字信號變換前后的波形。運行程序，觀察變換前后的波形有什么不同。3）添加繪圖語句，畫出變換后的頻譜圖，運行程序，仔細觀察頻譜圖，了解變換的意義。3.實驗程序及現(xiàn)象m=10;fori=1:m+1//系統(tǒng)函數(shù)中B為分子多項式的系數(shù)，用矩陣的形式表示ifi==1B(i)=1;elseifi==m+1B(i)=-1;elseB(i)=0;endendA=[1,-1];//系統(tǒng)函數(shù)中A為分母多項式的系數(shù)N=8192;[H,f]=freqz(B,A,N);plot(f*25/pi,abs(H));grid;figure,plot(f*25/pi,angle(H));grid;figure,zplane(B,A);k=0:N-1;f=2*k/N;load('RawData.mat');x=rawdata(1,1:N);w=filter(B,A,x);X=abs(fft(x,N));W=abs(fft(w,N));figure;plot(x);title('輸入信號');figure;plot(w);title('輸出信號');figure;plot(f,abs(X));title('輸入信號的幅頻響應(yīng)');figure;plot(f,abs(W));title('輸出信號的幅頻響應(yīng)');figure;plot(f,abs(X),'b',f,abs(W),'r');實驗現(xiàn)象如下圖圖形如下：（3）極點分布圖，一共有9個極點（4）輸入信號與（5）輸出信號的不同，是由于輸入信號經(jīng)過了低通濾波器的濾波，濾除了不必要的波形。（6）輸入頻響與（7）輸出頻響的區(qū)別，也與低通濾波器有關(guān)快速傅里葉變換的目的：快速傅里葉變化是為了為應(yīng)用各種信號的實時處理提供條件的方法，大大的提高了傅里葉變化的計算量。實驗四用雙線性變換法設(shè)計IIR數(shù)字濾波器1.實驗?zāi)康?1)熟悉用雙線性變換法設(shè)計IIR數(shù)字濾波器的原理與方法。(2)掌握用雙線性變換法設(shè)計數(shù)字濾波器的計算機仿真方法。(3)了解MATLAB有關(guān)雙線性變換法的子函數(shù)。2.實驗涉及的MATLAB子函數(shù)bilinear功能：雙線性變換——將s域映射到z域的標(biāo)準方法，將模擬濾波器變換成離散等效濾波器。調(diào)用格式：[numd,dend]=bilinear(num,den,Fs)；將模擬域傳遞函數(shù)變換為數(shù)字域傳遞函數(shù)，F(xiàn)s為取樣頻率。3.實驗原理(1)雙線性變換法是將整個s平面映射到整個z平面，其映射關(guān)系為雙線性變換法克服了脈沖響應(yīng)不變法從s平面到z平面的多值映射的缺點，消除了頻譜混疊現(xiàn)象。但其在變換過程中產(chǎn)生了非線性的畸變，在設(shè)計IIR數(shù)字濾波器的過程中需要進行一定的預(yù)修正。(2)雙線性變換法設(shè)計IIR數(shù)字濾波器的步驟：①輸入給定的數(shù)字濾波器設(shè)計指標(biāo)；②根據(jù)公式 進行預(yù)修正，將數(shù)字濾波器指標(biāo)轉(zhuǎn)換成模擬濾波器設(shè)計指標(biāo)；③確定模擬濾波器的最小階數(shù)和截止頻率；④計算模擬低通原型濾波器的系統(tǒng)傳遞函數(shù)；⑤利用模擬域頻率變換法，求解實際模擬濾波器的系統(tǒng)傳遞函數(shù)；⑥用雙線性變換法將模擬濾波器轉(zhuǎn)換為數(shù)字濾波器。 4.實驗內(nèi)容采用雙線性變換法設(shè)計一個巴特沃斯數(shù)字低通濾波器，要求：濾波器采樣頻率思考題用雙線性變換法設(shè)計數(shù)字濾波器過程中，變換公式中T的取值，對設(shè)計結(jié)果有無影響?為什么?答：沒有，一般取2/T=1，方便計算實驗程序現(xiàn)象程序如下：wp=0.25*pi;%濾波器的通帶截止頻率ws=0.4*pi;%濾波器的阻帶截止頻率Rp=1;As=15;%輸入濾波器的通阻帶衰減指標(biāo)ripple=10^(-Rp/20);Attn=10^(-As/20);Fs=100;T=1/Fs;Omgp=(2/T)*tan(wp/2);%原型通帶頻率預(yù)修正Omgs=(2/T)*tan(ws/2);%原型阻帶頻率預(yù)修正[n,Omgc]=butter(Omgp,Omgs,Rp,As,'s');%計算階數(shù)n和3dB截止頻率[z0,p0,k0]=buttap(n);%歸一化原型設(shè)計ba=k0*real(poly(z0));%求原型濾波器系數(shù)baa=real(poly(p0));%求原型濾波器系數(shù)a[ba1,aa1]=lp2lp(ba,aa,Omgc);%變換為模擬低通濾波器系數(shù)b,a[bd,ad]=bilinear(ba1,aa1,Fs)%用雙線性變換法求數(shù)字濾波器系數(shù)b,a[sos,g]=tf2sos(bd,ad);%由直接型變換為級聯(lián)型[H,w]=freqz(bd,ad);dbH=20*log10((abs(H)+eps)/max(abs(H)));%化為分貝值subplot(2,2,1),plot(w/pi,abs(H));ylabel('|H|');title('幅度響應(yīng)');axis([0,1,0,1.1]);set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);gridsubplot(2,2,2),plot(w/pi,angle(H)/pi);ylabel('\phi');title('相位響應(yīng)');axis([0,1,-1,1]);set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);set(gca,'YTickMode','manual','YTick',[-1,0,1]);gridsubplot(2,2,3),plot(w/pi,dbH);title('幅度響應(yīng)(dB)');ylabel('dB');xlabel('頻率(\pi)');axis([0,1,-40,5]);set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);set(gca,'YTickMode','manual','YTick',[-50,-15,-1,0]);gridsubplot(2,2,4),zplane(bd,ad);axis([-1.1,1.1,-1.1,