数字信号处理matlab实验第三章报告.docx
《数字信号处理matlab实验第三章报告.docx》由会员分享，可在线阅读，更多相关《数字信号处理matlab实验第三章报告.docx（26页珍藏版）》请在冰豆网上搜索。
数字信号处理matlab实验第三章报告
第三章离散时间信号的频域分析
一、实验目的
通过matlab仿真进一步研究离散时间系统的时域特性。
二、实验要求
1．学习并调试本章所给的例子。
2．回答书后给出的问题。
三、实验过程
Q.2运行程序P3.1,求离散时间傅里叶变换的实部、虚部以及幅度和相位谱。
离散时间傅里叶变换是w的周期函数吗?
若是周期是多少？
描述这四个图形表示的对称性。
clf;
w=-4*pi:
8*pi/511:
4*pi;
num=[21];den=[1-0.6];
h=freqz（num,den,w）;
%PlottheDTFT
subplot（2,2,1）
plot（w/pi,real（h））;grid
title（'H（e^{j\omega}）的实部'）
xlabel（'\omega/\pi'）;
ylabel（'振幅'）;
subplot（2,2,2）
plot（w/pi,imag（h））;grid
title（'H（e^{j\omega}）的虚部'）
xlabel（'\omega/\pi'）;
ylabel（'振幅'）;
subplot（2,2,3）
plot（w/pi,abs（h））;grid
title（'|H（e^{j\omega}）|幅度谱'）
xlabel（'\omega/\pi'）;
ylabel（'振幅'）;
subplot（2,2,4）
plot（w/pi,angle（h））;grid
title（'[H（e^{j\omega}）]相位谱'）
xlabel（'\omega/\pi'）;
ylabel（'以弧度为单位的相位'）;
由图得出是周期函数，周期为2，实部是和幅度谱是实偶对称，虚部和相位谱是实奇对称
Q.4修改程序P3.1，计算如下有限长序列的离散时间傅里叶变换：
g[n]=[1357911131517]
并重做习题Q.2。
讨论你的结果。
你能解释相位谱中的跳变吗？
clf;
w=-4*pi:
8*pi/511:
4*pi;
num=[1357911131517];den=[1];
h=freqz（num,den,w）;
%PlottheDTFT
subplot（2,2,1）
plot（w/pi,real（h））;grid
title（'H（e^{j\omega}）的实部'）
xlabel（'\omega/\pi'）;
ylabel（'振幅'）;
subplot（2,2,2）
plot（w/pi,imag（h））;grid
title（'H（e^{j\omega}）的虚部'）
xlabel（'\omega/\pi'）;
ylabel（'振幅'）;
subplot（2,2,3）
plot（w/pi,abs（h））;grid
title（'|H（e^{j\omega}）|幅度谱'）
xlabel（'\omega/\pi'）;
ylabel（'振幅'）;
subplot（2,2,4）
plot（w/pi,angle（h））;grid
title（'[H（e^{j\omega}）]相位谱'）
xlabel（'\omega/\pi'）;
ylabel（'以弧度为单位的相位'）;
1.由图得出是周期函数，周期为2，实部是和幅度谱是实偶对称，虚部和相位谱是实奇对称。
2.要计算一个系统相频特性，就要用到反正切函数，计算机中反正切函数规定，在一、二象限中的角度为0～pi，三四象限的角度为0～-pi。
若一个角度从0变到2pi，但实际得到的结果是0～pi，再由-pi～0，在w=pi处发生跳变。
Q.6通过加入合适的注释语句和程序语句，修改程序P3.2，对程序生成的图形的两个轴加标记。
哪个参数控制时移量？
clf;
w=-pi:
2*pi/255:
pi;wo=0.4*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））;grid
title（'原序列的幅度谱'）
xlabel（'\omega/\pi'）;
ylabel（'振幅'）;
subplot（2,2,2）
plot（w/pi,abs（h2））;grid
title（'时移后序列的幅度谱'）
xlabel（'\omega/\pi'）;
ylabel（'振幅'）;
subplot（2,2,3）
plot（w/pi,angle（h1））;grid
title（'原序列的相位谱'）
xlabel（'\omega/\pi'）;
ylabel（'以弧度为单位的相位'）;
subplot（2,2,4）
plot（w/pi,angle（h2））;grid
title（'时移后序列的相位谱'）
xlabel（'\omega/\pi'）;
ylabel（'以弧度为单位的相位'）;
参数D控制位移量
Q.8选取不同的时移值重做习题Q.6
令Q.6中的位移量D为20，程序结果如图：
Q.10通过加入合适的注释语句和程序语句，修改程序P3.3，对程序生成的图形的两个轴加标记。
哪个参数控制频移量？
clf;
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））;grid
title（'原序列的幅度谱'）
xlabel（'\omega/\pi'）;
ylabel（'振幅'）;
subplot（2,2,2）
plot（w/pi,abs（h2））;grid
title（'频移后序列的幅度谱'）
xlabel（'\omega/\pi'）;
ylabel（'振幅'）;
subplot（2,2,3）
plot（w/pi,angle（h1））;grid
title（'原序列的相位谱'）
xlabel（'\omega/\pi'）;
ylabel（'以弧度为单位的相位'）;
subplot（2,2,4）
plot（w/pi,angle（h2））;grid
title（'频移后序列的相位谱'）
xlabel（'\omega/\pi'）;
ylabel（'以弧度为单位的相位'）;
参数wo控制频移量
Q.12选取不同的频移值重做习题Q.10
令Q.10中wo值为π，程序结果如图
Q.14通过加入合适的注释语句和程序语句，修改程序P3.4，对程序生成的图形的两个轴加标记。
clf;
w=-pi:
2*pi/255:
pi;
x1=[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））;grid
title（'幅度谱的乘积'）
xlabel（'\omega/\pi'）;
ylabel（'振幅'）;
subplot（2,2,2）
plot（w/pi,abs（h3））;grid
title（'卷积后序列的幅度谱'）
xlabel（'\omega/\pi'）;
ylabel（'振幅'）;
subplot（2,2,3）
plot（w/pi,angle（hp））;grid
title（'相位谱的和'）
xlabel（'\omega/\pi'）;
ylabel（'以弧度为单位的相位'）;
subplot（2,2,4）
plot（w/pi,angle（h3））;grid
title（'卷积后序列的相位谱'）
xlabel（'\omega/\pi'）;
ylabel（'以弧度为单位的相位'）;
Q.16选取两个改变了长度的序列，重做习题Q.14
更改Q.14中的两个序列为X1=[1124],X2=[4231],程序结果如图：
Q.18通过加入合适的注释语句和程序语句，修改程序P3.5，对程序生成的图形的两个轴加标记。
clf;
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））;grid
title（'第一个序列的幅度谱'）
xlabel（'\omega/\pi'）;
ylabel（'振幅'）;
subplot（3,1,2）
plot（w/pi,abs（h2））;grid
title（'第二个序列的幅度谱'）
xlabel（'\omega/\pi'）;
ylabel（'振幅'）;
subplot（3,1,3）
plot（w/pi,abs（h3））;grid
title（'乘积序列的幅度谱'）
xlabel（'\omega/\pi'）;
ylabel（'振幅'）;
Q.20通过加入合适的注释语句和程序语句，修改程序P3.6，对程序生成的图形的两个轴加标记。
试解释程序怎样进行时间反转运算。
clf;
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））;grid
title（'原序列的幅度谱'）
xlabel（'\omega/\pi'）;
ylabel（'振幅'）;
subplot（2,2,2）
plot（w/pi,abs（h3））;grid
title（'时间反转后序列的幅度谱'）
xlabel（'\omega/\pi'）;
ylabel（'振幅'）;
subplot（2,2,3）
plot（w/pi,angle（h1））;grid
title（'原序列的相位谱'）
xlabel（'\omega/\pi'）;
ylabel（'以弧度为单位的相位'）;
subplot（2,2,4）
plot（w/pi,angle（h3））;grid
title（'时间反转后序列的相位谱'）
xlabel（'\omega/\pi'）;
ylabel（'以弧度为单位的相位'）;
程序先将num沿时间原点翻转，再进行离散时间傅里叶变换，然后变换后的函数乘以时移分量exp（w*L*i）。
Q.22选取两个不同长度的序列，重做习题Q.20
1.将Q.20中的序列改为num=[11245],程序结果如图：
2.将Q.20中的序列改为num=[168312],程序结果如图：
Q.24编写一个MATLAB程序，用一个N点复数离散傅里叶变换计算两个长度为N的实数序列的N点离散傅里叶变换，并将结果同直接使用两个N点离散傅里叶变换得到的结果进行比较。
clf;
g=[1124];
h=[2321];
x=g+i*h;
N=length（x）-1;n=0:
N;
gk=fft（g）;
hk=fft（h）;
xk=fft（x）;
xk1=fft（conj（x））;
gk1=（xk+xk1）/2;
hk1=（xk-xk1）/2i;
subplot（4,2,1）
stem（n,abs（gk））;grid
title（'实部序列gk的离散傅里叶变换的幅度'）
xlabel（'时间序号n'）;ylabel（'振幅'）;
subplot（4,2,2）
stem（n,abs（hk））;grid
title（'虚部序列gk的离散傅里叶变换的幅度'）
xlabel（'时间序号n'）;ylabel（'振幅'）;
subplot（4,2,3）
stem（n,abs（gk1））;grid
title（'通过xk得到的gk1的离散傅里叶变换的幅度'）
xlabel（'时间序号n'）;ylabel（'振幅'）;
subplot（4,2,4）
stem（n,abs（hk1））;grid
title（'通过xk得到的hk1的离散傅里叶变换的幅度'）
xlabel（'时间序号n'）;ylabel（'振幅'）;
subplot（4,2,5）
stem（n,angle（gk））;grid
title（'实部序列gk的离散傅里叶变换的相位'）
xlabel（'时间序号n'）;ylabel（'以弧度为单位的相位'）;
subplot（4,2,6）
stem（n,angle（hk））;grid
title（'虚部序列hk的离散傅里叶变换的相位'）
xlabel（'时间序号n'）;ylabel（'以弧度为单位的相位'）;
subplot（4,2,7）
stem（n,angle（gk1））;grid
title（'通过xk得到的gk1的离散傅里叶变换的相位'）
xlabel（'时间序号n'）;ylabel（'以弧度为单位的相位'）;
subplot（4,2,8）
stem（n,angle（hk1））;grid
title（'通过xk得到的hk1的离散傅里叶变换的相位'）
xlabel（'时间序号n'）;ylabel（'以弧度为单位的相位'）;
程序运行结果如下，说明两种方法得到的结果一样：
Q.26在函数circshift中，命令rem的作用是什么？
函数cirshift中rem是用来确定要补的零的个数。
Q.28在函数circonv中，运算符~=的作用是什么？
函数circonv中~=是用来确定两个序列的长度是否相等。
Q.30通过加入合适的注释语句和程序语句，修改程序P3.7，对程序生成的图形的两个轴加标记。
哪一个参数决定时移量？
若时移量大于序列长度，将会发生什么？
clf;
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（'原序列'）
xlabel（'时间序号n'）;ylabel（'a'）;
subplot（2,1,2）;
stem（n,b）;axis（[0,L,min（a）,max（a）]）;
title（['圆周移位',num2str（M）,'个样本得到的序列']）
xlabel（'时间序号n'）;ylabel（'a'）;
M控制时移量，若时移量大于序列长度会得到时移M1=M-序列长度的序列
下图是M=6时的结果，此时时移量M小于序列长度：
下图是M=13时的结果，此时时移量M大于序列长度：
Q.32通过加入合适的注释语句和程序语句，修改程序P3.8，对程序生成的图形的两个轴加标记。
clf;
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（'原序列的离散傅里叶变换的幅度'）;
xlabel（'时间序号n'）;ylabel（'振幅'）;
subplot（2,2,2）
stem（n,abs（YF））;grid
title（'圆周移位后的序列的离散傅里叶变换的幅度'）;
xlabel（'时间序号n'）;ylabel（'振幅'）;
subplot（2,2,3）
stem（n,angle（XF））;grid
title（'原序列的离散傅里叶变换的相位'）;
xlabel（'时间序号n'）;ylabel（'以弧度为单位的相位'）;
subplot（2,2,4）
stem（n,angle（YF））;grid
title（'圆周移位后的序列的离散傅里叶变换的相位'）;
xlabel（'时间序号n'）;ylabel（'以弧度为单位的相位'）;
时移量是5.结果如下图：
Q.34选取两个不同的时移量，重做习题Q.32
1.时移量为10，结果如下图：
2.时移量为15，结果如下图：
Q.36运行程序P3.9并验证离散傅里叶变换的圆周卷积性质
clf;
g1=[123456];g2=[1-233-21];
ycir=circonv（g1,g2）;
disp（'圆周卷积的结果='）;disp（ycir）
G1=fft（g1）;G2=fft（g2）;
yc=real（ifft（G1.*G2））;
disp（'离散傅里叶变换乘积的离散傅里叶逆变换的结果='）;disp（yc）
结果如下图：
两个结果一致，说明两个序列循环卷积等于这两个序列分别的DFT相乘后的IDFT
Q.38运行程序P3.10并验证线性卷积可通过圆周卷积得到
clf;
g1=[12345];g2=[22011];
g1e=[g1zeros（1,length（g2）-1）];
g2e=[g2zeros（1,length（g1）-1）];
ylin=circonv（g1e,g2e）;
disp（'通过圆周卷积的线性卷积='）;disp（ylin）;
y=conv（g1,g2）;
disp（'直接线性卷积='）;disp（y）
结果如下图：
两个结果一致，说明线性卷积可通过圆周卷积得到
Q.40编写一个MATLAB程序，对两个序列分别做离散傅里叶变换，以生成他们的线性卷积。
用此程序验证习题Q.38的结果
clf;
g1=[12345];g2=[22011];
g1e=[g1zeros（1,length（g2）-1）];
g2e=[g2zeros（1,length（g1）-1）];
G1e=fft（g1e）;G2e=fft（g2e）;
ylin=real（ifft（G1e.*G2e））;
disp（'通过离散傅里叶变换得到的线性卷积='）;disp（ylin）;
y=conv（g1,g2）;
disp（'直接线性卷积='）;disp（y）
结果如下图：
两者结果一致，验证了习题Q.38的结果
Q.42运行程序P3.11。
由于周期序列的偶数部分的离散傅里叶变换是原序列的XEF的实数部分，XEF的虚部应该为零。
你能验证它们吗？
你怎样解释仿真结果？
x=[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（'振幅'）;
title（'Re（DFT\{x[n]\}）'）;
subplot（2,2,2）;
plot（k/128,imag（XF））;grid;
ylabel（'振幅'）;
title（'Im（DFT\{x[n]\}）'）;
subplot（2,2,3）
plot（k/128,real（XEF））;grid;
xlabel（'时间序号n'）;ylabel（'振幅'）;
title（'Re（DFT\{x_{e}[n]\}）'）;
subplot（2,2,4）;
plot（k/128,imag（XEF））;grid;
xlabel（'时间序号n'）;ylabel（'振幅'）;
title（'Im（DFT\{x_{e}[n]\}）'）;
Q.44运行程序P3.12。
你得到的a和b的值相等吗？
clf;
x=[（1:
128）（128:
-1:
1）];
XF=fft（x）;
a=sum（x.*x）
b=round（sum（abs（XF）.^2）/256）
运行结果如图：
a和b的值相等
Q.46使用程序P3.1在单位圆上求下面的z变换：
clf;
w=-pi:
2*pi/511:
pi;
num=[25953];den=[545211];
h=freqz（num,den,w）;
%PlottheDTFT
subplot（2,2,1）
plot（w/pi,real（h））;grid
title（'H（e^{j\omega}）的实部'）
xlabel（'\omega/\pi'）;
ylabel（'振幅'）;
subplot（2,2,2）
plot（w/pi,imag（h））;grid
title（'H（e^{j\omega}）的虚部'）
xlabel（'\omega/\pi'）;
ylabel（'振幅'）;
subplot（2,2,3）
plot（w/pi,abs（h））;grid
title（'|H（e^{j\omega}）|幅度谱'）
xlabel（'\omega/\pi'）;
ylabel（'振幅'）;
subplot（2,2,4）
plot（w/pi,angle（h））;grid
title（'[H（e^{j\omega}）]相位谱'）
xlabel（'\omega/\pi'）;
ylabel（'以弧度为单位的相位'）;
Q.48通过下面程序产生的极零点图，求出G（z）的收敛域的数目。
清楚地显示出所有的收敛域。
有极零点图，说明离散时间傅里叶变换是否存在
clf;
w=-pi:
2*pi/511:
pi;
num=[25953];den=[545211];
[z,p,k]=tf2zp（num,den）;
sos=zp2sos（z,p,k）;
zplane（z,p）;
disp（'G（z）的零点='）;disp（z）
disp（'G（z）的极点='）;disp（p）
disp（'G（z）的因式形式='）;disp（sos）
极零点图如下:
有一个收敛域，收敛域包含单位圆，则离散时间傅里叶变换存在
Q.50编写一个MATLAB程序，计算一个有理逆Z变换的前L个
升级会员 