数字信号处理大作业.docx
《数字信号处理大作业.docx》由会员分享，可在线阅读，更多相关《数字信号处理大作业.docx（24页珍藏版）》请在冰豆网上搜索。
数字信号处理大作业
数字信号处理大作业
电子工程学院
M2.2ThesquarewaveandthesawtoothwavearetwoperiodicsequenceassketchedinFigureP2.1.UsingthefunctionssawtoothandsquarewriteaMATLABprogramtogeneratetheabovetwosequencesandplotthemandtheperiodN.Forthesquarewavesequenceanadditionaluser-specifiedparameteristhedutycycle,whileisthepercentoftheperiodforwhichthesignalispositive.Usingthisprogramgeneratethefirst100samplesofeachoftheabovesequenceswithasamplingrateof20khz,apeakvalueof7,aperiodof13andadutycycleof60%forthesquarewave.
如下为matlab程序代码：
%getuserinputs
a=input（'thepeakvalue='）;
l=input（'lengthofsequence='）;
n=input（'theperiodofsequence'）;
ft=input（'thedesiredsamplingfrequency='）;
dc=input（'thesquarewavedutycycle='）;
%creatsingals
t=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（t,x）
ylabel（'Amplitude'）;
xlabel（['Timein',num2str（t）,'sec']）;
subplot（212）
stem（t,y）;
ylabel（'Amplitude'）;
xlabel（['Timein',num2str（t）,'sec']）;
在仿真之前，分析题目中所给出的数据，需要注意此时的周期为序列的个数13，而非时间周期13s，否则100个取样最终时间才花费0.05s，在图中就无法显示出来。
同时占空比很明确的影响方波的取样正负所占比例，这也与占空比所本意相吻合。
周期的多少也直接影响一个周期内取样的个数。
M2.4（a）WriteaMatlabprogramtogenerateasinusoidalsequencex[n]=Acos（ω0n+φ）andplotthesequenceusingthestemfunction.TheinputdataspecifiedbytheuserarethedesiredlengthL,amplitudeA,theangularfrequencyω0,andthephaseφwhere0<ω0<πand0≤φ≤2π.Usingtheprogramgeneratethesinusoidalsequenceshowninfigure2.15.
（b）Generatesinusoidalsequenceswiththeangularfrequenciesgiveninproblem2.2.Determinetheperiodofeachsequencefromtheplotandverifytheresulttheoretically.
Answer:
如下为matlab代码：
>>L=input（'Desiredlength='）;
A=input（'Amplitude='）;
omega=input（'Angularfrequency='）;
phi=input（'Phase='）;
n=0:
L-1;
x=A*cos（omega*n+phi）;
stem（n,x）;
xlabel（'Timeindex'）;ylabel（'Amplitude'）;
title（['\omega_{o}=',num2str（omega）]）;
%%example：
Desiredlength=100
Amplitude=1.5
Angularfrequency=0.1*pi
Phase=0
ω0=0.1π
ω0=0.14π
由理论计算可知其周期t=2pi/0.14pi=100/7，故时间周期T=100，由图中也可明显看出，其波形以100开始循坏，故其周期也为100，相符。
.
ω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，一致。
有这么多的取样可以看出，取样频率的不同，会大大的影响得到的采样点的图形。
部分图像中会出现类似于几个正弦波叠加的情形，我的理解应该是因为周期太大，以至于一个周期内的取样点很多，压缩在一起就近似的恢复道了原来的图形，就得到了类似几个正弦序列叠加的情形，而且随着周期的进一步增大，其叠加的近似数目也会增多。
M2.5Generatethesequencesorproblem2.21（b）to2.21（b）usingmatlab.
代码同上M2.4：
b.x[n]=sin（0.6*pi*n+0.6pi）
c.
d.
e.
这个题目是熟悉matlab应用，熟悉一下它的函数多输入功能，对于多个函数可以多几个函数，加几个变量，如果频率或者振幅一致，还可以取巧让频率振幅参量表示为相同字母，减少变量的输入.
M2.6Writeamatlabprogramtoplotacontinuous-timesinusoidalsignalanditssampledversionandverifyfigure2.19.Youneedtousetheholdfunctiontokeepbothplots.
Answer：
t=0:
0.001:
1;
fo=input（'FrequencyofsinusoidinHz='）;
FT=input（'SamplefrequencyinHz='）;
g1=cos（2*pi*fo*t）;
plot（t,g1,'-'）
xlabel（'time'）;ylabel（'Amplitude'）
hold
n=0:
1:
FT;
gs=cos（2*pi*fo*n/FT）;
title（'g1（t）=cos（0.6*pi*t）,实线表示，抽样序列表示为圆：
'）;
plot（n/FT,gs,'o'）;holdoff
取样后，如图所示，频率为3,7,13Hz的函数采样结果都相同，都为cos（0.6*pi*t），其实所有频率为10K±3的余弦波形，采样的结果都一样。
这种混叠现象，会让我们在进行回复原来的x（t）时出现不确定性，因此也说明，采样和恢复不是一一对应的关系。
要恢复原来的确定的信号，也就必须要在附加一些很必要的条件才行。
M3.1Usingprogram3.1determineandplottherealandimaginarypartsandthemagnitudeandphasespectraofthefollowingDTFTforvariousvaluesofrandθ：
（a）程序及结果如下：
%读入离散时间傅里叶变换所需的长度
k=input（'频率点数量='）;
%读入分子和分母系数
num=input（'分子系数='）;
den=input（'分母系数'）;
%计算频率响应
w=0:
pi/（k-1）:
pi;
h=freqz（num,den,w）;
%绘制频率响应
subplot（2,2,1）
plot（w/pi,real（h））;grid
title（'实部'）
xlabel（'\omega/\pi'）;ylabel（'振幅'）
subplot（2,2,2）
plot（w/pi,imag（h））;grid
title（'虚部'）
xlabel（'\omega/\pi'）;ylabel（'振幅'）
subplot（2,2,3）
plot（w/pi,abs（h））;grid
title（'幅度谱'）
xlabel（'\omega/\pi'）;ylabel（'幅度'）
subplot（2,2,4）
plot（w/pi,angle（h））;grid
title（'相位谱'）
xlabel（'\omega/\pi'）;ylabel（'相位，弧度'）
取r=0.5，θ=60，
（b）取r=0.8，θ=90，如下：
又这两幅图中可以看出，幅度谱的变化基本和实部的变化一致，也是由于虚部的变化转折点基本和实部一致的原因，所以通过这两幅图也可看出，实部虚部与幅度的关系。
，而通过相位谱和幅度谱的关系，可以看出，幅度谱的转折点也几乎就是相位谱的转折点，同时也与虚部的变化几乎一致，它们之间的联系是非常紧密的，
，因此分析一个复数的时候，就应该把这四个方面都考虑进去，以便于更好的理解和学习。
M3.4UsingmatlabverifythefollowinggeneralpropertiesoftheDTFTaslistedinTable3.2:
（a）Linearity（b）time-shifting,（c）frequency-shifting（d）.differentiation-in-frequency.（e）convolution.（f）.modulation.and（g）Parseval`srelation.Sincealldatainmatlabhavetobefinite-lengthvectors,thesequencestobeusedtoverifythepropertiesarethusrestrictedtobeoffinitelength.
Answer:
（a）线性：
代码如下
N=input（'Thelengthofthesequence='）;
k=0:
N-1;
gamma=-0.5;
g=exp（gamma*k）;
%gisanexponentialsequence
h=sin（2*pi*k/（N/2））;
%hisasinusoidalsequencewithperiod=N/2
[G,W]=freqz（g,1,512）;
[H,W]=freqz（h,1,512）;
%property1
alpha=0.5;
beta=0.25;
y=alpha*g+beta*h;
[Y,W]=freqz（y,1,512）;
比较Y和alpha*G+beta*H即可：
可以看出两者完全一样，说明线性是完全成立的。
（b）时移
代码如下：
%property2
n0=5;%ssequenceshiftedby5simples
y2=[zeros（[1,n0]）g];
[Y2,W]=freqz（y2,1,512）;
G0=exp（-j*W*n0）.*G;
比较G0和deY2即可：
可以看出，两种方法得出来的频谱一样，说明时移性质正确。
（c）频移
代码如下:
%property3
w0=pi/2;%thevalueofcmega0=pi/2
r=256;%thevalueofomega0intermsofnumberofsamples
y3=g.*exp（j*w0*k）;
[Y3,w]=freqz（y3,1,512）;
k=0:
511;
w=-w0+pi*k/512;%creatingG（exp（w-w0））
G1=freqz（g,1,w'）;
比较G1和Y3即可：
由此可见，频移也满足。
（d）频域微分代码如下：
%property4
k=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）卷积代码如下：
%property5
y5=conv（g,h）;
[Y5,w]=freqz（y5,1,512）;
比较Y5和G.*H即可：
（f）调制
%property6
y6=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）Parseval`srelation
%property7
val1=sum（g.*conj（h））;
val2=sum（G0.*conj（H0））/512;
比较val1和val2可得如下：
由此可见，Parseval`srelation也完全符合。
这个题目，也通过仿真验证了傅里叶变换的这七个性质，也从实验验证了这几个性质的正确性。
通过6调制也可以很便利的计算个别点的调制特性。
这也便于我们更好的来理解傅里叶变换的这几个性质，也更利于我们深入的理解和学习。
其中，性质4没有理解程序意义，没有明白其原因结构，因此结果可能只是样子。
但是这个并不影响性质的正确与否。
M3.8UsingMATLABcomputetheN-pointDFTsofthelength-NsequenceofProblem3.12forN=3,5,7,and10.CompareyourresultswiththatobtainedbyevaluatingtheDTFTscomputedinProblem3.12atω=2πk/N.k=0,1,2,3,,,,N-1.
（a）.y1[n]=1,-N≤n≤N，其他为0.
代码如下：
N=input（'ThevalueofN='）;
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+1）,abs（Y1dft）,'o'）;
xlabel（'Normalizedfrequency'）;ylabel（'Amplitude'）;title（'3点DFT'）;
（b）.y2[n]=1-｜n｜/N,-N≤n≤N，其他为0.
（c）.y3[n]=cos（pi*n/2N）,-N≤n≤N，其他为0.
本道题中列举了3种序列的不同点的离散傅里叶变换的变换后的图像。
可以明确的看出，傅里叶变换的长度越大，其涵盖的内容就相应的越多，当然也就越来越精确，也就更符合原来函数的性质。
而且相应的，长度约长，其变换得到的图形高度两边的越来越高，中间的点也就越低，可能是由于点多的缘故，分散了原来所等量的能量，因此，点变多了，其幅度相应的也就要有所改变。
同时也可以看出，离散傅里叶变换具有对称特性，因此计算傅里叶的时候，可以计算一半的系数，另一半对照类比即可。
这其实也为后面的蝶形运算提供了依据。
M3.19UsingProgram3.10determinethez-transformasaratiooftwopolynomialsinz-1fromeachofthepartfractionexpansionslistedbelow.
（d）.X4（z）=4+10/（5+2z-1）+z-1/（6+5z-1+z-2）,｜z｜>0.5,
Answer:
（a）:
（b）程序代码如下：
r=input（'请输入留数='）;
请输入留数=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）看做一个系统的传递函数，这就为我们提供了计算极点和零点的快捷方法，也便于我们来将整合后的复式来分析系统的快速性和稳定性。
升级会员 