數(shù)字信號(hào)處理第三章

1、數(shù)字信號(hào)處理第三章實(shí)驗(yàn)程序3.1計(jì)算離散時(shí)間傅里葉變換% Program P3_1% Evaluation of the DTFT clf;% Compute the frequency samples of the DTFTw = -4*pi:8*pi/511:4*pi;num = 2 1;den = 1 -0.6;h = freqz(num, den, w);% Plot the DTFTsubplot(2,1,1)plot(w/pi,real(h);gridtitle('Real part of H(ejomega)')xlabel('omega /pi'

2、);ylabel('Amplitude');subplot(2,1,2)plot(w/pi,imag(h);gridtitle('Imaginary part of H(ejomega)')xlabel('omega /pi');ylabel('Amplitude');pausesubplot(2,1,1)plot(w/pi,abs(h);gridtitle('Magnitude Spectrum |H(ejomega)|')xlabel('omega /pi');ylabel('Ampli

3、tude');subplot(2,1,2)plot(w/pi,angle(h);gridtitle('Phase Spectrum argH(ejomega)')xlabel('omega /pi');ylabel('Phase in radians');Q3.1離散時(shí)間傅里葉變換的原始序列是H(ejw)=(2+z-1)/(1-0.6z-1)。Pause的作用是暫停等待用戶輸入任意鍵后接著執(zhí)行以下命令。Q3.2 是周期函數(shù)，周期是2。實(shí)部和幅度譜是關(guān)于y軸對(duì)稱，是偶函數(shù)；虛部和相位譜是關(guān)于原點(diǎn)對(duì)稱，是奇函數(shù)。Q3.3clf;N = 512

4、;num = 0.7 -0.5 0.3 1;den = 1 0.3 -0.5 0.7;h,w = freqz(num, den, N);subplot(2,1,1)plot(w/pi,real(h);gridtitle('Real part of H(ejomega)')xlabel('omega /pi');ylabel('Amplitude');subplot(2,1,2)plot(w/pi,imag(h);gridtitle('Imaginary part of H(ejomega)')xlabel('omega /

5、pi');ylabel('Amplitude');pausesubplot(2,1,1)plot(w/pi,abs(h);gridtitle('Magnitude Spectrum |H(ejomega)|')xlabel('omega /pi');ylabel('Amplitude');subplot(2,1,2)plot(w/pi,angle(h);gridtitle('Phase Spectrum argH(ejomega)')xlabel('omega /pi');ylabel(&#

6、39;Phase in radians'); 還是周期函數(shù)，周期是2。相位譜的跳變的原因是：在利用反正切函數(shù)計(jì)算角度的時(shí)候，其中的一個(gè)分支出現(xiàn)了衰減，造成了跳變。clf;N = 512;num = 0.7 -0.5 0.3 1;den = 1 0.3 -0.5 0.7;h,w = freqz(num, den, N);subplot(2,1,1)plot(w/pi,unwrap(angle(h);gridtitle('Phase Spectrum argH(ejomega)')xlabel('omega /pi');ylabel('Phase i

7、n radians');Q3.4 修改后的程序?yàn)閏lf;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,1,1)plot(w/pi,real(h);gridtitle('Real part of H(ejomega)')xlabel('omega /pi');ylabel('Amplitude');subplot(2,1,2)plot(w/pi,imag(h);gridti

8、tle('Imaginary part of H(ejomega)')xlabel('omega /pi');ylabel('Amplitude');pausesubplot(2,1,1)plot(w/pi,abs(h);gridtitle('Magnitude Spectrum |H(ejomega)|')xlabel('omega /pi');ylabel('Amplitude');subplot(2,1,2)plot(w/pi,angle(h);gridtitle('Phase Spe

9、ctrum argH(ejomega)')xlabel('omega /pi');ylabel('Phase in radians');w 是周期函數(shù)，周期是2。實(shí)部和幅度譜是關(guān)于y軸對(duì)稱，是偶函數(shù)；虛部和相位譜是關(guān)于原點(diǎn)對(duì)稱，是奇函數(shù)。Q3.5若要改為以度為單位，則將程序中的第二個(gè)圖的程序改為subplot(2,1,2)plot(w/pi,180*angle(h)/pi);gridtitle('Phase Spectrum argH(ejomega)')xlabel('omega /pi');ylabel('Ph

10、ase in degrees');就可以了。-3.2離散時(shí)間傅里葉變換的性質(zhì)1. 時(shí)移特性clf;w = -pi:2*pi/255:pi; D = 10; num = 1 2 3 4 5 6 7 8 9;h1 = freqz(num, 1, w);h2 = freqz(zeros(1,D) num, 1, w);subplot(2,2,1)plot(w/pi,abs(h1);gridtitle('Magnitude Spectrum of Original Sequence','FontSize',8)xlabel('omega /pi'

11、);ylabel('Amplitude');subplot(2,2,2)plot(w/pi,abs(h2);gridtitle('Magnitude Spectrum of Time-Shifted Sequence','FontSize',8)xlabel('omega /pi');ylabel('Amplitude');subplot(2,2,3)plot(w/pi,angle(h1);gridtitle('Phase Spectrum of Original Sequence','Fo

12、ntSize',8)xlabel('omega /pi');ylabel('Phase in radians');subplot(2,2,4)plot(w/pi,angle(h2);gridtitle('Phase Spectrum of Time-Shifted Sequence','FontSize',8)xlabel('omega /pi');ylabel('Phase in radians');Q3.6參數(shù)D控制時(shí)移量。Q3.7 D=10 D=50 時(shí)移特性：信號(hào)在時(shí)域移動(dòng)某個(gè)距離，

13、則所得信號(hào)的幅度譜和原信號(hào)相同，而相位譜是原信號(hào)的相位譜再附加一個(gè)線性相移，由時(shí)移特性可以看到，信號(hào)的相位譜可以反映信號(hào)在時(shí)域中的位置信息，不同位置上的同一信號(hào)，它們具有不同的相頻特性，而幅頻特性相同。Q3.8如上圖所示Q3.9改變序列長(zhǎng)度num = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 2425 26 27 28 29;所得的圖像為 D=10 D=50 從上圖中可以看出，增加序列的長(zhǎng)度，使得幅度譜更加窄，而相位譜則更加密集和陡峭。2. 平移特性Q3.10clf;w = -pi:2*pi/255:pi; wo =

14、 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 = freqz(num2, 1, w);subplot(2,2,1)plot(w/pi,abs(h1);gridtitle('Magnitude Spectrum of Original Sequence','FontSize',8)xlabel('omega /pi');ylabel('Amplitude&#

15、39;);subplot(2,2,2)plot(w/pi,abs(h2);gridtitle('Magnitude Spectrum of Frequency-Shifted Sequence','FontSize',8)xlabel('omega /pi');ylabel('Amplitude');subplot(2,2,3)plot(w/pi,angle(h1);gridtitle('Phase Spectrum of Original Sequence','FontSize',8)xlabel

16、('omega /pi');ylabel('Phase in radians');subplot(2,2,4)plot(w/pi,angle(h2);gridtitle('Phase Spectrum of Frequency-Shifted Sequence','FontSize',8)xlabel('omega /pi');ylabel('Phase in radians');Wo控制平移量。Q3.11由結(jié)果圖Q3.11可得出在參數(shù)wo的控制下，離散時(shí)間傅里葉變換的幅度譜和相位譜都隨著控制參數(shù)右

17、移k個(gè)單位（wo=k*pi）。 k=0.4 k=-0.4 Q3.12將k改為-0.4得到的運(yùn)行結(jié)果如上圖。Q3.13改變序列長(zhǎng)度序列：num1=1 3 5 7 9 11 13 15 17 19 21 23 25 27 29序列：num2=11 13 15 17 19 21 23 25 27 29 31 33 35 37 39;3. 卷積性質(zhì)Q3.14clf;w = -pi:2*pi/255:pi; % freqency vector for evaluating DTFTx1 = 1 3 5 7 9 11 13 15 17; x2 = 1 -2 3 -2 1;y = conv(x1,x2);h

18、1 = freqz(x1, 1, w);h2 = freqz(x2, 1, w);hp = h1.*h2;h3 = freqz(y,1,w);subplot(2,2,1)plot(w/pi,abs(hp);gridtitle('Product of Magnitude Spectra','FontSize',8)xlabel('omega /pi');ylabel('Amplitude');subplot(2,2,2)plot(w/pi,abs(h3);gridtitle('Magnitude Spectrum of Co

19、nvolved Sequence','FontSize',8)xlabel('omega /pi');ylabel('Amplitude');subplot(2,2,3)plot(w/pi,angle(hp);gridtitle('Sum of Phase Spectra','FontSize',8)xlabel('omega /pi');ylabel('Phase in radians');subplot(2,2,4)plot(w/pi,angle(h3);gridtitl

20、e('Phase Spectrum of Convolved Sequence','FontSize',8)xlabel('omega /pi');ylabel('Phase in radians');Q3.15分析結(jié)果圖可以得出幅度譜的乘積和卷積后的幅度譜相同，相位譜的乘積和卷積后的相位譜相同。 Q3.16 x1 = 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33; x2 = 1 -2 3 -2 1 -5 2 -3 1;運(yùn)行結(jié)果如上邊第二個(gè)圖所示。4. 調(diào)制性質(zhì)Q3.17clf;w =

21、 -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); subplot(3,1,1)plot(w/pi,abs(h1);gridtitle('Magnitude Spectrum of First Sequence')xlabel('omega /pi');ylabel('Amplitude');subplot(

22、3,1,2)plot(w/pi,abs(h2);gridtitle('Magnitude Spectrum of Second Sequence')xlabel('omega /pi');ylabel('Amplitude');subplot(3,1,3)plot(w/pi,abs(h3);gridtitle('Magnitude Spectrum of Product Sequence')xlabel('omega /pi');ylabel('Amplitude');Q3.18分析圖得出乘積序列的

23、幅度譜近似等于兩序列的幅度譜的和. Q3.19將序列改變?yōu)閤1 = 1 3 5 7 9 11 13 15 17 19 21 23 25 27; x2 = 1 -1 1 -1 1 -1 1 -1 1 0 2 -4 7 -1得到的運(yùn)行結(jié)果為上右圖。乘積序列的幅度譜依然近似等于兩序列的幅度譜的和.5. 時(shí)間反轉(zhuǎn)性質(zhì)Q3.20clf;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

24、,1)plot(w/pi,abs(h1);gridtitle('Magnitude Spectrum of Original Sequence','FontSize',8)xlabel('omega /pi');ylabel('Amplitude');subplot(2,2,2)plot(w/pi,abs(h3);gridtitle('Magnitude Spectrum of Time-Reversed Sequence','FontSize',8)xlabel('omega /pi

25、9;);ylabel('Amplitude');subplot(2,2,3)plot(w/pi,angle(h1);gridtitle('Phase Spectrum of Original Sequence','FontSize',8)xlabel('omega /pi');ylabel('Phase in radians');subplot(2,2,4)plot(w/pi,angle(h3);gridtitle('Phase Spectrum of Time-Reversed Sequence'

26、,'FontSize',8)xlabel('omega /pi');ylabel('Phase in radians');Q3.21分析圖得出序列的幅度譜隨時(shí)間反轉(zhuǎn)不發(fā)生變化，序列相位譜隨時(shí)間反轉(zhuǎn)而反轉(zhuǎn)180。 Q3.22改變序列長(zhǎng)度num = 1 -2 3 -4 5 -6 7 -8 ;得到的運(yùn)行結(jié)果為上右，結(jié)果依然是序列的幅度譜隨時(shí)間反轉(zhuǎn)不發(fā)生變化，序列相位譜隨時(shí)間反轉(zhuǎn)而反轉(zhuǎn)180。3.5離散傅里葉變換和離散傅里葉逆變換的運(yùn)算Q3.23clf;N=200; L=256; nn = 0:N-1;kk = 0:L-1;xR = 0.1*(1:100

27、) zeros(1,N-100); xI = zeros(1,N); x = xR + i*xI;XF = fft(x,L);subplot(3,2,1);grid;plot(nn,xR);grid;title('Rexn');xlabel('Time index n');ylabel('Amplitude');subplot(3,2,2);plot(nn,xI);grid;title('Imxn');xlabel('Time index n');ylabel('Amplitude');subplo

28、t(3,2,3);plot(kk,real(XF);grid;title('ReXk');xlabel('Frequency index k');ylabel('Amplitude');subplot(3,2,4);plot(kk,imag(XF);grid;title('ImXk');xlabel('Frequency index k');ylabel('Amplitude');xx = ifft(XF,L);subplot(3,2,5);plot(kk,real(xx);grid;title(&

29、#39;Real part of IDFTXk');xlabel('Time index n');ylabel('Amplitude');subplot(3,2,6);plot(kk,imag(xx);grid;title('Imag part of IDFTXk');xlabel('Time index n');ylabel('Amplitude');Q3.24clf;N=256;nn = 0:N-1;ntime = -N/2:N/2-1;g = (0.75).abs(ntime); h = (-0.9)

30、.ntime; GF = fft(g);HF = fft(h);x = g + i*h; XF = fft(x);XFstar = conj(XF);XFstarmod = XFstar(1) fliplr(XFstar(2:N);GF2 = 0.5*(XF + XFstarmod);HF2 = -i*0.5*(XF - XFstarmod);abs(max(GF-GF2)abs(max(HF-HF2)figure(1);clf;subplot(2,2,1);grid;plot(nn,real(GF);grid;title('Two N-point DFT''s'

31、;);xlabel('Frequency index k');ylabel('ReGk');subplot(2,2,2);plot(nn,imag(GF);grid;title('Two N-point DFT''s');xlabel('Frequency index k');ylabel('ImGk');subplot(2,2,3);grid;plot(nn,real(GF2);grid;title('Single N-point DFT');xlabel('Frequen

32、cy index k');ylabel('ReGk');subplot(2,2,4);plot(nn,imag(GF2);grid;title('Single N-point DFT');xlabel('Frequency index k');ylabel('ImGk');figure(2);clf;subplot(2,2,1);grid;plot(nn,real(HF);grid;title('Two N-point DFT''s');xlabel('Freq index k

33、9;);ylabel('ReHk');subplot(2,2,2);plot(nn,imag(HF);grid;title('Two N-point DFT''s');xlabel('Freq index k');ylabel('ImHk');subplot(2,2,3);grid;plot(nn,real(HF2);grid;title('Single N-point DFT');xlabel('Freq index k');ylabel('ReHk');subpl

34、ot(2,2,4);plot(nn,imag(HF2);grid;title('Single N-point DFT');xlabel('Freq index k');ylabel('ImHk'); Q3.25clf;N = 128; TwoN = 2*N;W2N = exp(-i*pi/N);k = 0:TwoN-1;v = (-0.7.k);g = downsample(v,2); h = downsample(v,2,1); x = g + i*h;XF = fft(x); XFstar = conj(XF);XFstarmod = XFs

35、tar(1) fliplr(XFstar(2:N);GF = 0.5*(XF + XFstarmod);HF = -i*0.5*(XF - XFstarmod);VF = GF GF + (W2N.k).*HF HF;VF2 = fft(v);abs(max(VF-VF2)subplot(2,2,1);plot(k,real(VF);grid;title('Complex N-point DFT');xlabel('Frequency index k');ylabel('ReVk');subplot(2,2,2);plot(k,imag(VF);

36、grid;title('Complex N-point DFT');xlabel('Frequency index k');ylabel('ImVk');subplot(2,2,3);plot(k,real(VF2);grid;title('Real 2N-point DFT');xlabel('Frequency index k');ylabel('ReVk');subplot(2,2,4);plot(k,imag(VF2);grid;title('Real 2N-point DFT

37、9;);xlabel('Frequency index k');ylabel('ImVk'); 3.4離散傅里葉函數(shù)的性質(zhì)Q3.26rem(x,y),x是除y以后剩余部分。Q3.27輸入序列x循環(huán)移位留下的位置。如果M > 0,那么circshift刪除左邊的元素向量x和附加他們右側(cè)獲得剩下的元素循環(huán)轉(zhuǎn)移序列。如果如果M < 0,然后circshift第一次補(bǔ)充的x的長(zhǎng)度,即。,最右邊的長(zhǎng)度(x)- m樣品從x和附加右邊的樣品得到循環(huán)轉(zhuǎn)移序列。Q3.28這是二元關(guān)系不等于操作符。 = B返回值1如果A和B是不平等的值0如果A和B都是平等的。Q3.29

38、輸入是平等的兩個(gè)向量x1和x2長(zhǎng)度l .理解circonv是如何工作的,它是有用的定期x2的延伸。讓x2p x2的無(wú)限長(zhǎng)的周期延長(zhǎng)。從概念上講,常規(guī)時(shí)間逆轉(zhuǎn)x2p和集x2tr 1到L等于元素的時(shí)間逆轉(zhuǎn)x2p版本。元素1通過(guò)y L的輸出向量然后通過(guò)x1和長(zhǎng)度之間的內(nèi)積向量sh循環(huán)變化對(duì)時(shí)間逆轉(zhuǎn)向量x2tr。對(duì)于輸出樣例yn,1nL、正確的循環(huán)移位是n - 1的位置。Q3.30clf;M = 6;a = 0 1 2 3 4 5 6 7 8 9;b = circshift(a,M); L = length(a)-1;n = 0:L;subplot(2,1,1);stem(n,a);axis(0,L,

39、min(a),max(a);title('Original Sequence');xlabel('time index n');ylabel('an');subplot(2,1,2);stem(n,b);axis(0,L,min(a),max(a);title('Sequence Obtained by Circularly Shifting by ',num2str(M),'Samples');xlabel('time index n');ylabel('bn')M值決定時(shí)移量。Q

40、3.31Q3.32clf;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);grid;title('Magnitude of DFT of Original Sequence');xlabel('Frequency index k');ylabel('|Xk|');subplot(2,2,2);stem(n,abs(YF);grid;title(

41、'Magnitude of DFT of Circularly Shifted Sequence');xlabel('Frequency index k');ylabel('|Yk|');subplot(2,2,3);stem(n,angle(XF);grid;title('Phase of DFT of Original Sequence');xlabel('Frequency index k');ylabel('arg(Xk)');subplot(2,2,4);stem(n,angle(YF);

42、grid;title('Phase of DFT of Circularly Shifted Sequence');xlabel('Frequency index k');ylabel('arg(Yk)');時(shí)移量是8.Q3.33 Q3.34 M=5 運(yùn)行結(jié)果如上右圖所示。Q3.35Length = 13 Length = 20 Q3.36 g1 = 1 2 3 4 5 6; g2 = 1 -2 3 3 -2 1;ycir = circonv(g1,g2);disp('Result of circular convolution = &#

43、39;);disp(ycir)G1 = fft(g1); G2 = fft(g2);yc = real(ifft(G1.*G2);disp('Result of IDFT of the DFT products = ');disp(yc)運(yùn)行結(jié)果Q3.37結(jié)果如下：Q3.38g1 = 1 2 3 4 5;g2 = 2 2 0 1 1;g1e = g1 zeros(1,length(g2)-1);g2e = g2 zeros(1,length(g1)-1);ylin = circonv(g1e,g2e);disp('Linear convolution via circu

44、lar convolution = ');disp(ylin);y = conv(g1, g2);disp('Direct linear convolution = ');disp(y)結(jié)果如下：Q3.39 g1 = 3 1 4 1 5 9 2;g2 = 1 1 1 0 0; g1 = 5 4 3 2 1 0;g2 = -2 1 2 3 4;Q3.40g1 = 1 2 3 4 5;g2 = 2 2 0 1 1;g1e = g1 zeros(1,length(g2)-1);g2e = g2 zeros(1,length(g1)-1);G1EF = fft(g1e);G2EF = fft(g2e);ylin = real(ifft(G1EF.*G2EF);disp('直線線性卷積 = ' );disp(ylin)；Q3.41x = 1 2 4 2 6 32 6 4 2 zeros(1,247);