DSP实验Matlab程序实例.docx
《DSP实验Matlab程序实例.docx》由会员分享，可在线阅读，更多相关《DSP实验Matlab程序实例.docx（30页珍藏版）》请在冰豆网上搜索。
DSP实验Matlab程序实例
2．离散时间信号（序列）的产生
利用MATLAB产生和绘制下列有限长序列：
单位脉冲序列
，单位阶跃序列
，矩形序列
x=-10:
10;
y=[];
fori=1:
21
ifx（i）==0
y（i）=1;
elsey（i）=0;
end
end
subplot（3,1,1）;
%在子图画单位冲激序列
stem（x,y,'*'）;
title（'冲激'）;
xlabel（'n'）;
holdon;
fori=1:
21
ifx（i）<0
y（i）=0;
elsey（i）=1;
end
end
subplot（3,1,2）;
%在子图画单位阶跃序列
stem（x,y,'*'）;
title（'阶跃'）;
xlabel（'n'）;
holdon;
fori=1:
21
ifx（i）<-4
y（i）=0;
elseifx（i）>4
y（i）=0;
elsey（i）=1;
end
end
end
subplot（3,1,3）;
%在子图画矩形脉冲序列
stem（x,y,'*'）;
title（'矩形'）;
xlabel（'n'）;
holdon
在一幅图上绘出曲线
，
，A=2,α=0.5,f=2Hz
，A1=1,A2=0.5,A3=0.2，f=2Hz。
t=[0:
0.5:
720]*pi/180;
x1=5*sin（pi*t/5+pi/3）;
x2=5*cos（t-pi/3）;
x3=2*sin（3*t）.*cos（2*t-pi/3）;
subplot（3,1,1）;%在子图中画出一系列正弦余弦曲线
plot（t,x1,'-red',t,x2,'-gr',t,x3,'-bl'）;
title（'正弦余弦曲线'）;
xlabel（'t'）;
holdon;
A=2;
a=0.5;
f=2;
x4=A*exp（-a*t）.*sin（2*pi*f*t）;
subplot（3,1,2）;%在子图中画出正弦衰减信号
plot（t,x4）
title（'正弦衰减信号'）;
xlabel（'t'）;
holdon;
A1=1;
A2=0.5;
A3=0.2;
f=2;
x5=A1*sin（2*pi*f*t）+A2*sin（2*pi*2*f*t）+A3*sin（2*pi*3*f*t）;
subplot（3,1,3）;%在子图中画出谐波信号
plot（t,x5）
title（'谐波信号'）;
xlabel（'t'）;
holdon;
3．序列的运算
生成下列序列：
（1）利用MATLAB编程完成上述两序列的卷积，并绘制运算后序列的波形。
x=[12345];
h=[1212];
n1=[01234];
n2=[2345];
y=conv（x,h）;%计算序列x与h的卷积和序列y
n0=n1
（1）+n2
（1）;%计算卷积和序列y的起点位置
n3=length（x）+length（h）-2+n0;%计算卷积和序列y的终止位置
n=n0:
1:
n3;%确定卷积和y非零样值的时间向量
subplot（3,1,1）;
stem（n1,x,'.'）;%在子图绘x（n）
title（'x（n）'）;
xlabel（'n'）;
ylabel（'x（n）'）;
subplot（3,1,2）;
stem（n2,h,'.'）;%在子图绘h（n）
title（'h（n）'）;
xlabel（'n'）;
ylabel（'h（n）'）;
subplot（3,1,3）;
stem（n,y,'.'）;%在子图绘卷积和y（n）
title（'y（n）=x（n）*h（n）'）;
xlabel（'n'）;
ylabel（'y（n）'）;
（2）编写一个完成卷积的函数，输入是两个序列，输出是结果。
function[y,n]=jsjuanji（x,n1,h,n2）
%计算序列卷积和y（n）=x（n）*h（n）
%y:
卷积和y（n）对应的非零样值向量
%n：
卷积和y（n）对应的时间向量
%x:
x（n）对应的非零样值向量
%n1:
x（n）对应的时间向量
%h:
h（n）对应的非零样值向量
%n2:
h（n）对应的时间向量
n3=n1
（1）+n2
（1）;%卷积和序列y的起始位置
n4=length（x）+length（h）-1;%卷积和序列y的长度
m=0;
fori=1:
n4%求卷积和y
ifi>length（x）
x（i）=0;
end
forj=1:
i
ifj>length（h）
h（j）=0;
end
m=m+x（i+1-j）*h（j）;
end
y（i）=m;
m=0;
end
n1=n1
（1）:
1:
n1
（1）+n4-1;n1:
%x（n）对应的新的时间向量
n2=n2
（1）:
1:
n2
（1）+n4-1;%h（n）对应的新的时间向量
n=n3:
1:
n3+n4-1;%y（n）对应的时间向量
subplot（3,1,1）;
stem（n1,x,'.'）;%在子图绘x（n）
title（'x（n）'）;
xlabel（'n'）;
ylabel（'x（n）'）;
subplot（3,1,2）;
stem（n2,h,'.'）;%在子图绘h（n）
title（'h（n）'）;
xlabel（'n'）;
ylabel（'h（n）'）;
subplot（3,1,3）;
stem（n,y,'.'）;%在子图绘卷积和y（n）
title（'y（n）=x（n）*h（n）'）;
xlabel（'n'）;
ylabel（'y（n）'）;
4．采样定理的研究
分别令采样周期Ts为不同值，绘出不同采样周期下x（t）=sin（t）的频谱（直接用FFT函数），观察频谱混叠现象和防止混叠的采样周期。
fs=2;%设定采样频率
T=1/fs;%采样周期
N=128;%采样点数
n=0:
N-1;
t=n*T;%采样时间点
x=sin（t）;%生成正弦信号
subplot（4,2,1）;
plot（t,x）;%画正弦信号的时域波形
xlabel（'t'）;
ylabel（'x'）;
title（'正弦信号x=sin（t）时域波形'）;
y=fft（x,N）;%进行fft变换
mag=sqrt（y.*conj（y））;%求幅值
f=（0:
N-1）*fs/N;%进行对应的频率转换
subplot（4,2,2）;
plot（f,mag）;%画频谱图
axis（[-10,50,0,80]）;
xlabel（'频率'）;
ylabel（'幅值'）;
title（'正弦信号x=sin（t）幅频谱图fs=2'）;
holdon;
fs=10;%设定采样频率
T=1/fs;%采样周期
N=1024;%采样点数
n=0:
N-1;
t=n*T;%采样时间点
x=sin（t）;%生成正弦信号
subplot（4,2,3）;
plot（t,x）;%画正弦信号的时域波形
xlabel（'t'）;
ylabel（'x'）;
title（'正弦信号x=sin（t）时域波形'）;
y=fft（x,N）;%进行fft变换
mag=sqrt（y.*conj（y））;%求幅值
f=（0:
N-1）*fs/N;%进行对应的频率转换
subplot（4,2,4）;
plot（f,mag）;%画频谱图
axis（[-10,50,0,400]）;
xlabel（'频率'）;
ylabel（'幅值'）;
title（'正弦信号x=sin（t）幅频谱图fs=10'）;
holdon;
fs=20;%设定采样频率
T=1/fs;%采样周期
N=1024;%采样点数
n=0:
N-1;
t=n*T;%采样时间点
x=sin（t）;%生成正弦信号
subplot（4,2,5）;
plot（t,x）;%画正弦信号的时域波形
xlabel（'t'）;
ylabel（'x'）;
title（'正弦信号x=sin（t）时域波形'）;
y=fft（x,N）;%进行fft变换
mag=sqrt（y.*conj（y））;%求幅值
f=（0:
N-1）*fs/N;%进行对应的频率转换
subplot（4,2,6）;
plot（f,mag）;%画频谱图
axis（[-10,50,0,400]）;
xlabel（'频率'）;
ylabel（'幅值'）;
title（'正弦信号x=sin（t）幅频谱图fs=20'）;
holdon;
fs=30;%设定采样频率
T=1/fs;%采样周期
N=1024;%采样点数
n=0:
N-1;
t=n*T;%采样时间点
x=sin（t）;%生成正弦信号
subplot（4,2,7）;
plot（t,x）;%画正弦信号的时域波形
xlabel（'t'）;
ylabel（'x'）;
title（'正弦信号x=sin（t）时域波形'）;
y=fft（x,N）;%进行fft变换
mag=sqrt（y.*conj（y））;%求幅值
f=（0:
N-1）*fs/N;%进行对应的频率转换
subplot（4,2,8）;
plot（f,mag）;%画频谱图
axis（[-10,50,0,400]）;
xlabel（'频率'）;
ylabel（'幅值'）;
title（'正弦信号x=sin（t）幅频谱图fs=30'）;
五、思考题
1．如何产生方波信号序列和锯齿波信号序列？
答：
方波信号由square函数可得，锯齿波信号由循环语句产生。
%方波
T=0:
0.001:
2*pi;%方波的时间向量
y=square（6*T）;%方波w=6,周期为T=2*pi/6
subplot（2,1,1）;
plot（T,y）;%画方波
title（'方波'）;
axis（[T
（1）-1T（end）+1-22]）;%坐标轴区域
holdon
%锯齿波
x=-10:
10;
y=[];
fori=1:
21%rem为取余函数，此循环使得y（i）交替取值1和-1
ifrem（i,2）==0
y（i）=-1;
else
y（i）=1;
end
end
subplot（2,1,2）;
plot（x,y）;
title（'锯齿波'）;
1编写一个DFT计算程序，要求该程序具有正变换及反变换的双重功能。
令x（n）=sin（n）n（-1）n,n=0,1,….,31,这是一个震荡信号，由它来验证所编程序是否正确（即：
对x（n）作DFT得X（k），再对X（k）作反变换，检验是否得到同样的x（n））。
function[X,x2]=bianhuan（x,N）
%DFT
X=zeros（1,N）;
form=1:
N
forn=1:
N
X（m）=X（m）+x（n）*exp（-j*2*pi*n*m/N）;
end
end
%IDFT
x2=zeros（1,N）;
form=1:
N
forn=1:
N
x2（m）=x2（m）+X（n）*exp（j*2*pi*n*m/N）;
end
x2（m）=x2（m）/N;
end
clear
N=32;
x=[];
forp=1:
N
x（p）=sin（p）*p*（-1）^p
end
[X,x2]=bianhuan（x,N）
n=0:
1:
N-1;
subplot（3,1,1）;stem（n,x,'*'）;title（'original'）;
subplot（3,1,2）;stem（n,X,'*'）;title（'afterDFT'）;
subplot（3,1,3）;stem（n,x2,'*'）;title（'afterIDFT'）;
2.序列x（n）的N点DTFT的物理意义是对X（ejω）的在[0，2π]上进行N点等间隔采样
clear;
N1=1024;
w1=（0:
2*pi/1024:
2*pi）;
X1=（1-exp（-j*4*w1））./（1-exp（-j*w1））;
N2=8;
fori=1:
N2
X2（i）=X1（i*N1/N2）;
w2（i）=2*pi*i/N2;
end
N3=16;
fori=1:
N3
X3（i）=X1（i*N1/N3）;
w3（i）=2*pi*i/N3;
end
subplot（3,2,1）;plot（w1,abs（X1））;title
subplot（3,2,2）,plot（w1,angle（X1））;title
subplot（3,2,3）;stem（w2,abs（X2））;title
subplot（3,2,4）,stem（w2,angle（X2））;title
subplot（3,2,5）;stem（w3,abs（X3））;title
subplot（3,2,6）,stem（w3,angle（X3））;title
4.混叠现象
研究衰减正弦信号x（t）=[1+sin（7πf0t）]cos（2πf0t）,f0=120Hz,fs=200Hz,采样点数N=64，用DFT计算信号的幅频和相频，观察混叠现象
clear;
f0=120;
N=64;
fs1=200;
t1=（0:
1:
N-1）/fs1;
x1=（1+sin（7*pi*f0*t1））.*cos（2*pi*f0*t1）;
y1=fft（x1）;
fs2=1500;
t2=（0:
1:
N-1）/fs2;
x2=（1+sin（7*pi*f0*t2））.*cos（2*pi*f0*t2）;
y2=fft（x2）;
subplot（2,1,1）;stem（abs（y1））;title（'fs=200'）;
subplot（2,1,2）;stem（abs（y2））;title（'fs=1500'）;
5.泄漏现象
令余弦信号x（t）=cos（2πf0t）,f0=60Hz,fs=200Hz,采样点数N=64
要求：
先用DFT计算采样信号x（n）的频谱X（k），显示结果图形。
然后将采样信号x（n）乘上哈明窗函数：
x4（t）=0.54-0.46cos（2πn/N）,n=0,1,2,……,N-1
用DFT计算加窗后信号的频谱XW（k），分别在屏幕上显示信号幅频X（k），窗函数幅频X4（k）及加窗后信号的频谱XW（k），比较频谱X（k）与XW（k）的不同，并给出解释。
clear;
f0=60;
fs=200;
N=64;
t=（0:
1:
N-1）/fs;
x=cos（2*pi*f0*t）;
X=fft（x）;
n=0:
1:
N-1;
x1=0.54-0.46.*cos（2*pi*n/N）;
X1=fft（x1）;
x2=x.*x1;
X2=fft（x2）;
subplot（3,1,1）;stem（abs（X））;title；
subplot（3,1,2）;stem（abs（X1））;title;
subplot（3,1,3）;stem（abs（X2））;title;
6.栅栏效应
将内容1中的信号补零至L,L分别为16，32，64，计算x′（n）的频谱X′（k），并与X（k）比较，将频谱图X′（k）与X（k）显示出来，观察补零的效果。
clear;
N=8;
x=[];
forp=1:
N
x（p）=sin（p）*p*（-1）^p
end
X=zeros（1,N）;
forp=1:
N
forq=1:
N
X（p）=X（p）+x（q）*exp（-j*2*pi*p*q/N）;
end
end
L1=16;
x1=[];
forp=1:
L1
ifp<=8
x1（p）=sin（p）*p*（-1）^p
elsex1（p）=0
end
end
X1=zeros（1,L1）;
forp=1:
L1
forq=1:
L1
X1（p）=X1（p）+x1（q）*exp（-j*2*pi*p*q/L1）;
end
end
L2=32;
x2=[];
forp=1:
L2
ifp<=8
x2（p）=sin（p）*p*（-1）^p
elsex2（p）=0
end
end
X2=zeros（1,L2）;
forp=1:
L2
forq=1:
L2
X2（p）=X2（p）+x2（q）*exp（-j*2*pi*p*q/L2）;
end
end
L3=64;
x3=[];
forp=1:
L3
ifp<=8
x3（p）=sin（p）*p*（-1）^p
elsex3（p）=0
end
end
X3=zeros（1,L3）;
forp=1:
L3
forq=1:
L3
X3（p）=X3（p）+x3（q）*exp（-j*2*pi*p*q/L3）;
end
end
n=0:
1:
N-1;
n1=0:
1:
L1-1;
n2=0:
1:
L2-1;
n3=0:
1:
L3-1;
subplot（4,1,1）;stem（n,X,'*'）;title;
subplot（4,1,2）;stem（n1,X1,'*'）;title（'X1:
L=16'）;
subplot（4,1,3）;stem（n2,X2,'*'）;title（'X2:
L=32'）;
subplot（4,1,4）;stem（n3,X3,'*'）;title（'X3:
L=64'）;
1.基于FFT的卷积计算
对上述序列做基于FFT的卷积计算，做FFT时注意延拓序列长度，将结果显示出来，并和直接卷积结果比较.
clear;
N1=40;
N2=30;
N=41+31-1;
x1=cos（0.025*pi*（0:
N1）.^2）;
x=[x1zeros（1,N-N1）];
h1=sin（2.5*pi*（0:
N2））;
h=[h1zeros（1,N-N2）];
y1=conv（x1,h1）;
X=fft（x）;
H=fft（h）
Y=X.*H;
y=ifft（Y）;
subplot（2,1,1）;stem（y,'.'）;title（'基于FFT卷积计算'）;
subplot（2,1,2）;stem（y1,'.'）;title（直接卷积计算）;
2.基于FFT的相关计算
请计算h（n）和x（n）的基于FFT的相关函数，将结果显示出来，并和直接相关结果比较。
clear;
n1=0:
128;
n2=129:
192;
n3=193:
256;
n4=0:
31;
N=257+32-1;
x1=cos（0.25*pi*n1.^2）+sin（0.5*pi*n1）;
x2=cos（pi*n2）+sin（0.5*pi*n2）/6;
x3=sin（0.5*pi*n3）/6+（-1）.^n3;
x=[x1x2x3zeros（1,（N-257））];
h1=cos（pi*n4）+sin（0.5*pi*n4）/6;
h=[h1zeros（1,（N-32））];
r1=xcorr（[x1x2x3],h1）;
X=fft（x）;
H=fft（h）;
R=X.*conj（H）;
r=ifft（R）;
subplot（2,1,1）;stem（r,'.'）;title（'FFT相关计算'）;
subplot（2,1,2）;stem（r1,'.'）;title（'直接计算'）;
3.功率谱计算
信号设为：
f1=42Hz,f2=45Hz,f3=47Hz,采样频率fs=300Hz,采样点数N=256。
计算信号的功率谱，将功率谱图显示出来。
功率谱定义为P（k）=|X（k）|^2
clear;
f1=42;
f2=45;
f3=47;
fs=300;
N=64;
t=（0:
1:
N-1）/fs;
x=sin（2*pi*f1*t）+cos（2*pi*f2*t）+cos（2*pi*f3*t）;
X=abs（fft（x））;
P=X.^2;
stem（P,'.'）;title（功率谱）
1．给定差分方程y（n）-y（n-1）+0.9y（n-2）=x（n）
（1）计算并画出单位脉冲响应h（n），n=0,1,……,127；
（2）计算并画出单位阶跃响应s（n），n=0,1,……,127；
（3）判断系统是否稳定。
clear;
num=[100];
den=[1-10.9];
subplot（2,1,1）;dimpulse（num,den,128）;title（'单位冲激响应'）;%绘
subplot（2,1,2）;dstep（num,den,128）;title（'单位阶跃响应'）;
y1=dimpulse（num,den,128）%求
y2=dstep（num,den,128）
以下为定义法
%x1=zeros（1,130）;%单位冲激函数
%x1（3）=1;
%y1
（1）=0;y1
（2）=0;
%fork=1:
128
%y1（k+2）=x1（k+2）+y1（k+1）-0.9*y1（k）;
%end
%
%x2=ones（1,130）;%单位阶跃函数
%y2
（1）=0;y2
（2）=0;
%fork=1:
128
%y2（k+2）=x2（k+2）+y2（k+1）-0.9*y2（k）;
%end
%
%figure
%subplot（2,1,1）;plot（y1）;title（'单位脉冲响应'）
%subplot（2,1,2）;plot（y2）;title（'单位阶跃响应'）
%y1
%y2
2．给定系统函数
（1）画出H（z）零-极点示意图；
（2）求出并绘出H（z）的幅频响应、相频响应；
（3）求出并绘出该系统的单位抽样响应h（n）；
（4）求系统的单位阶跃响应s（n）。
clear;
num=[-0.20];
den=[100.7];
sys=tf（num,den）;
figure;pzmap（sys）;title（'零极点分布图'）
figure;freqz（num,den）;title（'频率响应'）
figure;
subplot（2,1,1）;dimpulse（num,den）;title（'单位冲激响应'）;
subplot（2,1,2）;dstep（num,den）;title（'单位阶跃响应'）;
y1=dimpulse（num,den,128）%求
y2=d
升级会员 