数字信号处理第三章.doc

数字信号处理第三章实验程序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);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(Amplitude);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;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 /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(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 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 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);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(Amplitude);subplot(2,1,2)plot(w/pi,angle(h);gridtitle(Phase Spectrum 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(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);gridtitle(Magnitude Spectrum of Original Sequence,FontSize,8)xlabel(omega /pi);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,FontSize,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 时移特性：信号在时域移动某个距离，则所得信号的幅度谱和原信号相同，而相位谱是原信号的相位谱再附加一个线性相移，由时移特性可以看到，信号的相位谱可以反映信号在时域中的位置信息，不同位置上的同一信号，它们具有不同的相频特性，而幅频特性相同。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 = 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);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(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的控制下，离散时间傅里叶变换的幅度谱和相位谱都随着控制参数右移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);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);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 Convolved 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);gridtitle(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 = -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(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分析图得出乘积序列的幅度谱近似等于两序列的幅度谱的和. 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,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);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,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) 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);subplot(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(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).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 DFTs);xlabel(Frequency index k);ylabel(ReGk);subplot(2,2,2);plot(nn,imag(GF);grid;title(Two N-point DFTs);xlabel(Frequency index k);ylabel(ImGk);subplot(2,2,3);grid;plot(nn,real(GF2);grid;title(Single N-point DFT);xlabel(Frequency 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 DFTs);xlabel(Freq index k);ylabel(ReHk);subplot(2,2,2);plot(nn,imag(HF);grid;title(Two N-point DFTs);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);subplot(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 = 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(ReVk);subplot(2,2,2);plot(k,imag(VF);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);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输入是平等的两个向量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,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值决定时移量。Q3.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(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);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 = );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 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.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);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);plo