實(shí)驗(yàn)一、數(shù)字信號(hào)處理基礎(chǔ)

實(shí)驗(yàn)一：數(shù)字信號(hào)處理基礎(chǔ)一、實(shí)驗(yàn)?zāi)康?、熟悉并掌握離散系統(tǒng)的差分方程表示法；2、加深對(duì)沖激響應(yīng)和卷積分析方法的理解；3、加深對(duì)離散信號(hào)的DFT的理解；4、熟悉離散系統(tǒng)的頻率響應(yīng)分析方法、加深對(duì)零、極點(diǎn)分布的概念理解；5、掌握Matlab常用函數(shù)的使用方法。二、實(shí)驗(yàn)原理1.LTI系統(tǒng)的表示在時(shí)域中，離散時(shí)間系統(tǒng)對(duì)輸入信號(hào)或者延遲信號(hào)進(jìn)行運(yùn)算處理，生成具有所需特性的輸出信號(hào)，具體框圖如下:其輸入、輸出關(guān)系可用以下差分方程描述：輸入信號(hào)分解為沖激信號(hào),記系統(tǒng)單位沖激響應(yīng)，則系統(tǒng)響應(yīng)為如下的卷積計(jì)算式：當(dāng)時(shí)，h[n]是有限長(zhǎng)度的(n=1:M)，稱系統(tǒng)為FIR系統(tǒng)；反之，稱系統(tǒng)為IIR系統(tǒng)。2.離散傅里葉變換(DFT)的定義N點(diǎn)序列的DFT和IDFT變換定義式如下：, ,利用旋轉(zhuǎn)因子具有周期性，可以得到快速算法（FFT）。3.離散時(shí)間系統(tǒng)的變換域分析方法離散系統(tǒng)的時(shí)域方程為其變換域分析方法如下：頻域:系統(tǒng)的頻率響應(yīng)為:Z域:系統(tǒng)的轉(zhuǎn)移函數(shù)為:分解因式:，其中和稱為零、極點(diǎn)。三、預(yù)習(xí)要求1、在MATLAB中，熟悉利用函數(shù)y=filter(b,a,x)實(shí)現(xiàn)差分方程的使用方法；2、在MATLAB中，熟悉用函數(shù)y=conv(x,h)計(jì)算卷積的使用方法；3、在MATLAB中，熟悉用函數(shù)y=impz(b,a)求系統(tǒng)沖激響應(yīng)的使用方法。4、在MATLAB中，熟悉用函數(shù)y=stepz(b,a)求系統(tǒng)沖激響應(yīng)的使用方法。5、在MATLAB中，函數(shù)和計(jì)算N點(diǎn)序列的DFT正、反變換的使用方法。6、在MATLAB中，熟悉函數(shù)tf2zp、zplane、freqz、residuez、zp2sos的使用。函數(shù)的使用方法：1、filter:這是一個(gè)一維數(shù)字濾波器；也可以用作卷積運(yùn)算，結(jié)果與conv函數(shù)一樣。y=filter(A,B,X)，x為待處理序列，H=B/A為一個(gè)系統(tǒng)函數(shù)，y的結(jié)果是H與X相乘后作傅氏反變換的結(jié)果，如果A=1，則y為兩個(gè)序列的卷積。2、conv:求卷積函數(shù)，y=conv(x,h)，y是x與h卷積的結(jié)果。3、impz:求系統(tǒng)的沖激響應(yīng)。[h,t]=impz(b,a):b、a分別為系統(tǒng)的傳遞函數(shù)的分子和分母系數(shù)向量。返回系統(tǒng)(b,a)的沖激響應(yīng)h和相應(yīng)的時(shí)間軸向量。4、stepz:求系統(tǒng)的單位階躍響應(yīng)。[s,n]=stepz(b,a):求解離散系統(tǒng)的單位階躍響應(yīng)，其中b、a為向量，n為時(shí)間序列[s,n]=stepz(b,a,N)：求解離散系統(tǒng)的單位階躍響應(yīng)，N為采樣點(diǎn)數(shù)stepz(b,a)：在當(dāng)前窗口，用stem(n,s)繪出圖形。5、求離散時(shí)間序列的傅立葉變換和反變換：X=fft(x,N)，x=ifft(X,N)其中，X為x的DFT變換，而x則是X的傅立葉反變換，N為計(jì)算的點(diǎn)數(shù)，可以能過(guò)上面兩個(gè)函數(shù)實(shí)現(xiàn)DFT。6、tf2ss：由傳遞函數(shù)到狀態(tài)空間的轉(zhuǎn)換。[A,B,C,D]=tf2ss(NUM,DEN),NUM,DEM為向量H(s)=NUM(s)/DEN(s)到x=Ax+Bu,y=Cx+Du7、zplane(b,a)：繪制由行向量b和a構(gòu)成的系統(tǒng)函數(shù)的零極點(diǎn)分布圖；zplane(z,p)：繪制由列向量z確定的零點(diǎn)、列向量p確定的極點(diǎn)構(gòu)成的零極點(diǎn)分布圖。8、freqz：求離散時(shí)間系統(tǒng)的頻率響應(yīng)函數(shù)；[h,w]=freqz(b,a,n)。可以得到數(shù)字濾波器的n點(diǎn)復(fù)頻響應(yīng)值，這n個(gè)點(diǎn)均勻地分布這[0,pi]上，并將這n個(gè)頻點(diǎn)的頻率記錄在w中，相應(yīng)的頻響值記錄在h中，n缺省時(shí)取512點(diǎn)。[h,f]=freqz(b,a,n,Fs)，用于對(duì)H(jw)在[0,Fs/2]上等間隔采樣n點(diǎn)，采樣點(diǎn)頻率及相應(yīng)頻響值分別記錄在f和h中，由用戶指定Fs(以Hz為單位)的值。9、[r,p,c]=residuez(b,a)：把b(z)/a(z)展開成部分分式；

