数字信号处理第三章

1、数字信号处理第三章实验程序3.1计算离散时间傅里叶变换% 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离散时间傅里叶变换的原始序列是H(ejw)=(2+z-1)/(1-0.6z-1)。Pause的作用是暂停等待用户输入任意键后接着执行以下命令。Q3.2 是周期函数，周期是2。实部和幅度谱是关于y轴对称，是偶函数；虚部和相位谱是关于原点对称，是奇函数。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'); 还是周期函数，周期是2。相位谱的跳变的原因是：在利用反正切函数计算角度的时候，其中的一个分支出现了衰减，造成了跳变。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 修改后的程序为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,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 是周期函数，周期是2。实部和幅度谱是关于y轴对称，是偶函数；虚部和相位谱是关于原点对称，是奇函数。Q3.5若要改为以度为单位，则将程序中的第二个图的程序改为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离散时间傅里叶变换的性质1. 时移特性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参数D控制时移量。Q3.7 D=10 D=50 时移特性：信号在时域移动某个距离，

13、则所得信号的幅度谱和原信号相同，而相位谱是原信号的相位谱再附加一个线性相移，由时移特性可以看到，信号的相位谱可以反映信号在时域中的位置信息，不同位置上的同一信号，它们具有不同的相频特性，而幅频特性相同。Q3.8如上图所示Q3.9改变序列长度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 从上图中可以看出，增加序列的长度，使得幅度谱更加窄，而相位谱则更加密集和陡峭。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由结果图Q3.11可得出在参数wo的控制下，离散时间傅里叶变换的幅度谱和相位谱都随着控制参数右

17、移k个单位（wo=k*pi）。 k=0.4 k=-0.4 Q3.12将k改为-0.4得到的运行结果如上图。Q3.13改变序列长度序列：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. 卷积性质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分析结果图可以得出幅度谱的乘积和卷积后的幅度谱相同，相位谱的乘积和卷积后的相位谱相同。 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;运行结果如上边第二个图所示。4. 调制性质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将序列改变为x1 = 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得到的运行结果为上右图。乘积序列的幅度谱依然近似等于两序列的幅度谱的和.5. 时间反转性质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分析图得出序列的幅度谱随时间反转不发生变化，序列相位谱随时间反转而反转180。 Q3.22改变序列长度num = 1 -2 3 -4 5 -6 7 -8 ;得到的运行结果为上右，结果依然是序列的幅度谱随时间反转不发生变化，序列相位谱随时间反转而反转180。3.5离散傅里叶变换和离散傅里叶逆变换的运算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离散傅里叶函数的性质Q3.26rem(x,y),x是除y以后剩余部分。Q3.27输入序列x循环移位留下的位置。如果M > 0,那么circshift删除左边的元素向量x和附加他们右侧获得剩下的元素循环转移序列。如果如果M < 0,然后circshift第一次补充的x的长度,即。,最右边的长度(x)- m样品从x和附加右边的样品得到循环转移序列。Q3.28这是二元关系不等于操作符。 = B返回值1如果A和B是不平等的值0如果A和B都是平等的。Q3.29

38、输入是平等的两个向量x1和x2长度l .理解circonv是如何工作的,它是有用的定期x2的延伸。让x2p x2的无限长的周期延长。从概念上讲,常规时间逆转x2p和集x2tr 1到L等于元素的时间逆转x2p版本。元素1通过y L的输出向量然后通过x1和长度之间的内积向量sh循环变化对时间逆转向量x2tr。对于输出样例yn,1nL、正确的循环移位是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值决定时移量。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)');时移量是8.Q3.33 Q3.34 M=5 运行结果如上右图所示。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)运行结果Q3.37结果如下：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)结果如下：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);