数字信号处理第三章

数字信号处理第三章数字信号处理第三章实验程序3.1计算离散时间傅里叶变换%ProgramP3_1%EvaluationoftheDTFTclf;%ComputethefrequencysamplesoftheDTFTw=-4*pi:8*pi/511:4*pi;num=[21];den=[1-0.6];h=freqz(num,den,w);%PlottheDTFTsubplot(2,1,1)plot(w/pi,real(h));gridtitle('RealpartofH(e^{j\omega})')xlabel('\omega/\pi');ylabel('Amplitude');subplot(2,1,2)plot(w/pi,imag(h));gridtitle('ImaginarypartofH(e^{j\omega})')xlabel('\omega/\pi');ylabel('Amplitude');pausesubplot(2,1,1)plot(w/pi,abs(h));gridtitle('MagnitudeSpectrum|H(e^{j\omega})|')xlabel('\omega/\pi');ylabel('Amplitude');subplot(2,1,2)plot(w/pi,angle(h));gridtitle('PhaseSpectrumarg[H(e^{j\omega})]')xlabel('\omega/\pi');ylabel('Phaseinradians');Q3.1离散时间傅里叶变换的原始序列是H(e^jw)=(2+z^-1)/(1-0.6z^-1)。Pause的作用是暂停等待用户输入任意键后接着执行以下命令。Q3.2是周期函数，周期是2π。实部和幅度谱是关于y轴对称，是偶函数；虚部和相位谱是关于原点对称，是奇函数。Q3.3clf;N=512;num=[0.7-0.50.31];den=[10.3-0.50.7];[h,w]=freqz(num,den,N);subplot(2,1,1)plot(w/pi,real(h));gridtitle('RealpartofH(e^{j\omega})')xlabel('\omega/\pi');ylabel('Amplitude');subplot(2,1,2)plot(w/pi,imag(h));gridtitle('ImaginarypartofH(e^{j\omega})')xlabel('\omega/\pi');ylabel('Amplitude');pausesubplot(2,1,1)plot(w/pi,abs(h));gridtitle('MagnitudeSpectrum|H(e^{j\omega})|')xlabel('\omega/\pi');ylabel('Amplitude');subplot(2,1,2)plot(w/pi,angle(h));gridtitle('PhaseSpectrumarg[H(e^{j\omega})]')xlabel('\omega/\pi');ylabel('Phaseinradians');还是周期函数，周期是2π。相位谱的跳变的原因是：在利用反正切函数计算角度的时候，其中的一个分支出现了衰减，造成了跳变。clf;N=512;num=[0.7-0.50.31];den=[10.3-0.50.7];[h,w]=freqz(num,den,N);subplot(2,1,1)plot(w/pi,unwrap(angle(h)));gridtitle('PhaseSpectrumarg[H(e^{j\omega})]')xlabel('\omega/\pi');ylabel('Phaseinradians');Q3.4修改后的程序为clf;w=-4*pi:8*pi/511:4*pi;num=[1357911131517];den=1;h=freqz(num,den,w);%PlottheDTFTsubplot(2,1,1)plot(w/pi,real(h));gridtitle('RealpartofH(e^{j\omega})')xlabel('\omega/\pi');ylabel('Amplitude');subplot(2,1,2)plot(w/pi,imag(h));gridtitle('ImaginarypartofH(e^{j\omega})')xlabel('\omega/\pi');ylabel('Amplitude');pausesubplot(2,1,1)plot(w/pi,abs(h));gridtitle('MagnitudeSpectrum|H(e^{j\omega})|')xlabel('\omega/\pi');ylabel('Amplitude');subplot(2,1,2)plot(w/pi,angle(h));gridtitle('PhaseSpectrumarg[H(e^{j\omega})]')xlabel('\omega/\pi');ylabel('Phaseinradians');是周期函数，周期是2π。实部和幅度谱是关于y轴对称，是偶函数；虚部和相位谱是关于原点对称，是奇函数。Q3.5若要改为以度为单位，则将程序中的第二个图的程序改为subplot(2,1,2)plot(w/pi,180*angle(h)/pi);gridtitle('PhaseSpectrumarg[H(e^{j\omega})]')xlabel('\omega/\pi');ylabel('Phaseindegrees');就可以了。3.2离散时间傅里叶变换的性质时移特性clf;w=-pi:2*pi/255:pi;D=10;num=[123456789];h1=freqz(num,1,w);h2=freqz([zeros(1,D)num],1,w);subplot(2,2,1)plot(w/pi,abs(h1));gridtitle('MagnitudeSpectrumofOriginalSequence','FontSize',8)xlabel('\omega/\pi');ylabel('Amplitude');subplot(2,2,2)plot(w/pi,abs(h2));gridtitle('MagnitudeSpectrumofTime-ShiftedSequence','FontSize',8)xlabel('\omega/\pi');ylabel('Amplitude');subplot(2,2,3)plot(w/pi,angle(h1));gridtitle('PhaseSpectrumofOriginalSequence','FontSize',8)xlabel('\omega/\pi');ylabel('Phaseinradians');subplot(2,2,4)plot(w/pi,angle(h2));gridtitle('PhaseSpectrumofTime-ShiftedSequence','FontSize',8)xlabel('\omega/\pi');ylabel('Phaseinradians');Q3.6参数D控制时移量。Q3.7D=10D=50时移特性：信号在时域移动某个距离，则所得信号的幅度谱和原信号相同，而相位谱是原信号的相位谱再附加一个线性相移，由时移特性可以看到，信号的相位谱可以反映信号在时域中的位置信息，不同位置上的同一信号，它们具有不同的相频特性，而幅频特性相同。Q3.8如上图所示Q3.9改变序列长度num=[1234567891011121314151617181920212223242526272829];所得的图像为D=10D=50从上图中可以看出，增加序列的长度，使得幅度谱更加窄，而相位谱则更加密集和陡峭。平移特性Q3.10clf;w=-pi:2*pi/255:pi;wo=0.4*pi;num1=[1357911131517];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('MagnitudeSpectrumofOriginalSequence','FontSize',8)xlabel('\omega/\pi');ylabel('Amplitude');subplot(2,2,2)plot(w/pi,abs(h2));gridtitle('MagnitudeSpectrumofFrequency-ShiftedSequence','FontSize',8)xlabel('\omega/\pi');ylabel('Amplitude');subplot(2,2,3)plot(w/pi,angle(h1));gridtitle('PhaseSpectrumofOriginalSequence','FontSize',8)xlabel('\omega/\pi');ylabel('Phaseinradians');subplot(2,2,4)plot(w/pi,angle(h2));gridtitle('PhaseSpectrumofFrequency-ShiftedSequence','FontSize',8)xlabel('\omega/\pi');ylabel('Phaseinradians');Wo控制平移量。Q3.11由结果图Q3.11可得出在参数wo的控制下，离散时间傅里叶变换的幅度谱和相位谱都随着控制参数右移k个单位（wo=k*pi）。k=0.4k=-0.4Q3.12将k改为-0.4得到的运行结果如上图。Q3.13改变序列长度序列：num1=[1357911131517192123252729]序列：num2=[111315171921232527293133353739];卷积性质Q3.14clf;w=-pi:2*pi/255:pi;%freqencyvectorforevaluatingDTFTx1=[1357911131517];x2=[1-23-21];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('ProductofMagnitudeSpectra','FontSize',8)xlabel('\omega/\pi');ylabel('Amplitude');subplot(2,2,2)plot(w/pi,abs(h3));gridtitle('MagnitudeSpectrumofConvolvedSequence','FontSize',8)xlabel('\omega/\pi');ylabel('Amplitude');subplot(2,2,3)plot(w/pi,angle(hp));gridtitle('SumofPhaseSpectra','FontSize',8)xlabel('\omega/\pi');ylabel('Phaseinradians');subplot(2,2,4)plot(w/pi,angle(h3));gridtitle('PhaseSpectrumofConvolvedSequence','FontSize',8)xlabel('\omega/\pi');ylabel('Phaseinradians');Q3.15分析结果图可以得出幅度谱的乘积和卷积后的幅度谱相同，相位谱的乘积和卷积后的相位谱相同。Q3.16x1=[13579111315171921232527293133];x2=[1-23-21-52-31];运行结果如上边第二个图所示。调制性质Q3.17clf;w=-pi:2*pi/255:pi;x1=[1357911131517];x2=[1-11-11-11-11];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('MagnitudeSpectrumofFirstSequence')xlabel('\omega/\pi');ylabel('Amplitude');subplot(3,1,2)plot(w/pi,abs(h2));gridtitle('MagnitudeSpectrumofSecondSequence')xlabel('\omega/\pi');ylabel('Amplitude');subplot(3,1,3)plot(w/pi,abs(h3));gridtitle('MagnitudeSpectrumofProductSequence')xlabel('\omega/\pi');ylabel('Amplitude');Q3.18分析图得出乘积序列的幅度谱近似等于两序列的幅度谱的和.Q3.19将序列改变为x1=[13579111315171921232527];x2=[1-11-11-11-1102-47-1]得到的运行结果为上右图。乘积序列的幅度谱依然近似等于两序列的幅度谱的和.时间反转性质Q3.20clf;w=-pi:2*pi/255:pi;num=[1234];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('MagnitudeSpectrumofOriginalSequence','FontSize',8)xlabel('\omega/\pi');ylabel('Amplitude');subplot(2,2,2)plot(w/pi,abs(h3));gridtitle('MagnitudeSpectrumofTime-ReversedSequence','FontSize',8)xlabel('\omega/\pi');ylabel('Amplitude');subplot(2,2,3)plot(w/pi,angle(h1));gridtitle('PhaseSpectrumofOriginalSequence','FontSize',8)xlabel('\omega/\pi');ylabel('Phaseinradians');subplot(2,2,4)plot(w/pi,angle(h3));gridtitle('PhaseSpectrumofTime-ReversedSequence','FontSize',8)xlabel('\omega/\pi');ylabel('Phaseinradians');Q3.21分析图得出序列的幅度谱随时间反转不发生变化，序列相位谱随时间反转而反转180。Q3.22改变序列长度num=[1-23-45-67-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('Re\{x[n]\}');xlabel('Timeindexn');ylabel('Amplitude');subplot(3,2,2);plot(nn,xI);grid;title('Im\{x[n]\}');xlabel('Timeindexn');ylabel('Amplitude');subplot(3,2,3);plot(kk,real(XF));grid;title('Re\{X[k]\}');xlabel('Frequencyindexk');ylabel('Amplitude');subplot(3,2,4);plot(kk,imag(XF));grid;title('Im\{X[k]\}');xlabel('Frequencyindexk');ylabel('Amplitude');xx=ifft(XF,L);subplot(3,2,5);plot(kk,real(xx));grid;title('RealpartofIDFT\{X[k]\}');xlabel('Timeindexn');ylabel('Amplitude');subplot(3,2,6);plot(kk,imag(xx));grid;title('ImagpartofIDFT\{X[k]\}');xlabel('Timeindexn');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('TwoN-pointDFT''s');xlabel('Frequencyindexk');ylabel('Re\{G[k]\}');subplot(2,2,2);plot(nn,imag(GF));grid;title('TwoN-pointDFT''s');xlabel('Frequencyindexk');ylabel('Im\{G[k]\}');subplot(2,2,3);grid;plot(nn,real(GF2));grid;title('SingleN-pointDFT');xlabel('Frequencyindexk');ylabel('Re\{G[k]\}');subplot(2,2,4);plot(nn,imag(GF2));grid;title('SingleN-pointDFT');xlabel('Frequencyindexk');ylabel('Im\{G[k]\}');figure(2);clf;subplot(2,2,1);grid;plot(nn,real(HF));grid;title('TwoN-pointDFT''s');xlabel('Freqindexk');ylabel('Re\{H[k]\}');subplot(2,2,2);plot(nn,imag(HF));grid;title('TwoN-pointDFT''s');xlabel('Freqindexk');ylabel('Im\{H[k]\}');subplot(2,2,3);grid;plot(nn,real(HF2));grid;title('SingleN-pointDFT');xlabel('Freqindexk');ylabel('Re\{H[k]\}');subplot(2,2,4);plot(nn,imag(HF2));grid;title('SingleN-pointDFT');xlabel('Freqindexk');ylabel('Im\{H[k]\}');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=[GFGF]+(W2N.^k).*[HFHF];VF2=fft(v);abs(max(VF-VF2))subplot(2,2,1);plot(k,real(VF));grid;title('ComplexN-pointDFT');xlabel('Frequencyindexk');ylabel('Re\{V[k]\}');subplot(2,2,2);plot(k,imag(VF));grid;title('ComplexN-pointDFT');xlabel('Frequencyindexk');ylabel('Im\{V[k]\}');subplot(2,2,3);plot(k,real(VF2));grid;title('Real2N-pointDFT');xlabel('Frequencyindexk');ylabel('Re\{V[k]\}');subplot(2,2,4);plot(k,imag(VF2));grid;title('Real2N-pointDFT');xlabel('Frequencyindexk');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的延伸。让x2px2的无限长的周期延长。从概念上讲,常规时间逆转x2p和集x2tr1到L等于元素的时间逆转x2p版本。元素1通过yL的输出向量然后通过x1和长度之间的内积向量sh循环变化对时间逆转向量x2tr。对于输出样例y[n],1≤n≤L、正确的循环移位是n-1的位置。Q3.30clf;M=6;a=[0123456789];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('OriginalSequence');xlabel('timeindexn');ylabel('a[n]');subplot(2,1,2);stem(n,b);axis([0,L,min(a),max(a)]);title(['SequenceObtainedbyCircularlyShiftingby',num2str(M),'Samples']);xlabel('timeindexn');ylabel('b[n]')M值决定时移量。Q3.31Q3.32clf;x=[0246810121416];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('MagnitudeofDFTofOriginalSequence');xlabel('Frequencyindexk');ylabel('|X[k]|');subplot(2,2,2);stem(n,abs(YF));grid;title('MagnitudeofDFTofCircularlyShiftedSequence');xlabel('Frequencyindexk');ylabel('|Y[k]|');subplot(2,2,3);stem(n,angle(XF));grid;title('PhaseofDFTofOriginalSequence');xlabel('Frequencyindexk');ylabel('arg(X[k])');subplot(2,2,4);stem(n,angle(YF));grid;title('PhaseofDFTofCircularlyShiftedSequence');xlabel('Frequencyindexk');ylabel('arg(Y[k])');时移量是8.Q3.33Q3.34M=5运行结果如上右图所示。Q3.35Length=13Length=20Q3.36g1=[123456];g2=[1-233-21];ycir=circonv(g1,g2);disp('Resultofcircularconvolution=');disp(ycir)G1=fft(g1);G2=fft(g2);yc=real(ifft(G1.*G2));disp('ResultofIDFToftheDFTproducts=');disp(yc)运行结果Q3.37结果如下：Q3.38g1=[12345];g2=[22011];g1e=[g1zeros(1,length(g2)-1)];g2e=[g2zeros(1,length(g1)-1)];ylin=circonv(g1e,g2e);disp('Linearconvolutionviacircularconvolution=');disp(ylin);y=conv(g1,g2);disp('Directlinearconvolution=');disp(y)结果如下：Q3.39g1=[3141592];g2=[11100];g1=[543210];g2=[-21234];Q3.40g1=[12345];g2=[22011];g1e=[g1zeros(1,length(g2)-1)];g2e=[g2zeros(1,length(g1)-1)];G1EF=fft(g1e);G2EF=fft(g2e);ylin=real(ifft(G1EF.*G2EF));disp('直线线性卷积=');disp(ylin)；Q3.41x=[1242632642zeros(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('Timeindexn');ylabel('Amplitude');title('Re(DFT\{x_{e}[n]\})');subplot(2,2,4);plot(k/128,imag(XEF));grid;xlabel('Timeindexn');ylabel('Amplitude');title('Im(DFT\{x_{e}[n]\})');X1[n]是X[n]的逆转序列Q3.42XEF等于零的虚部在浮点精度。这个结果可以解释如下:真正的转换的一部分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.45x=[(1:128)(128:-1:1)];XF=fft(x);a=sum(x.*x)b=round(sum(XF.*conj(XF))/256)Q3.46Q3.47运行结果为Zeros:-1.0000+1.4142i-1.0000-1.4142i-0.2500+0.6614i-0.2500-0.6614iPoles:-8.9576-0.27180.1147+0.2627i0.1147-0.2627isos=1.00002.00003.00001.00009.22932.43441.00000.50000.50001.0000-0.22930.0822k=0.4000Q3.48结果为：Q3.49clf;z=[0.32.5-0.2+i*0.4-0.2-i*0.4]';p=[0.5-0.750.6+i*0.70.6-i*0.7]';k=3.9;[numden]=zp2tf(z,p,k)num=3.9000-9.3600-0.6630-1.01400.5850den=1.0000-0.95000.17500.6625-0.3187Q3.50L=50Q3.51clf;num=[25953];den=[545211];[rpk]=residuez(num,den)r=0.3109-1.0254-0.3547i-1.0254+0.3547i-0.8601p=-8.95760.1147+0.2627i0.1147-0.2627i-0.2718k=