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