数字信号处理大作业.docx

1、数字信号处理大作业 数字信号处理大作业 电子工程学院 M2.2 The square wave and the sawtooth wave are two periodic sequence as sketched in Figure P2.1. Using the functions sawtooth and square write a MATLAB program to generate the above two sequences and plot them and the period N. For the square wave sequence an additional us

2、er-specified parameter is the duty cycle,while is the percent of the period for which the signal is positive. Using this program generate the first 100 samples of each of the above sequences with a sampling rate of 20 khz,a peak value of 7, a period of 13 and a duty cycle of 60% for the square wave.

3、如下为matlab程序代码： %get user inputsa=input(the peak value=);l=input(length of sequence=);n=input (the period of sequence);ft=input(the desired sampling frequency=);dc=input(the square wave duty cycle =);%creat singalst=1/ft;t=0:l-1;x=a*sawtooth(2*pi*t/n);y=a*square(2*pi*(t/n),dc);%polt subplot(211)stem(

4、t,x)ylabel(Amplitude);xlabel(Time in ,num2str(t),sec);subplot(212)stem(t,y);ylabel(Amplitude);xlabel(Time in ,num2str(t),sec);在仿真之前，分析题目中所给出的数据，需要注意此时的周期为序列的个数13，而非时间周期13s，否则100个取样最终时间才花费0.05s，在图中就无法显示出来。同时占空比很明确的影响方波的取样正负所占比例，这也与占空比所本意相吻合。周期的多少也直接影响一个周期内取样的个数。M 2.4 (a) Write a Matlab program to gen

