数字信号处理第三章

数字信号处理第三章实验程序 3.1计算离散时间傅里叶变换 % Program P3_1 % Evaluation of the DTFT clf; % Compute the frequency samples of the DTFT w = -4*pi:8*pi/511:4*pi; num = [2 1];den = [1 -0.6]; h = freqz(num, den, w); % Plot the DTFT subplot(2,1,1) plot(w/pi,real(h));grid title('Real part of H(e^{j\omega})') xlabel('\omega /\pi'); ylabel('Amplitude'); subplot(2,1,2) plot(w/pi,imag(h));grid title('Imaginary part of H(e^{j\omega})') xlabel('\omega /\pi'); ylabel('Amplitude'); pause subplot(2,1,1) plot(w/pi,abs(h));grid title('Magnitude Spectrum |H(e^{j\omega})|') xlabel('\omega /\pi'); ylabel('Amplitude'); subplot(2,1,2) plot(w/pi,angle(h));grid title('Phase Spectrum arg[H(e^{j\omega})]') xlabel('\omega /\pi'); ylabel('Phase in radians'); Q3.1离散时间傅里叶变换的原始序列是H(e^jw)=(2+z^-1)/(1-0.6z^-1)。Pause的作用是暂停等待用户输入任意键后接着执行以下命令。 Q3.2 是周期函数，周期是2π。实部和幅度谱是关于y轴对称，是偶函数；虚部和相位谱是关于原点对称，是奇函数。 Q3.3 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,real(h));grid title('Real part of H(e^{j\omega})') xlabel('\omega /\pi'); ylabel('Amplitude'); subplot(2,1,2) plot(w/pi,imag(h));grid title('Imaginary part of H(e^{j\omega})') xlabel('\omega /\pi'); ylabel('Amplitude'); pause subplot(2,1,1) plot(w/pi,abs(h));grid title('Magnitude Spectrum |H(e^{j\omega})|') xlabel('\omega /\pi'); ylabel('Amplitude'); subplot(2,1,2) plot(w/pi,angle(h));grid title('Phase Spectrum arg[H(e^{j\omega})]') xlabel('\omega /\pi'); ylabel('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)));grid title('Phase Spectrum arg[H(e^{j\omega})]') xlabel('\omega /\pi'); ylabel('Phase in 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 DTFT subplot(2,1,1) plot(w/pi,real(h));grid title('Real part of H(e^{j\omega})') xlabel('\omega /\pi'); ylabel('Amplitude'); subplot(2,1,2) plot(w/pi,imag(h));grid title('Imaginary part of H(e^{j\omega})') xlabel('\omega /\pi'); ylabel('Amplitude'); pause subplot(2,1,1) plot(w/pi,abs(h));grid title('Magnitude Spectrum |H(e^{j\omega})|') xlabel('\omega /\pi'); ylabel('Amplitude'); subplot(2,1,2) plot(w/pi,angle(h));grid title('Phase Spectrum arg[H(e^{j\omega})]') xlabel('\omega /\pi'); ylabel('Phase in radians'); 是周期函数，周期是2π。实部和幅度谱是关于y轴对称，是偶函数；虚部和相位谱是关于原点对称，是奇函数。 Q3.5若要改为以度为单位，则将程序中的第二个图的程序改为 subplot(2,1,2) plot(w/pi,180*angle(h)/pi);grid title('Phase Spectrum arg[H(e^{j\omega})]') xlabel('\omega /\pi'); ylabel('Phase 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));grid title('Magnitude Spectrum of Original Sequence','FontSize',8) xlabel('\omega /\pi'); ylabel('Amplitude'); subplot(2,2,2) plot(w/pi,abs(h2));grid title('Magnitude Spectrum of Time-Shifted Sequence','FontSize',8) xlabel('\omega /\pi'); ylabel('Amplitude'); subplot(2,2,3) plot(w/pi,angle(h1));grid title('Phase Spectrum of Original Sequence','FontSize',8) xlabel('\omega /\pi'); ylabel('Phase in radians'); subplot(2,2,4) plot(w/pi,angle(h2));grid title('Phase Spectrum of Time-Shifted Sequence','FontSize',8) xlabel('\omega /\pi'); ylabel('Phase in radians'); Q3.6参数D控制时移量。 Q3.7 D=10                                      D=50 时移特性：信号在时域移动某个距离，则所得信号的幅度谱和原信号相同，而相位谱是原信号的相位谱再附加一个线性相移，由时移特性可以看到，信号的相位谱可以反映信号在时域中的位置信息，不同位置上的同一信号，它们具有不同的相频特性，而幅频特性相同。 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.10 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 = freqz(num2, 1, w); subplot(2,2,1) plot(w/pi,abs(h1));grid title('Magnitude Spectrum of Original Sequence','FontSize',8) xlabel('\omega /\pi'); ylabel('Amplitude'); subplot(2,2,2) plot(w/pi,abs(h2));grid title('Magnitude Spectrum of Frequency-Shifted Sequence','FontSize',8) xlabel('\omega /\pi'); ylabel('Amplitude'); subplot(2,2,3) plot(w/pi,angle(h1));grid title('Phase Spectrum of Original Sequence','FontSize',8) xlabel('\omega /\pi'); ylabel('Phase in radians'); subplot(2,2,4) plot(w/pi,angle(h2));grid title('Phase Spectrum of Frequency-Shifted Sequence','FontSize',8) xlabel('\omega /\pi'); ylabel('Phase in radians'); Wo控制平移量。 Q3.11由结果图Q3.11可得出在参数wo的控制下，离散时间傅里叶变换的幅度谱和相位谱都随着控制参数右移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.14 clf; w = -pi:2*pi/255:pi; % freqency vector for evaluating DTFT 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,w); subplot(2,2,1) plot(w/pi,abs(hp));grid title('Product of Magnitude Spectra','FontSize',8) xlabel('\omega /\pi'); ylabel('Amplitude'); subplot(2,2,2) plot(w/pi,abs(h3));grid title('Magnitude Spectrum of Convolved Sequence','FontSize',8) xlabel('\omega /\pi'); ylabel('Amplitude'); subplot(2,2,3) plot(w/pi,angle(hp));grid title('Sum of Phase Spectra','FontSize',8) xlabel('\omega /\pi'); ylabel('Phase in radians'); subplot(2,2,4) plot(w/pi,angle(h3));grid title('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.17 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); subplot(3,1,1) plot(w/pi,abs(h1));grid title('Magnitude Spectrum of First Sequence') xlabel('\omega /\pi'); ylabel('Amplitude'); subplot(3,1,2) plot(w/pi,abs(h2));grid title('Magnitude Spectrum of Second Sequence') xlabel('\omega /\pi'); ylabel('Amplitude'); subplot(3,1,3) plot(w/pi,abs(h3));grid title('Magnitude Spectrum of Product Sequence') xlabel('\omega /\pi'); ylabel('Amplitude'); Q3.18分析图得出乘积序列的幅度谱近似等于两序列的幅度谱的和. 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.20 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));grid title('Magnitude Spectrum of Original Sequence','FontSize',8) xlabel('\omega /\pi'); ylabel('Amplitude'); subplot(2,2,2) plot(w/pi,abs(h3));grid title('Magnitude Spectrum of Time-Reversed Sequence','FontSize',8) xlabel('\omega /\pi'); ylabel('Amplitude'); subplot(2,2,3) plot(w/pi,angle(h1));grid title('Phase Spectrum of Original Sequence','FontSize',8) xlabel('\omega /\pi'); ylabel('Phase in radians'); subplot(2,2,4) plot(w/pi,angle(h3));grid title('Phase Spectrum of Time-Reversed Sequence','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.23 clf; N=200; L=256; nn = [0:N-1]; kk = [0:L-1]; xR = [0.1*(1:100) 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('Re\{x[n]\}'); xlabel('Time index n'); ylabel('Amplitude'); subplot(3,2,2); plot(nn,xI);grid; title('Im\{x[n]\}'); xlabel('Time index n'); ylabel('Amplitude'); subplot(3,2,3); plot(kk,real(XF));grid; title('Re\{X[k]\}'); xlabel('Frequency index k'); ylabel('Amplitude'); subplot(3,2,4); plot(kk,imag(XF));grid; title('Im\{X[k]\}'); xlabel('Frequency index k'); ylabel('Amplitude'); xx = ifft(XF,L); subplot(3,2,5); plot(kk,real(xx));grid; title('Real part of IDFT\{X[k]\}'); xlabel('Time index n'); ylabel('Amplitude'); subplot(3,2,6); plot(kk,imag(xx));grid; title('Imag part of IDFT\{X[k]\}'); xlabel('Time index n'); ylabel('Amplitude'); Q3.24 clf; N=256; nn = [0:N-1]; ntime = [-N/2:N/2-1]; g = (0.75).^abs(ntime); h = (-0.9).^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'); xlabel('Frequency index k'); ylabel('Re\{G[k]\}'); subplot(2,2,2); plot(nn,imag(GF));grid; title('Two N-point DFT''s'); xlabel('Frequency index k'); ylabel('Im\{G[k]\}'); subplot(2,2,3);grid; plot(nn,real(GF2));grid; title('Single N-point DFT'); xlabel('Frequency index k'); ylabel('Re\{G[k]\}'); subplot(2,2,4); plot(nn,imag(GF2));grid; title('Single N-point DFT'); xlabel('Frequency index k'); ylabel('Im\{G[k]\}'); figure(2);clf; subplot(2,2,1);grid; plot(nn,real(HF));grid; title('Two N-point DFT''s'); xlabel('Freq index k'); ylabel('Re\{H[k]\}'); subplot(2,2,2); plot(nn,imag(HF));grid; title('Two N-point DFT''s'); xlabel('Freq index k'); ylabel('Im\{H[k]\}'); subplot(2,2,3);grid; plot(nn,real(HF2));grid; title('Single N-point DFT'); xlabel('Freq index k'); ylabel('Re\{H[k]\}'); subplot(2,2,4); plot(nn,imag(HF2));grid; title('Single N-point DFT'); xlabel('Freq index k'); ylabel('Im\{H[k]\}'); Q3.25 clf; 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 = [XFstar(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('Re\{V[k]\}'); subplot(2,2,2); plot(k,imag(VF));grid; title('Complex N-point DFT'); xlabel('Frequency index k'); ylabel('Im\{V[k]\}'); subplot(2,2,3); plot(k,real(VF2));grid; title('Real 2N-point DFT'); xlabel('Frequency index k'); ylabel('Re\{V[k]\}'); subplot(2,2,4); plot(k,imag(VF2));grid; title('Real 2N-point DFT'); xlabel('Frequency index k'); ylabel('Im\{V[k]\}'); 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输入是平等的两个向量x1和x2长度l .理解circonv是如何工作的,它是有用的定期x2的延伸。让x2p x2的无限长的周期延长。从概念上讲,常规时间逆转x2p和集x2tr 1到L等于元素的时间逆转x2p版本。元素1通过y L的输出向量然后通过x1和长度之间的内积向量sh循环变化对时间逆转向量x2tr。对于输出样例y[n],1≤n≤L、正确的循环移位是n - 1的位置。 Q3.30 clf; 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,min(a),max(a)]); title('Original Sequence'); xlabel('time index n'); ylabel('a[n]'); 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('b[n]') M值决定时移量。 Q3.31 Q3.32 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));grid; title('Magnitude of DFT of Original Sequence'); xlabel('Frequency index k'); ylabel('|X[k]|'); subplot(2,2,2); stem(n,abs(YF));grid; title('Magnitude of DFT of Circularly Shifted Sequence'); xlabel('Frequency index k'); ylabel('|Y[k]|'); subplot(2,2,3); stem(n,angle(XF));grid; title('Phase of DFT of Original Sequence'); xlabel('Frequency index k'); ylabel('arg(X[k])'); subplot(2,2,4); stem(n,angle(YF));grid; title('Phase of DFT of Circularly Shifted Sequence'); xlabel('Frequency index k'); ylabel('arg(Y[k])'); 时移量是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 = ');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.38 g1 = [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 circular 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.40 g1 = [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.41 x = [1 2 4 2 6 32 6 4 2 zeros(1,247)]; 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('Amplitude'); title('Re(DFT\{x[n]\})'); subplot(2,2,2); plot(k/128,imag(XF)); grid; ylabel('Amplitude'); title('Im(DFT\{x[n]\})'); subplot(2,2,3); plot(k/128,real(XEF)); grid; xlabel('Time index n');ylabel('Amplitude'); title('Re(DFT\{x_{e}[n]\})'); subplot(2,2,4); plot(k/128,imag(XEF)); grid; xlabel('Time index n');ylabel('Amplitude'); title('Im(DFT\{x_{e}[n]\})'); X1[n]是X[n]的逆转序列 Q3.42 XEF等于零的虚部在浮点精度。这个结果可以解释如下:真正的转换的一部分x[n]的转换定期甚至x[n]的一部分。因此,定期的DFT甚至x[n]的一部分，真正的一部分,正是真正的X[k]的一部分,一个虚部为零。 Q3.43仿真图为上右图。 Q3.44a和b的值相等。 x = [(1:128) (128:-1:1)]; XF = fft(x); a = sum(x.*x) b = round(sum(abs(XF).^2)/256) Q3.45 x = [(1:128) (128:-1:1)]; XF = fft(x); a = sum(x.*x) b = round(sum(XF.*conj(XF))/256) Q3.46 Q3.47运行结果为 Zeros: -1.0000 + 1.4142i -1.0000 - 1.4142i -0.2500 + 0.6614i -0.2500 - 0.6614i Poles: -8.9576 -0.2718 0.1147 + 0.2627i 0.1147 - 0.2627i sos = 1.0000 2.0000 3.0000 1.0000 9.2293 2.4344 1.0000 0.5000 0.5000 1.0000 -0.2293 0.0822 k = 0.4000 Q3.48  结果为： Q3.49 clf; z = [0.3 2.5 -0.2+i*0.4 -0.2-i*0.4]'; p = [0.5 -0.75 0.6+i*0.7 0.6-i*0.7]'; k = 3.9; [num den] = zp2tf(z,p,k) num = 3.9000 -9.3600 -0.6630 -1.0140 0.5850 den = 1.0000 -0.9500 0.1750 0.6625 -0.3187 Q3.50    L=50 Q3.51 clf; num = [2 5 9 5 3]; den = [5 45 2 1 1]; [r p k] = residuez(num,den) r = 0.3109 -1.0254 - 0.3547i -1.0254 + 0.3547i -0.8601 p = -8.9576 0.1147 + 0.2627i 0.1147 - 0.2627i -0.2718 k =