[b,a]=residuez(r,p,c):根據(jù)部分分式的r、p、c數(shù)組，返回有理多項(xiàng)式。10、zp2sos:變系統(tǒng)零極點(diǎn)形式為二階分割形式[SOS,G]=ZP2SOS(Z,P,K)；SOS為矩陣形式，G為增益?！?.四、實(shí)驗(yàn)內(nèi)容1.Matlab函數(shù)conv和filter的使用1.1實(shí)驗(yàn)要求以下程序中分別使用conv和filter函數(shù)計(jì)算h和x的卷積y和y1，運(yùn)行程序，并分析y和y1是否有差別，為什么要使用x[n]補(bǔ)零后的x1來(lái)產(chǎn)生y1；具體分析當(dāng)h[n]有M個(gè)值，x[n]有M個(gè)值，使用filter完成卷積功能，需要如何補(bǔ)零？%ProgramP2_7clf;h=[321-210-403];%impulseresponsex=[1-23-4321];%inputsequencey=conv(h,x);n=0:14;subplot(2,1,1);stem(n,y);xlabel('Timeindexn');ylabel('Amplitude');title('OutputObtainedbyConvolution');grid;x1=[xzeros(1,8)];y1=filter(h,1,x1);subplot(2,1,2);stem(n,y1);xlabel('Timeindexn');ylabel('Amplitude');title('OutputGeneratedbyFiltering');grid;1.2實(shí)驗(yàn)結(jié)果與分析上面程序的運(yùn)行結(jié)果如下：將程序作如下改動(dòng)：h=[321-210-403];%impulseresponsex=[1-23-4321];%inputsequencey=conv(h,x)y1=filter(h,1,x)y2=filter(x,1,h)結(jié)果輸出如下：y=Columns1through103-46-1090172-910Columns11through15-2-20563y1=3-46-109017y2=3-46-1090172-9從結(jié)果分析可知，conv(h,x)是專門用來(lái)做卷積運(yùn)算的，得出的結(jié)果長(zhǎng)度為(nh+nx-1)，而對(duì)于filter(h,1,x)，得出的結(jié)果長(zhǎng)度與第三個(gè)參數(shù)x相同，而對(duì)于前面的運(yùn)算結(jié)果與conv的結(jié)果一致，所以，用filter函數(shù)進(jìn)行卷積運(yùn)算時(shí)要對(duì)x進(jìn)行補(bǔ)0操作。當(dāng)h[n]有M個(gè)值，x[n]有M個(gè)值時(shí)，需要令x=[xzeros(M-1,1)]。………………..2.Matlab函數(shù)impz和stepz的使用2.1實(shí)驗(yàn)要求編制程序用impz和stepz分別求解下列兩個(gè)系統(tǒng)的單位沖激響應(yīng)和階躍響應(yīng)：y[n]+0.75y[n-1]+0.125y[n-2]=x[n]–x[n-1]y[n]=0.25*(x[n-1]+x[n-2]+x[n-3]+x[n-4])給出理論計(jì)算結(jié)果和程序計(jì)算結(jié)果并討論。2.2程序代碼N=16;n=0:N-1;y1=[1,.75,0.125];x1=[1,-1];y2=1;x2=[0.25,0.25,0.25,0.25];h1=6*(-0.5).^n.*(n>=0)-5*(-0.25).^n.*(n>=0);%式1的單位沖激響應(yīng)理論計(jì)算值h2=2*(-0.5).^n.*(n>=0)-(-0.25).^n.*(n>=0);%式1的單位階躍響應(yīng)理論計(jì)算值h3=[0,0.25,0.25,0.25,0.25,zeros(1,11)];h4=0.25*((n>=1)+(n>=2)+(n>=3)+(n>=4));subplot(4,2,1);impz(x1,y1,N);%求系統(tǒng)的沖激響應(yīng)title('式1的沖激響應(yīng)matlab結(jié)果');subplot(4,2,2);stem(n,h1);title('式1的沖激響應(yīng)理論計(jì)算結(jié)果');xlabel('n(samples)');ylabel('Amplitude');subplot(4,2,3);stepz(x1,y1,N);%求系統(tǒng)的階躍響應(yīng)title('式1的階躍響應(yīng)matlab結(jié)果');subplot(4,2,4);stem(n,h2);title('式1的階躍響應(yīng)理論計(jì)算結(jié)果');xlabel('n(samples)');ylabel('Amplitude');subplot(4,2,5);impz(x2,y2,N);title('式2的沖激響應(yīng)matlab結(jié)果');subplot(4,2,6);stem(n,h3);title('式2的沖激響應(yīng)理論計(jì)算結(jié)果');xlabel('n(samples)');ylabel('Amplitude');subplot(4,2,7);stepz(x2,y2,N);title('式2的階躍響應(yīng)matlab結(jié)果');subplot(4,2,8);stem(n,h4);title('式2的階躍響應(yīng)理論計(jì)算結(jié)果');xlabel('n(samples)');ylabel('Amplitude');………………..2.3實(shí)驗(yàn)結(jié)果與分析理論計(jì)算：對(duì)式1：y[n]+0.75y[n-1]+0.125y[n-2]=x[n]–x[n-1]單位沖擊響應(yīng)為：階躍響應(yīng)為：對(duì)式2：y[n]=0.25*(x[n-1]+x[n-2]+x[n-3]+x[n-4])單位沖擊響應(yīng)為：階躍響應(yīng)為：由結(jié)果圖像可得，通過(guò)impz和stepz計(jì)算的單位沖激響應(yīng)和階躍響應(yīng)與理論計(jì)算結(jié)果基本一致?！?.3.用N點(diǎn)DFT計(jì)算2N點(diǎn)實(shí)數(shù)的DFT3.1實(shí)驗(yàn)要求2N點(diǎn)實(shí)數(shù)序列N=64。用一個(gè)64點(diǎn)的復(fù)數(shù)FFT程序，一次算出，并繪出的圖形。3.2實(shí)現(xiàn)原理v[n]是長(zhǎng)度為2N的實(shí)序列，V[k]表示該實(shí)序列的2N點(diǎn)DFT。定義兩個(gè)長(zhǎng)度均為N的實(shí)序列g(shù)[n]和h[n]為g[n]=v[2n]，h[n]=v[2n+1]，。G[k]和H[k]表示它們的N點(diǎn)DFT。通過(guò)下式可以算出v[n]的2N點(diǎn)傅立葉變換。在matlab中，可以運(yùn)用函數(shù)X=fft(x,N)直接計(jì)算出x[n]的N點(diǎn)傅立葉變換?！?.3.3實(shí)現(xiàn)代碼N=64;n=0:N-1;k=0:2*N-1;x=cos(14*pi*k/N)+0.5*cos(38*pi*k/N);g=cos(14*pi*2*n/N)+0.5*cos(38*pi*2*n/N);%g[n]=v[2n]h=cos(14*pi*(2*n+1)/N)+0.5*cos(38*pi*(2*n+1)/N);%h[n]=v[2n+1]G=fft(g,N);H=fft(h,N);form=1:N%利用G,K求X，2N點(diǎn)的傅氏變換X(m)=G(m)+exp(-j*2*pi*(m-1)/(2*N))*H(m);endform=(N+1):2*NX(m)=G(m-N)+exp(-j*2*pi*(m-1)/(2*N))*H(m-N);endX1=fft(x,2*N);%直接計(jì)算2N點(diǎn)傅區(qū)變換dX=X1-X;%作差看兩種算法的區(qū)別subplot(3,1,1);stem(k,abs(X));title('利用N點(diǎn)DFT求2N點(diǎn)DFT');xlabel('k');ylabel('|X|');subplot(3,1,2);stem(k,abs(X1));title('直接求2N點(diǎn)DFT');xlabel('k');ylabel('|X1|');subplot(3,1,3);stem(k,abs(dX));title('兩種算法的差值');xlabel('k');ylabel('|cX|');………………..3.4實(shí)驗(yàn)結(jié)果與分析從結(jié)果圖像可以看出，運(yùn)用N點(diǎn)DFT求2N點(diǎn)DFT的方法求得的結(jié)果與直接求2N點(diǎn)DFT的結(jié)果基乎一樣，雖然其中存在一定的差別，但是從工程上講，這種差別可以忽略不計(jì)?！?.4.用DFT驗(yàn)證Z變換4.1實(shí)驗(yàn)要求已知某序列x(n)的Z變換在單位圓上的N=64等分樣點(diǎn)為：。用N點(diǎn)IFFT程序計(jì)算出和。4.2實(shí)現(xiàn)原理在matlab中，可以運(yùn)用函數(shù)x=ifft(X,N)來(lái)計(jì)算頻域函數(shù)X的N點(diǎn)傅立葉反變換IDFT。而對(duì)于理論計(jì)算，可以運(yùn)用傅立葉變換對(duì)進(jìn)行變換：由于DFT相當(dāng)于無(wú)限長(zhǎng)序列的傅立葉變換的一個(gè)抽樣，所以對(duì)于理論值的DFT只要取DTFT的有限個(gè)點(diǎn)即可。………………..4.3實(shí)現(xiàn)代碼N=64;k=0:N-1;X=1./(1-0.8*exp(-j*2*pi*k/N));%原系統(tǒng)函數(shù)x1=ifft(X,N);%運(yùn)用matlab求解傅立葉反變換x2=0.8.^k;%理論計(jì)算結(jié)果dx=x2-x1;%理論結(jié)果與數(shù)值計(jì)算結(jié)果作差對(duì)比subplot(3,1,1);stem(k,x1);axis([07001.01]);%對(duì)圖像的x軸修正title('matlab計(jì)算結(jié)果');xlabel('n');ylabel('x1');subplot(3,1,2);stem(k,x2);title('理論計(jì)算結(jié)果');xlabel('n');ylabel('x2');subplot(3,1,3);stem(k,dx);title('matlab計(jì)算結(jié)果與理論結(jié)果的差值');xlabel('n');ylabel('dx');………………..4.4實(shí)驗(yàn)結(jié)果與分析 理論計(jì)算結(jié)果為：。從結(jié)果圖中可以看出，matlab計(jì)算結(jié)果與理論計(jì)算結(jié)果基本一致，但是，從誤差分析來(lái)看，matlab計(jì)算結(jié)果比理論計(jì)算結(jié)果要大一些，不過(guò)，這個(gè)誤差非常小，所以，從工程上說(shuō)，這種誤差應(yīng)該是可以忽略不計(jì)的。………………..5離散時(shí)間系統(tǒng)的變換域分析5.1實(shí)驗(yàn)要求用matlab函數(shù)，求系統(tǒng)的零、極點(diǎn)和幅度頻率響應(yīng)和相位響應(yīng)。5.2實(shí)現(xiàn)代碼b=[0.05280.07970.12950.12950.7970.0528];a=[1-1.81072.4947-1.88010.9537-0.2336];[z,p]=tf2zp(b,a);%求系統(tǒng)的