5、erate a sinusoidal sequence xn=Acos(0 n+）and plot the sequence using the stem function. The input data specified by the user are the desired length L,amplitude A,the angular frequency 0,and the phase where 00 L = input(Desired length = );A = input(Amplitude = );omega = input(Angular frequency = );ph

6、i = input(Phase = );n = 0:L-1;x = A*cos(omega*n + phi);stem(n,x);xlabel(Time index);ylabel(Amplitude);title(omega_o = ,num2str(omega);%example：Desired length = 100Amplitude = 1.5Angular frequency = 0.1*piPhase = 00=0.10=0.14由理论计算可知其周期t=2pi/0.14pi=100/7，故时间周期T=100，由图中也可明显看出，其波形以100开始循坏，故其周期也为100，相符。.

7、0=0.24理论计算t=2pi/0.24pi=25/3，故其周期为25，由图中也可看出，其图像每25一循坏，周期为25，与理论也相符。w0=0.34pi理论计算t=2pi/0.34pi=100/17，故其周期为100，由图中，基本上也看不出来周期，但是可以通过理论计算值来从图中看出隔100是大致高度图像相同的0=0.68这次可以明显看出，周期应该是在50附近，理论计算t=2pi/0.68pi=50/17，周期T=50，符合。0=0.75明显看出，在第八个取样点出出现了重复现象，因此其周期为8，而通过理论计算的出的结果t=2pi/0.5pi=8/3，周期T=8，一致。有这么多的取样可以看出，取样

8、频率的不同，会大大的影响得到的采样点的图形。部分图像中会出现类似于几个正弦波叠加的情形，我的理解应该是因为周期太大，以至于一个周期内的取样点很多，压缩在一起就近似的恢复道了原来的图形，就得到了类似几个正弦序列叠加的情形，而且随着周期的进一步增大，其叠加的近似数目也会增多。M2.5 Generate the sequences or problem 2.21(b) to 2.21(b) using matlab.代码同上M2.4：b.xn=sin(0.6*pi*n+0.6pi)c.d.e.这个题目是熟悉matlab应用，熟悉一下它的函数多输入功能，对于多个函数可以多几个函数，加几个变量，如果频率

9、或者振幅一致，还可以取巧让频率振幅参量表示为相同字母，减少变量的输入.M2.6 Write a matlab program to plot a continuous-time sinusoidal signal and its sampled version and verify figure 2.19.You need to use the hold function to keep both plots.Answer：t = 0:0.001:1;fo = input(Frequency of sinusoid in Hz = );FT = input(Sample frequency i

10、n Hz = );g1 = cos(2*pi*fo*t);plot(t,g1,-)xlabel(time);ylabel(Amplitude)holdn = 0:1:FT;gs = cos(2*pi*fo*n/FT);title(g1(t)=cos(0.6*pi*t),实线表示，抽样序列表示为圆：);plot(n/FT,gs,o);hold off取样后，如图所示，频率为3,7,13Hz的函数采样结果都相同，都为cos（0.6*pi*t），其实所有频率为10K3的余弦波形，采样的结果都一样。这种混叠现象，会让我们在进行回复原来的x（t）时出现不确定性，因此也说明，采样和恢复不是一一对应的关系。

11、要恢复原来的确定的信号，也就必须要在附加一些很必要的条件才行。M3.1 Using program 3.1 determine and plot the real and imaginary parts and the magnitude and phase spectra of the following DTFT for various values of r and ：(a) 程序及结果如下： %读入离散时间傅里叶变换所需的长度k=input(频率点数量=);%读入分子和分母系数num=input(分子系数=);den=input(分母系数);%计算频率响应w=0:pi/(k-1):pi

12、;h=freqz(num,den,w);%绘制频率响应subplot(2,2,1)plot(w/pi,real(h);gridtitle(实部)xlabel(omega/pi);ylabel(振幅)subplot(2,2,2)plot(w/pi,imag(h);gridtitle(虚部)xlabel(omega/pi);ylabel(振幅)subplot(2,2,3)plot(w/pi,abs(h);gridtitle(幅度谱)xlabel(omega/pi);ylabel(幅度)subplot(2,2,4)plot(w/pi,angle(h);gridtitle(相位谱)xlabel(ome

13、ga/pi);ylabel(相位，弧度)取r=0.5，=60，（b）取r=0.8，=90，如下：又这两幅图中可以看出，幅度谱的变化基本和实部的变化一致，也是由于虚部的变化转折点基本和实部一致的原因，所以通过这两幅图也可看出，实部虚部与幅度的关系。，而通过相位谱和幅度谱的关系，可以看出，幅度谱的转折点也几乎就是相位谱的转折点，同时也与虚部的变化几乎一致，它们之间的联系是非常紧密的，因此分析一个复数的时候，就应该把这四个方面都考虑进去，以便于更好的理解和学习。M3.4 Using matlab verify the following general properties of the DTFT

14、as listed in Table 3.2:(a)Linearity(b)time-shifting,(c)frequency-shifting(d).differentiation-in-frequency.(e)convolution.(f).modulation.and (g)Parsevals relation. Since all data in matlab have to be finite-length vectors, the sequences to be used to verify the properties are thus restricted to be of

15、 finite length.Answer:(a)线性：代码如下N=input(The length of the sequence=);k=0:N-1;gamma=-0.5;g=exp(gamma*k);%g is an exponential sequenceh=sin(2*pi*k/(N/2);%h is a sinusoidal sequence with period =N/2G,W=freqz(g,1,512);H,W=freqz(h,1,512);%property 1alpha=0.5;beta=0.25;y=alpha*g+beta*h;Y,W=freqz(y,1,512);

16、比较Y和alpha*G+beta*H即可：可以看出两者完全一样，说明线性是完全成立的。(b)时移代码如下：%property 2n0=5;% s sequence shifted by 5 simplesy2=zeros(1,n0) g;Y2,W=freqz(y2,1,512);G0=exp(-j*W*n0).*G;比较G0和deY2即可：可以看出，两种方法得出来的频谱一样，说明时移性质正确。（c）频移代码如下:%property 3w0=pi/2;% the value of cmega0=pi/2r=256;%the value of omega0 in terms of number o

17、f samplesy3=g.*exp(j*w0*k);Y3,w=freqz(y3,1,512);k=0:511;w=-w0+pi*k/512;% creating G(exp(w-w0)G1=freqz(g,1,w);比较G1和Y3即可：由此可见，频移也满足。（d）频域微分代码如下：%property 4k=0:N-1;y4=k.*g;Y4,w=freqz(y4,1,512);y0=(-1).k).*g;G2=G(2:512) sum(y0);delG=(G2-G)*512/pi;比较Y4和delG即可：可见，频域微分性质也成立。(e) 卷积代码如下：%property 5y5=conv(g,

18、h);Y5,w=freqz(y5,1,512);比较Y5和G.*H即可：（f） 调制%property 6y6=g.*h;Y6,w=freqz(y6,1,512,whole);G0,w=freqz(g,1,512,whole);H0,w=freqz(h,1,512,whole); H1=fliplr(H0(1:129) fliplr(H0(130:512); val=1/(512)*sum(G0.*H1);比较Y6和val即可：调制也完全符合。（g） Parsevals relation%property 7val1=sum(g.*conj(h);val2=sum(G0.*conj(H0)/5

19、12;比较val1和val2可得如下：由此可见，Parsevals relation也完全符合。 这个题目，也通过仿真验证了傅里叶变换的这七个性质，也从实验验证了这几个性质的正确性。通过6调制也可以很便利的计算个别点的调制特性。这也便于我们更好的来理解傅里叶变换的这几个性质，也更利于我们深入的理解和学习。其中，性质4没有理解程序意义，没有明白其原因结构，因此结果可能只是样子。但是这个并不影响性质的正确与否。M3.8 Using MATLAB compute the N-point DFTs of the length -N sequence of Problem 3.12 for N=3,5,

20、7,and 10.Compare your results with that obtained by evaluating the DTFTs computed in Problem 3.12 at =2k/N. k=0,1,2,3,N-1.(a).y1n=1, -NnN，其他为0.代码如下：N = input(The value of N = );k = -N:N;y1 = ones(1,2*N+1);w = 0:2*pi/255:2*pi;Y1 = freqz(y1, 1, w);Y1dft = fft(y1);k = 0:1:2*N;plot(w/pi,abs(Y1),k*2/(2*N

21、+1),abs(Y1dft),o);xlabel(Normalized frequency);ylabel(Amplitude);title(3点DFT);(b).y2n=1-n/N,-NnN，其他为0.(c).y3n=cos(pi*n/2N), -NnN，其他为0. 本道题中列举了3种序列的不同点的离散傅里叶变换的变换后的图像。可以明确的看出，傅里叶变换的长度越大，其涵盖的内容就相应的越多，当然也就越来越精确，也就更符合原来函数的性质。而且相应的，长度约长，其变换得到的图形高度两边的越来越高，中间的点也就越低，可能是由于点多的缘故，分散了原来所等量的能量，因此，点变多了，其幅度相应的也就要有

22、所改变。同时也可以看出，离散傅里叶变换具有对称特性，因此计算傅里叶的时候，可以计算一半的系数，另一半对照类比即可。这其实也为后面的蝶形运算提供了依据。M3.19 Using Program 3.10 determine the z-transform as a ratio of two polynomials in z-1 from each of the part fraction expansions listed below. (d).X4(z)=4+10/(5+2z-1 )+z-1 /(6+5z-1+z -2), z0.5,Answer:（a）:（b）程序代码如下：r=input(请输

23、入留数=);请输入留数=0 r=input(请输入留数=);p=input(请输入极点=);k=input(请输入常量=);num,den=residuez(r,p,k);disp(分子多项式系数);disp(num)disp(分母多项式系数);disp(den)由MATLAB得出的结果如下：则X1(z)=-(3.5+1.25*z-1+0.25*z-2)/(1+0.75*z-1+0.125*z-2).(b).先将后面一项化成分母一次的单项，待定系数法可以得到-(3+z-1)/（1-0.25z-2）=0.5/(1+0.5z-1)-2.5/(1-0.5z-1),再化简，顾可得如下：X2(z)=-(0.5+2.5z-1+0.875z-2)/(1-0.25z-2).(c)同理，将最后一项化成两单项。最后结果如下：(d)同上，可得结果： 通过此题，我们掌握了一种用matlab快速转换单式与复式的方法，如果将x（z）看做一个系统的传递函数，这就为我们提供了计算极点和零点的快捷方法，也便于我们来将整合后的复式来分析系统的快速性和稳定性。