数字信号处理实验吴镇扬版matlab程序.docx
《数字信号处理实验吴镇扬版matlab程序.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验吴镇扬版matlab程序.docx（34页珍藏版）》请在冰豆网上搜索。
数字信号处理实验吴镇扬版matlab程序
（1）数组的加、减、乘、除和乘方运算。
输入A=[1234]，B=[3456]，求C=A+B，D=A-B，E=A.*B，F=A./B，G=A.^B并用stem语句画出A、B、C、D、E、F、G。
clearall;
a=[1234];
b=[3456];
c=a+b;
d=a-b;
e=a.*b;
f=a./b;
g=a.^b;
n=1:
4;
subplot（4,2,1）;stem（n,a）;
xlabel（'n'）;xlim（[05]）;ylabel（'A'）;
subplot（4,2,2）;stem（n,b）;
xlabel（'n'）;xlim（[05]）;ylabel（'B'）;
subplot（4,2,3）;stem（n,c）;
xlabel（'n'）;xlim（[05]）;ylabel（'C'）;
subplot（4,2,4）;stem（n,d）;
xlabel（'n'）;xlim（[05]）;ylabel（'D'）;
subplot（4,2,5）;stem（n,e）;
xlabel（'n'）;xlim（[05]）;ylabel（'E'）;
subplot（4,2,6）;stem（n,f）;
xlabel（'n'）;xlim（[05]）;ylabel（'F'）;
subplot（4,2,7）;stem（n,g）;
xlabel（'n'）;xlim（[05]）;ylabel（'G'）;
（2）用MATLAB实现下列序列：
a）x（n）=0.8n0≤n≤15
b）x（n）=e（0.2+3j）n0≤n≤15
c）x（n）=3cos（0.125πn+0.2π）+2sin（0.25πn+0.1π）0≤n≤15
d）将c）中的x（n）扩展为以16为周期的函数x16（n）=x（n+16）,绘出四个周期。
e）将c）中的x（n）扩展为以10为周期的函数x10（n）=x（n+10）,绘出四个周期。
clearall;
N=0:
15;
xa=0.8.^N;
figure;subplot（2,1,1）;stem（N,xa）;xlabel（'n'）;xlim（[016]）;ylabel（'xa'）;
xb=exp（（0.2+3*j）*N）;
subplot（2,1,2）;stem（N,xb）;
xlabel（'n'）;xlim（[016]）;ylabel（'xb'）;figure;
xc=3*cos（0.125*pi*N+0.2*pi）+2*sin（0.25*pi*N+0.1*pi）;
subplot（3,1,1）;stem（N,xc）;xlabel（'n'）;xlim（[016]）;ylabel（'xc'）;
k=0:
3;m=0;
fori=1:
4
forj=1:
16
m=m+1;
n（m）=N（j）+16*k（i）;
x16（m）=3*cos（0.125*pi*n（m）+0.2*pi）+2*sin（0.25*pi*n（m）+0.1*pi）;
end
end
subplot（3,1,2）;stem（n,x16）;xlabel（'n'）;ylabel（'x16'）;
forj=1:
10
x10（j）=x16（j）;
end
fori=1:
3
form=1:
10
x10（i*10+m）=x10（m）;
end
end
n=1:
40;
subplot（3,1,3）;stem（n,x10）;xlabel（'n'）;ylabel（'x10'）;
（3）x（n）=[1,-1,3,5],产生并绘出下列序列的样本:
a）x1（n）=2x（n+2）-x（n-1）-2x（n）
b）
clearall
n=1:
4;
T=4;
x=[1-135];
x（5:
8）=x（1:
4）;
subplot（2,1,1）;stem（1:
8,x）;grid;
fori=1:
4
ifi-1<0
x1（i）=2*x（i+2）-x（i-1）-2*x（i）;
else
x1（i）=2*x（i+2）-x（i-1+T）-2*x（i）;
end
end
x1（5:
8）=x1（1:
4）;
subplot（2,1,2）;stem（1:
8,x1）;grid;
（4）绘出下列时间函数的图形，对x轴、y轴以及图形上方均须加上适当的标注：
a）x（t）=sin（2πt）0≤t≤10s
b）x（t）=cos（100πt）sin（πt）0≤t≤4s
ta=0:
0.05:
10;
xa=sin（2*pi*ta）;
subplot（2,1,1）;plot（ta,xa）;
xlabel（'t'）;ylabel（'幅度'）;
tb=0:
0.01:
4;
xb=cos（100*pi*tb）.*sin（pi*tb）;
subplot（2,1,2）;plot（tb,xb）;
xlabel（'t'）;ylabel（'幅度'）;
（5）编写函数stepshift（n0,n1,n2）实现u（n-n0）,n1n0=5;ns=1;nf=10;%ns为起点;nf为终点；在=n=n0处生成单位阶跃序列
n=[ns:
nf];
x=[（n-n0）>=0];
stem（n,x）;
（6）给一定因果系统
求出并绘制H（z）的幅频响应与相频响应。
clearall;
b=[1,sqrt
（2）,1];
a=[1,-0.67,0.9];
[h,w]=freqz（b,a）;
am=20*log10（abs（h））;
subplot（2,1,1）;plot（w,am）;
ph=angle（h）;
subplot（2,1,2）;plot（w,ph）;
（7）计算序列{8-2-123}和序列{23-1-3}的离散卷积，并作图表示卷积结果。
clearall;
a=[8-2-123];
b=[23-1-3];
c=conv（a,b）;%计算卷积
M=length（c）-1;
n=0:
1:
M;
stem（n,c）;
xlabel（'n'）;ylabel（'幅度'）;
（8）求以下差分方程所描述系统的单位脉冲响应h（n）,0≤n≤50
y（n）+0.1y（n-1）-0.06y（n-2）=x（n）-2x（n-1）
clearall;
N=50;
a=[1-2];
b=[10.1-0.06];
x=[1zeros（1,N-1）];
k=0:
1:
N-1;
y=filter（a,b,x）;
stem（k,y）;
xlabel（'n'）;ylabel（'幅度'）;
（1）、观察高斯序列的时域和幅频特性，固定信号xa（n）中参数p=8，改变q的值，使q分别等于2，4，8，观察它们的时域和幅频特性，了解当q取不同值时，对信号序列的时域幅频特性的影响；固定q=8，改变p，使p分别等于8，13，14，观察参数p变化对信号序列的时域及幅频特性的影响，观察p等于多少时，会发生明显的泄漏现象，混叠是否也随之出现？
记录实验中观察到的现象，绘出相应的时域序列和幅频特性曲线。
解：
程序见附录程序一：
n=0:
1:
15;
%p=8不变，q变化（2,4,8）;
p=8;q=2;%p=8;q=2;
xa1=exp（-（（n-p）.^2）/q）;
subplot（5,2,1）;
plot（n,xa1,'-*'）;
xlabel（'t/T'）;
ylabel（'xa（n）'）;
title（'p=8q=2'）
xk1=abs（fft（xa1））;
subplot（5,2,2）;
stem（n,xk1）
xlabel（'k'）;
ylabel（'Xa（k）'）;
title（'p=8q=2'）
p=8;q=4;%p=8;q=4;
xa1=exp（-（（n-p）.^2）/q）;
subplot（5,2,3）;
plot（n,xa1,'-*'）;
xlabel（'t/T'）;
ylabel（'xa（n）'）;
title（'p=8q=4'）
xk1=abs（fft（xa1））;
subplot（5,2,4）;
stem（n,xk1）
xlabel（'k'）;
ylabel（'Xa（k）'）;
title（'p=8q=4'）
p=8;q=8;%p=8;q=8;
xa1=exp（-（（n-p）.^2）/q）;
subplot（5,2,5）;
plot（n,xa1,'-*'）;
xlabel（'t/T'）;
ylabel（'xa（n）'）;
xk1=abs（fft（xa1））;
title（'p=8q=8'）
subplot（5,2,6）;
stem（n,xk1）
xlabel（'k'）;
ylabel（'Xa（k）'）;
title（'p=8q=8'）
%q=8不变,p变化（8,13,14）;
p=8;q=8;%p=8;q=8;
xa1=exp（-（（n-p）.^2）/q）;
subplot（5,2,5）;
plot（n,xa1,'-*'）;
xlabel（'t/T'）;
ylabel（'xa（n）'）;
xk1=abs（fft（xa1））;
title（'p=8q=8'）
subplot（5,2,6）;
stem（n,xk1）
xlabel（'k'）;
ylabel（'Xa（k）'）;
title（'p=8q=8'）
p=13;q=8;%p=13;q=8;
xa1=exp（-（（n-p）.^2）/q）;
subplot（5,2,7）;
plot（n,xa1,'-*'）;
xlabel（'t/T'）;
ylabel（'xa（n）'）;
xk1=abs（fft（xa1））;
title（'p=13q=8'）
subplot（5,2,8）;
stem（n,xk1）
xlabel（'k'）;
ylabel（'Xa（k）'）;
title（'p=13q=8'）
p=14;q=8;%p=14;q=8;
xa1=exp（-（（n-p）.^2）/q）;
subplot（5,2,9）;
plot（n,xa1,'-*'）;
xlabel（'t/T'）;
ylabel（'xa（n）'）;
title（'p=14q=8'）
xk1=abs（fft（xa1））;
subplot（5,2,10）;
stem（n,xk1）
xlabel（'k'）;
ylabel（'Xa（k）'）;
title（'p=14q=8’）
分析：
由高斯序列表达式知n=p为期对称轴；
当p取固定值时，时域图都关于n=8对称截取长度为周期的整数倍，没有发生明显的泄漏现象；但存在混叠，当q由2增加至8过程中，时域图形变化越来越平缓，中间包络越来越大，可能函数周期开始增加，频率降低，渐渐小于fs/2,混叠减弱；
当q值固定不变，p变化时，时域对称中轴右移，截取的时域长度渐渐地不再是周期的整数倍，开始无法代表一个周期，泄漏现象也来越明显，因而图形越来越偏离真实值，
p=14时的泄漏现象最为明显，混叠可能也随之出现；
（2）、观察衰减正弦序列xb（n）的时域和幅频特性，a=0.1，f=0.0625，检查谱峰出现位置是否正确，注意频谱的形状，绘出幅频特性曲线，改变f，使f分别等于0.4375和0.5625，观察这两种情况下，频谱的形状和谱峰出现位置，有无混叠和泄漏现象？
说明产生现象的原因。
n1=0:
1:
15;
xb1=exp（-0.1*n1）.*sin（2*pi*0.0625*n1）;
subplot（3,2,1）;
plot（n1,xb1,'-*'）;
xlabel（'n'）;
ylabel（'x（n）'）;
title（'f=0.0625'）;
xk1=abs（fft（xb1））;
subplot（3,2,2）;
stem（n1,xk1）
xlabel（'k'）;
ylabel（'X（k）'）;
title（'f=0.0625'）;
n2=0:
1:
15;
xb2=exp（-0.1*n2）.*sin（2*pi*0.4375*n2）;
subplot（3,2,3）;
plot（n2,xb2,'-*'）;
xlabel（'n'）;
ylabel（'x（n）'）;
title（'f=0.4375'）;
xk2=abs（fft（xb2））;
subplot（3,2,4）;
stem（n2,xk2）
xlabel（'k'）;
ylabel（'X（k）'）;
title（'f=0.4375'）;
n3=0:
1:
15;
xb3=exp（-0.1*n3）.*sin（2*pi*0.5625*n3）;
subplot（3,2,5）;
plot（n3,xb3,'-*'）;
xlabel（'n'）;
ylabel（'x（n）'）;
title（'f=0.5625'）;
xk3=abs（fft（xb3））;
subplot（3,2,6）;
stem（n3,xk3）
xlabel（'k'）;
ylabel（'X（k）'）;
title（'f=0.5625'）;
分析：
当f=f1=0.0625时，谱峰位置出现正确，存在在混叠现象，时域采样为一周期，不满足采样定理。
当f=0.4375和0.5625时，时域图像关于Y轴对称，频域完全相同。
这是因为频域图是取绝对值的结果，所以完全相同。
另外由于时域采样为6个半周期，满足采样定理，无混叠；但由于截取长度不是周期整数倍，出现泄漏。
（3）、观察三角波和反三角波序列的时域和幅频特性，用N=8点FFT分析信号序列xc（n）和xd（n）的幅频特性，观察两者的序列形状和频谱曲线有什么异同？
绘出两序列及其幅频特性曲线。
在xc（n）和xd（n）末尾补零，用N=16点FFT分析这两个信号的幅频特性，观察幅频特性发生了什么变化？
两情况的FFT频谱还有相同之处吗？
这些变化说明了什么？
三角波序列:
反三角波序列：
clearall;
n=[0:
3];k=[1:
8];
%定义三角波序列
Xc（n+1）=n;Xc（n+5）=4-n;
%定义反三角波序列
Xd（n+1）=4-n;Xd（n+5）=n;
%%%%%%%%%%%三角波特性%%%%%%%%%%%%%
subplot（2,2,1）;plot（k-1,Xc）;
xlabel（'n'）;ylabel（'时域特性'）;text（1,3,'三角波'）;
subplot（2,2,2）;plot（k-1,abs（fft（Xc）））;
xlabel（'k'）;ylabel（'幅频特性'）;text（4,10,'三角波'）;
%%%%%%%%%%%反三角波特性%%%%%%%%%%
subplot（2,2,3）;plot（k-1,Xd）;
xlabel（'n'）;ylabel（'时域特性'）;text（3,3,'反三角波'）;
subplot（2,2,4）;plot（k-1,abs（fft（Xd）））;
xlabel（'k'）;ylabel（'幅频特性'）;text（4,10,'反三角波'）;
%末尾补0，计算32点FFT
Xc（9:
32）=0;Xd（9:
32）=0;k=1:
32;figure;
%%%%%%%%%%%三角波特性%%%%%%%%%%%%%
subplot（2,2,1）;plot（k-1,Xc）;
xlabel（'n'）;ylabel（'时域特性'）;text（1,3,'三角波'）;
subplot（2,2,2）;plot（k-1,abs（fft（Xc）））;
xlabel（'k'）;ylabel（'幅频特性'）;text（4,10,'三角波'）;
%%%%%%%%%%%fan三角波特性%%%%%%%%%%
subplot（2,2,3）;plot（k-1,Xd）;
xlabel（'n'）;ylabel（'时域特性'）;text（3,3,'反三角波'）;
subplot（2,2,4）;plot（k-1,abs（fft（Xd）））;
xlabel（'k'）;ylabel（'幅频特性'）;text（4,10,'反三角波'）;
N=8时域和幅度频谱图：
分析：
由图知，三角波序列和反三角波序列的时域图像成镜像关系，但频域图像完全一样，只是因为幅频图是对x（k）的值取绝对值。
N=32时域和幅度频谱图：
分析：
由实验所得的图形知，N=32点时
和
的幅频特性都更加密集，更多离散点的幅值显示，“栅栏效应”减小，分辨率提高，而对于
来说变化更加明显。
在原序列的末端填补零值，变动了DFT的点数，人为的改变了对真实频谱采样的点数和位置，相当于搬动了“尖桩栅栏”的位置，从而使得频谱的峰点和谷点暴露出来。
N=32时，
和
的频谱差别较大，但总体趋势仍然都是中间最小，两侧呈对称。
（4）、一个连续信号含两个频率分量，经采样得
x（n）=sin2π*0.125n+cos2π*（0.125+Δf）nn=0,1……,N-1
已知N=16，Δf分别为1/16和1/64，观察其频谱；当N=128时，Δf不变，其结果有何不同，为什么？
clearall;
%%%%%%%N=16%%%%%%%%%%%%
N=16;detf=1/16;n=[0:
N-1];
x1（n+1）=sin（2*pi*0.125.*n）+cos（2*pi*（0.125+detf）.*n）;
detf=1/64;x2（n+1）=sin（2*pi*0.125.*n）+cos（2*pi*（0.125+detf）.*n）;
%%%%%%%N=16,detf=1/16%%%%%%%%%%%%
subplot（2,2,1）;stem（n,x1）;hold;plot（n,x1）;
xlabel（'n'）;ylabel（'时域特性'）;text（6,1,'N=16,detf=1/16'）;
subplot（2,2,2）;stem（n,abs（fft（x1）））;
xlabel（'n'）;ylabel（'幅值特性'）;text（6,4,'N=16,detf=1/16'）;
%%%%%%%N=16,detf=1/64%%%%%%%%%%%%
subplot（2,2,3）;stem（n,x2）;
xlabel（'n'）;ylabel（'时域特性'）;text（6,1,'N=16,detf=1/64'）;
subplot（2,2,4）;stem（n,abs（fft（x2）））;
xlabel（'n'）;ylabel（'幅值特性'）;text（6,4,'N=16,detf=1/64'）;
%%%%%%%N=128%%%%%%%%%%%%
N=128;detf=1/16;n=[0:
N-1];
x3（n+1）=sin（2*pi*0.125.*n）+cos（2*pi*（0.125+detf）.*n）;
detf=1/64;
x4（n+1）=sin（2*pi*0.125.*n）+cos（2*pi*（0.125+detf）.*n）;
%%%%%%%N=128,detf=1/16%%%%%%%%%%%%
figure;subplot（2,2,1）;stem（n,x3）;
xlabel（'n'）;ylabel（'时域特性'）;axis（[0128-22]）;text（6,1.5,'N=128,detf=1/16'）;
subplot（2,2,2）;stem（n,abs（fft（x3）））;
xlabel（'n'）;ylabel（'幅值特性'）;axis（[0128-1070]）;text（40,60,'N=128,detf=1/16'）;
%%%%%%%N=128,detf=1/64%%%%%%%%%%%%
subplot（2,2,3）;stem（n,x3）;
xlabel（'n'）;ylabel（'时域特性'）;axis（[0128-22]）;text（6,1.5,'N=128,detf=1/16'）;
subplot（2,2,4）;stem（n,abs（fft（x4）））;
xlabel（'n'）;ylabel（'幅值特性'）;axis（[0128-1070]）;text（40,60,'N=128,detf=1/16'）;
（5）、用FFT分别实现xa（n）（p＝8，q＝2）和xb（n）（a＝0.1，f＝0.0625）的16点圆周卷积和线性卷积。
clearall;
N=16;
n=0:
N-1;
p=8;q=2;
Xa（n+1）=exp（-（n-p）.^2./q）;
a=0.1;f=0.0625;
Xb（n+1）=exp（-a.*n）.*sin（2*pi*f.*n）;
%16点循环卷积
Fa=fft（Xa）;Fb=fft（Xb）;
Fx=Fa.*Fb;
X51=ifft（Fx）;
stem（n,X51）;
%16点线性卷积
Xa（N+1:
2*N-1）=0;
Xb（N+1:
2*N-1）=0;
Fa=fft（Xa）;Fb=fft（Xb）;
Fc=Fa.*Fb;
X52=ifft（Fc）;
figure;stem（1:
2*N-1,X52）;
（7）用FFT分别计算xa（n）（p=8,q=2）和xb（n）（a=0.1,f=0.0625）的16点循环相关和线性相关，问一共有多少种结果，他们之间有何异同点。
clearall;
N=16;
n=0:
N-1;
p=8;q=2;
Xa（n+1）=exp（-（n-p）.^2./q）;
a=0.1;f=0.0625;
Xb（n+1）=exp（-a.*n）.*sin（2*pi*f.*n）;
N=length（Xa）;
%16点循环相关
Fa=fft（Xa,2*N）;Fb=fft（Xb,2*N）;
Fx=conj（Fa）.*Fb;
X71=real（ifft（Fx））;
X71=[X71（N+2:
2*N）X71（1:
N）];
n=（-N+1）:
（N-1）;stem（n,X71）;
%16点线性相关
Xa（N+1:
2*N-1）=0;
Xb（N+1:
2*N-1）=0;
Fa=fft（Xa）;Fb=fft（Xb）;
Fc=conj（Fa）.*Fb;
X72=real（ifft（Fc））;
figure;stem（1:
2*N-1,X72）;
（8）用FFT分别计算xa（n）（p=8,q=2）和xb（n）（a=0.1,f=0.0625）的自相关函数。
clearall;
N=16;
n=0:
N-1;
p=8;q=2;
Xa（n+1）=exp（-（n-p）.^2./q）;
a=0.1;f=0.0625;
Xb（n+1）=exp（-a.*n）.*sin（2*pi*f.*n）;%自然对数的底:
e=:
2.71828182845904523536
N=length（Xa）;
%Xa（n）16点自相关
Fa=fft（Xa,2*N）;Fb=fft（X
升级会员 