數(shù)字信號處理matlab實驗第三章報告

1、第三章 離散時間信號的頻域分析一、實驗目的通過matlab仿真進一步研究離散時間系統(tǒng)的時域特性。二、實驗要求1 學習并調試本章所給的例子。2 回答書后給出的問題。三、實驗過程Q.2 運行程序P3.1,求離散時間傅里葉變換的實部、虛部以及幅度和相位譜。離散時間傅里葉變換是w的周期函數(shù)嗎?若是周期是多少？描述這四個圖形表示的對稱性。clf;w = -4*pi:8*pi/511:4*pi;num = 2 1;den = 1 -0.6;h = freqz(num, den, w);% Plot the DTFTsubplot(2,2,1) plot(w/pi,real(h);gridtitle(H(e

2、jomega)的實部)xlabel(omega /pi);ylabel(振幅);subplot(2,2,2)plot(w/pi,imag(h);gridtitle(H(ejomega)的虛部)xlabel(omega /pi);ylabel(振幅);subplot(2,2,3)plot(w/pi,abs(h);gridtitle(|H(ejomega)|幅度譜)xlabel(omega /pi);ylabel(振幅);subplot(2,2,4)plot(w/pi,angle(h);gridtitle(H(ejomega)相位譜)xlabel(omega /pi);ylabel(以弧度為單位的

3、相位);由圖得出是周期函數(shù)，周期為2，實部是和幅度譜是實偶對稱，虛部和相位譜是實奇對稱Q.4修改程序P3.1，計算如下有限長序列的離散時間傅里葉變換：gn=1 3 5 7 9 11 13 15 17并重做習題Q.2。討論你的結果。你能解釋相位譜中的跳變嗎？clf;w = -4*pi:8*pi/511:4*pi;num = 1 3 5 7 9 11 13 15 17;den = 1;h = freqz(num, den, w);% Plot the DTFTsubplot(2,2,1)plot(w/pi,real(h);gridtitle(H(ejomega)的實部) xlabel(omega

4、/pi);ylabel(振幅);subplot(2,2,2)plot(w/pi,imag(h);gridtitle(H(ejomega)的虛部)xlabel(omega /pi);ylabel(振幅);subplot(2,2,3)plot(w/pi,abs(h);gridtitle(|H(ejomega)|幅度譜)xlabel(omega /pi);ylabel(振幅);subplot(2,2,4)plot(w/pi,angle(h);gridtitle(H(ejomega)相位譜)xlabel(omega /pi);ylabel(以弧度為單位的相位);1.由圖得出是周期函數(shù)，周期為2，實部是

5、和幅度譜是實偶對稱，虛部和相位譜是實奇對稱。2.要計算一個系統(tǒng)相頻特性，就要用到反正切函數(shù)，計算機中反正切函數(shù)規(guī)定，在一、二象限中的角度為0pi，三四象限的角度為0-pi。若一個角度從0變到2pi，但實際得到的結果是0pi，再由-pi0，在w=pi處發(fā)生跳變。Q.6通過加入合適的注釋語句和程序語句，修改程序P3.2，對程序生成的圖形的兩個軸加標記。哪個參數(shù)控制時移量？clf;w = -pi:2*pi/255:pi; wo = 0.4*pi; D = 10;num = 1 2 3 4 5 6 7 8 9;h1 = freqz(num, 1, w);h2 = freqz(zeros(1,D) nu

6、m, 1, w);subplot(2,2,1)plot(w/pi,abs(h1);gridtitle(原序列的幅度譜) xlabel(omega /pi);ylabel(振幅);subplot(2,2,2)plot(w/pi,abs(h2);gridtitle(時移后序列的幅度譜)xlabel(omega /pi);ylabel(振幅);subplot(2,2,3)plot(w/pi,angle(h1);gridtitle(原序列的相位譜)xlabel(omega /pi);ylabel(以弧度為單位的相位);subplot(2,2,4)plot(w/pi,angle(h2);gridtitl

7、e(時移后序列的相位譜)xlabel(omega /pi);ylabel(以弧度為單位的相位);參數(shù)D控制位移量Q.8選取不同的時移值重做習題Q.6令Q.6中的位移量D為20，程序結果如圖： Q.10通過加入合適的注釋語句和程序語句，修改程序P3.3，對程序生成的圖形的兩個軸加標記。哪個參數(shù)控制頻移量？clf;w = -pi:2*pi/255:pi; wo = 0.4*pi;num1 = 1 3 5 7 9 11 13 15 17;L = length(num1);h1 = freqz(num1, 1, w);n = 0:L-1; num2 = exp(wo*i*n).*num1;h2 = f

8、reqz(num2, 1, w);subplot(2,2,1)plot(w/pi,abs(h1);gridtitle(原序列的幅度譜)xlabel(omega /pi);ylabel(振幅);subplot(2,2,2)plot(w/pi,abs(h2);gridtitle(頻移后序列的幅度譜)xlabel(omega /pi);ylabel(振幅);subplot(2,2,3)plot(w/pi,angle(h1);gridtitle(原序列的相位譜)xlabel(omega /pi);ylabel(以弧度為單位的相位);subplot(2,2,4)plot(w/pi,angle(h2);g

9、ridtitle(頻移后序列的相位譜)xlabel(omega /pi);ylabel(以弧度為單位的相位);參數(shù)wo控制頻移量Q.12選取不同的頻移值重做習題Q.10令Q.10中wo值為，程序結果如圖Q.14通過加入合適的注釋語句和程序語句，修改程序P3.4，對程序生成的圖形的兩個軸加標記。clf;w = -pi:2*pi/255:pi;x1 = 1 3 5 7 9 11 13 15 17;x2 = 1 -2 3 -2 1;y = conv(x1,x2);h1 = freqz(x1, 1, w);h2 = freqz(x2, 1, w);hp = h1.*h2;h3 = freqz(y,1,

10、w); subplot(2,2,1)plot(w/pi,abs(hp);gridtitle(幅度譜的乘積)xlabel(omega /pi);ylabel(振幅);subplot(2,2,2)plot(w/pi,abs(h3);gridtitle(卷積后序列的幅度譜)xlabel(omega /pi);ylabel(振幅);subplot(2,2,3)plot(w/pi,angle(hp);gridtitle(相位譜的和)xlabel(omega /pi);ylabel(以弧度為單位的相位);subplot(2,2,4)plot(w/pi,angle(h3);gridtitle(卷積后序列的相

11、位譜)xlabel(omega /pi);ylabel(以弧度為單位的相位);Q.16選取兩個改變了長度的序列，重做習題Q.14更改Q.14中的兩個序列為X1=1 1 2 4,X2=4 2 3 1,程序結果如圖：Q.18通過加入合適的注釋語句和程序語句，修改程序P3.5，對程序生成的圖形的兩個軸加標記。clf;w = -pi:2*pi/255:pi;x1 = 1 3 5 7 9 11 13 15 17;x2 = 1 -1 1 -1 1 -1 1 -1 1;y = x1.*x2;h1 = freqz(x1, 1, w);h2 = freqz(x2, 1, w);h3 = freqz(y,1,w)

12、; subplot(3,1,1)plot(w/pi,abs(h1);gridtitle(第一個序列的幅度譜)xlabel(omega /pi);ylabel(振幅);subplot(3,1,2)plot(w/pi,abs(h2);gridtitle(第二個序列的幅度譜)xlabel(omega /pi);ylabel(振幅);subplot(3,1,3)plot(w/pi,abs(h3);gridtitle(乘積序列的幅度譜)xlabel(omega /pi);ylabel(振幅);Q.20通過加入合適的注釋語句和程序語句，修改程序P3.6，對程序生成的圖形的兩個軸加標記。試解釋程序怎樣進行時

13、間反轉運算。clf;w = -pi:2*pi/255:pi;num = 1 2 3 4;L = length(num)-1;h1 = freqz(num, 1, w);h2 = freqz(fliplr(num), 1, w); h3 = exp(w*L*i).*h2;subplot(2,2,1)plot(w/pi,abs(h1);gridtitle(原序列的幅度譜)xlabel(omega /pi);ylabel(振幅);subplot(2,2,2)plot(w/pi,abs(h3);gridtitle(時間反轉后序列的幅度譜)xlabel(omega /pi);ylabel(振幅);sub

14、plot(2,2,3)plot(w/pi,angle(h1);gridtitle(原序列的相位譜)xlabel(omega /pi);ylabel(以弧度為單位的相位);subplot(2,2,4)plot(w/pi,angle(h3);gridtitle(時間反轉后序列的相位譜)xlabel(omega /pi);ylabel(以弧度為單位的相位);程序先將num沿時間原點翻轉，再進行離散時間傅里葉變換，然后變換后的函數(shù)乘以時移分量exp(w*L*i)。Q.22選取兩個不同長度的序列，重做習題Q.201.將Q.20中的序列改為num=1 1 2 4 5,程序結果如圖：2. 將Q.20中的序列

15、改為num=1 6 8 3 1 2,程序結果如圖： Q.24編寫一個MATLAB程序，用一個N點復數(shù)離散傅里葉變換計算兩個長度為N的實數(shù)序列的N點離散傅里葉變換，并將結果同直接使用兩個N點離散傅里葉變換得到的結果進行比較。clf;g=1 1 2 4;h=2 3 2 1;x=g+i*h;N = length(x)-1; n = 0:N;gk=fft(g);hk=fft(h);xk=fft(x);xk1=fft(conj(x);gk1=(xk+xk1)/2;hk1=(xk-xk1)/2i;subplot(4,2,1)stem(n,abs(gk);gridtitle(實部序列gk的離散傅里葉變換的幅

16、度)xlabel(時間序號n);ylabel(振幅);subplot(4,2,2)stem(n,abs(hk);gridtitle(虛部序列gk的離散傅里葉變換的幅度)xlabel(時間序號n);ylabel(振幅);subplot(4,2,3)stem(n,abs(gk1);gridtitle(通過xk得到的gk1的離散傅里葉變換的幅度)xlabel(時間序號n);ylabel(振幅);subplot(4,2,4)stem(n,abs(hk1);gridtitle(通過xk得到的hk1的離散傅里葉變換的幅度)xlabel(時間序號n);ylabel(振幅);subplot(4,2,5)ste

17、m(n,angle(gk);gridtitle(實部序列gk的離散傅里葉變換的相位)xlabel(時間序號n);ylabel(以弧度為單位的相位);subplot(4,2,6)stem(n,angle(hk);gridtitle(虛部序列hk的離散傅里葉變換的相位)xlabel(時間序號n);ylabel(以弧度為單位的相位);subplot(4,2,7)stem(n,angle(gk1);gridtitle(通過xk得到的gk1的離散傅里葉變換的相位)xlabel(時間序號n);ylabel(以弧度為單位的相位);subplot(4,2,8)stem(n,angle(hk1);gridtit

18、le(通過xk得到的hk1的離散傅里葉變換的相位)xlabel(時間序號n);ylabel(以弧度為單位的相位);程序運行結果如下，說明兩種方法得到的結果一樣：Q.26在函數(shù)circshift中，命令rem的作用是什么？函數(shù)cirshift中rem是用來確定要補的零的個數(shù)。Q.28在函數(shù)circonv中，運算符=的作用是什么？函數(shù)circonv中=是用來確定兩個序列的長度是否相等。Q.30通過加入合適的注釋語句和程序語句，修改程序P3.7，對程序生成的圖形的兩個軸加標記。哪一個參數(shù)決定時移量？若時移量大于序列長度，將會發(fā)生什么？clf;M = 6;a = 0 1 2 3 4 5 6 7 8 9

19、;b = circshift(a,M);L = length(a)-1;n = 0:L;subplot(2,1,1);stem(n,a);axis(0,L,min(a),max(a);title(原序列)xlabel(時間序號n);ylabel(a);subplot(2,1,2);stem(n,b);axis(0,L,min(a),max(a);title(圓周移位,num2str(M),個樣本得到的序列)xlabel(時間序號n);ylabel(a);M控制時移量，若時移量大于序列長度會得到時移M1=M-序列長度 的序列 下圖是M=6時的結果，此時時移量M小于序列長度： 下圖是M=13時的結

20、果，此時時移量M大于序列長度： Q.32通過加入合適的注釋語句和程序語句，修改程序P3.8，對程序生成的圖形的兩個軸加標記。clf;x = 0 2 4 6 8 10 12 14 16;N = length(x)-1; n = 0:N;y = circshift(x,5);XF = fft(x);YF = fft(y);subplot(2,2,1)stem(n,abs(XF);gridtitle(原序列的離散傅里葉變換的幅度);xlabel(時間序號n);ylabel(振幅);subplot(2,2,2)stem(n,abs(YF);gridtitle(圓周移位后的序列的離散傅里葉變換的幅度);

21、xlabel(時間序號n);ylabel(振幅);subplot(2,2,3)stem(n,angle(XF);gridtitle(原序列的離散傅里葉變換的相位);xlabel(時間序號n);ylabel(以弧度為單位的相位);subplot(2,2,4)stem(n,angle(YF);gridtitle(圓周移位后的序列的離散傅里葉變換的相位);xlabel(時間序號n);ylabel(以弧度為單位的相位);時移量是5.結果如下圖：Q.34選取兩個不同的時移量，重做習題Q.321.時移量為10，結果如下圖：2.時移量為15，結果如下圖：Q.36運行程序P3.9并驗證離散傅里葉變換的圓周卷積

22、性質clf;g1 = 1 2 3 4 5 6; g2 = 1 -2 3 3 -2 1;ycir = circonv(g1,g2);disp(圓周卷積的結果 = );disp(ycir)G1 = fft(g1); G2 = fft(g2);yc = real(ifft(G1.*G2);disp(離散傅里葉變換乘積的離散傅里葉逆變換的結果 = );disp(yc)結果如下圖：兩個結果一致，說明兩個序列循環(huán)卷積等于這兩個序列分別的DFT相乘后的IDFTQ.38運行程序P3.10并驗證線性卷積可通過圓周卷積得到clf;g1 = 1 2 3 4 5;g2 = 2 2 0 1 1;g1e = g1 zer

23、os(1,length(g2)-1);g2e = g2 zeros(1,length(g1)-1);ylin = circonv(g1e,g2e);disp(通過圓周卷積的線性卷積 = );disp(ylin);y = conv(g1, g2);disp(直接線性卷積 = );disp(y)結果如下圖：兩個結果一致，說明線性卷積可通過圓周卷積得到Q.40編寫一個MATLAB程序，對兩個序列分別做離散傅里葉變換，以生成他們的線性卷積。用此程序驗證習題Q.38的結果clf;g1 = 1 2 3 4 5;g2 = 2 2 0 1 1;g1e = g1 zeros(1,length(g2)-1);g2

24、e = g2 zeros(1,length(g1)-1);G1e = fft(g1e); G2e = fft(g2e);ylin = real(ifft(G1e.*G2e);disp(通過離散傅里葉變換得到的線性卷積 = );disp(ylin);y = conv(g1, g2);disp(直接線性卷積 = );disp(y)結果如下圖：兩者結果一致，驗證了習題Q.38的結果Q.42運行程序P3.11。由于周期序列的偶數(shù)部分的離散傅里葉變換是原序列的XEF的實數(shù)部分，XEF的虛部應該為零。你能驗證它們嗎？你怎樣解釋仿真結果？x = 1 2 4 2 6 32 6 4 2 zeros(1,247)

25、;x1 = x(1) x(256:-1:2);xe = 0.5 *(x + x1);XF = fft(x); XEF = fft(xe);clf;k = 0:255;subplot(2,2,1);plot(k/128,real(XF); grid;ylabel(振幅);title(Re(DFTxn);subplot(2,2,2);plot(k/128,imag(XF); grid;ylabel(振幅);title(Im(DFTxn);subplot(2,2,3)plot(k/128,real(XEF); grid;xlabel(時間序號n);ylabel(振幅);title(Re(DFTx_e

26、n);subplot(2,2,4);plot(k/128,imag(XEF); grid;xlabel(時間序號n);ylabel(振幅);title(Im(DFTx_en);Q.44運行程序P3.12。你得到的a和b 的值相等嗎？clf;x = (1:128) (128:-1:1);XF = fft(x);a = sum(x.*x)b = round(sum(abs(XF).2)/256)運行結果如圖：a和b的值相等Q.46使用程序P3.1在單位圓上求下面的z變換：Gz=2+5z-1+9z-2+5z-3+3z-45+45z-1+2z-2+z-3+z-4clf;w = -pi:2*pi/511:pi;num = 2 5 9 5 3;den = 5 45 2 1 1;h = f