山东大学威海数字信号处理实验.docx
《山东大学威海数字信号处理实验.docx》由会员分享，可在线阅读，更多相关《山东大学威海数字信号处理实验.docx（22页珍藏版）》请在冰豆网上搜索。
山东大学威海数字信号处理实验
1.冲激信号
（1）x1[n]=0.9δ[n-5]1<=n<=20
程序：
nn=1:
20;
imp=zeros（20,1）;
imp（5）=1;x1=0.9.*imp;
stem（nn,x1）
ylabel（'x1[n]=0.9.*δ[n-5]'）
grid
gtext（'n'）
（2）x2[n]=0.8δ[n]-15<=n<=15
程序：
nn=-15:
15;
imp=zeros（31,1）;
imp（16）=1;x2=0.8.*imp;
stem（nn,x2）
ylabel（'x2[n]=0.8δ[n]'）
grid
gtext（'n'）
（3）x3[n]=1.5δ[n-333]300<=n<=350
程序：
nn=300:
350;
imp=zeros（51,1）;
imp（34）=1;x3=0.8.*imp;
stem（nn,x3）
ylabel（'x3[n]=1.5δ[n-333]'）
grid
gtext（'n'）
（4）x4[n]=4.5δ[n+7]-10<=n<=0
程序：
nn=-10:
0;
imp=zeros（11,1）;
imp（4）=1;x4=4.5.*imp;
stem（nn,x4）
ylabel（'x4[n]=4.5δ[n+7]'）
grid
gtext（'n'）
2、正弦信号
程序：
>>n=0:
25;
>>X1=sin（pi*n/17）;
>>stem（n,X1,'b'）
>>ylabel（'X1=sin（pi*n/17）'）
>>grid
>>gtext（'n'）
实验图示：
>>n=-15:
25;
>>X2=sin（pi*n/17）;
>>stem（n,X2,'b'）
>>ylabel（'X2=sin（pi*n/17）'）
>>grid
>>gtext（'n'）
实验图示：
>>n=-10:
10;
>>X3=sin（pi*3*n+pi/2）;
>>stem（n,X3,'b'）
>>ylabel（'X3=sin（pi*3*n+pi/2）'）
>>grid
>>gtext（'n'）
实验图示：
>>clear
>>n=0:
50;
>>x4=cos（pi*n/sqrt（23））;
>>stem（n,x4,'b'）
>>ylabel（'x4=cos（pi*n/sqrt（23））'）
>>grid
>>gtext（'n'）
实验图示：
3、指数信号
程序：
functiony=genexp（b,n0,L）
if（L<=0）
error（'GENEXP:
lengthnotpositive'）
end
nn=n0+[1:
L]'-1;
y=b.^nn;
end然后编写出x[n]=
x=genexp（0.9,1,20）;
stem（x）
程序：
nn=0:
20;
x=zeros（21,1）;
x
（1）=1;
a=[1,-0.9];
b=[1];
y=filter（b,a,x）;
stem（nn,y）
4、系统响应及系统的稳定性
A=[1-0.9];B=[0.050.05];
xn1=[ones（1,8）];
xn2=ones（1,30）;
xn3=[010];%单位冲击信号
yn1=filter（B,A,xn1）;
yn2=filter（B,A,xn2）;
yn3=filter（B,A,xn3）;%系统的单位脉冲响应
hn1=ones（1,10）;
hn2=[12.52.51];
y1n=conv（xn1,hn1）;%用系统响应求响应
y2n=conv（xn1,hn2）;
AA=[1-1.82370.9801];b=1/100.49;
BB=[b-b];
xn3=ones（1,120）;
y3n=filter（BB,AA,xn3）;%第四问的系统输出响应
n=0:
30;
xn=sin（0.014*n）+sin（0.4*n）;
ynn=filter（BB,AA,xn）;%xn的系统响应
subplot（3,3,1）;
m=0:
length（yn1）-1;
stem（m,yn1,'.'）;
title（'当输入信号为宽度为8的矩形信号时'）;
subplot（3,3,2）;
m=0:
length（yn2）-1;
stem（m,yn2,'.'）;
title（'当输入信号为阶跃序列信号时'）;
subplot（3,3,3）;
m=0:
length（yn3）-1;
stem（m,yn3,'.'）;
title（'当输入信号为单位冲击信号时'）;
subplot（3,3,4）;
m=0:
length（y1n）-1;
stem（m,y1n,'.'）;
title（'利用系统的单位脉冲响应hn1求xn1的响应'）;
subplot（3,3,5）;
m=0:
length（y2n）-1;
stem（m,y2n,'.'）;
title（'利用系统的单位脉冲响应hn2求xn1的响应'）;
subplot（3,3,6）;
m=0:
length（y3n）-1;
stem（m,y3n,'.'）;
title（'当输入信号为阶跃序列时求谐振器的响应'）;
subplot（3,3,7）;
m=0:
length（ynn）-1;
stem（m,ynn,'.'）;
title（'当输入信号为xn时求谐振器的响应’）
5、时域采样与频域采样
%时域采样理论验证程序exp2a.m
Tp=64/1000;
Fs=1000;T=1/Fs;
M=Tp*Fs;n=0:
M-1;
A=444.128;alph=pi*50*2^0.5;omega=pi*50*2^0.5;
xnt=A*exp（-alph*n*T）.*sin（omega*n*T）;
Xk=T*fft（xnt,M）;
yn='xa（nT）';subplot（3,2,1）;
tstem（xnt,yn）;
boxon;title（'（a）Fs=1000Hz'）;
k=0:
M-1;fk=k/Tp;
subplot（3,2,2）;plot（fk,abs（Xk））;title（'（a）T*FT[xa（nT）],Fs=1000Hz'）;
xlabel（'f（Hz）'）;ylabel（'幅度'）;axis（[0,Fs,0,1.2*max（abs（Xk））]）
%频域采样理论验证程序exp2b.m
M=27;N=32;n=0:
M;
xa=0:
floor（M/2）;xb=ceil（M/2）-1:
-1:
0;xn=[xa,xb];
Xk=fft（xn,1024）;
X32k=fft（xn,32）;
x32n=ifft（X32k）;
X16k=X32k（1:
2:
N）;
x16n=ifft（X16k,N/2）;
subplot（3,2,2）;stem（n,xn,'.'）;
title（'（b）三角波序列x（n）'）;xlabel（'n'）;ylabel（'x（n）'）;axis（[0,32,0,20]）
k=0:
1023;wk=2*k/1024;
subplot（3,2,1）;plot（wk,abs（Xk））;title（'（a）FT[x（n）]'）;
xlabel（'\omega/\pi'）;ylabel（'|X（e^j^\omega）|'）;
axis（[0,1,0,200]）
k=0:
N/2-1;
subplot（3,2,3）;stem（k,abs（X16k）,'.'）;
title（'（c）16点频域采样'）;xlabel（'k'）;ylabel（'|X_1_6（k）|'）;
axis（[0,8,0,200]）
n1=0:
N/2-1;
subplot（3,2,4）;stem（n1,x16n,'.'）;
title（'（d）16点IDFT[X_1_6（k）]'）;xlabel（'n'）;ylabel（'x_1_6（n）'）;
axis（[0,32,0,20]）
k=0:
N-1;
subplot（3,2,5）;stem（k,abs（X32k）,'.'）;
title（'（e）32点频域采样'）;xlabel（'k'）;ylabel（'|X_3_2（k）|'）;axis（[0,16,0,200]）
n1=0:
N-1;
subplot（3,2,6）;stem（n1,x32n,'.'）;
title（'（f）32点IDFT[X_3_2（k）]'）;xlabel（'n'）;ylabel（'x_3_2（n）'）;axis（[0,32,0,20]）
2时域采样理论的验证程序exp2b.m运行结果如图10.3.3所示。
图10.3.3
6、自相关与互相关函数
clc;
k1=3;
k2=3;
k=k1+k2-1;
f1=[1,1,1];
f2=[0,1,2,3];
f=conv（f1,f2）;
nf=0:
k;
stem（nf,f,'*r'）;xlabel（'n'）;
ylabel（'f（n）'）;
gridon
clear;
N=500;
p1=1;
p2=0.1;
f=1/8;
Mlag=50;
u=randn（1,N）;
n=[0:
N-1];
s=sin（2*pi*f*n）;
u1=u*sqrt（p1）;
x1=u1（1:
N）+s;
rx1=xcorr（x1,Mlag,'biased'）;
subplot（221）;
plot（x1（1:
Mlag））;
xlabel（'n'）;
ylabel（'x1（n）'）;gridon;
subplot（223）;
plot（（-Mlag:
Mlag）,rx1）;gridon;
xlabel（'m'）;ylabel（'rx1（m）'）;
u2=u*sqrt（p2）;
x2=u2（1:
N）+s;
rx2=xcorr（x2,Mlag,'biased'）;
subplot（222）;
plot（x2（1:
Mlag））;
xlabel（'n'）;ylabel（'x2（n）'）;gridon;
subplot（224）;
plot（（-Mlag:
Mlag）,rx2）;
gridon;xlabel（'m'）;ylabel（'rx2（m）'）;
7、离散傅立叶变换
clearall;closeall;clc;
n=[0:
1:
99];x=cos（0.48*pi*n）+cos（0.52*pi*n）;
n1=[0:
1:
9];y1=x（1:
1:
10）;
subplot（2,2,1）;stem（n1,y1）;title（'信号x（n），0<=n<=9'）;
axis（[0,10,-2.5,2.5]）;text（10.2,-2.5,'n'）
Y1=dft（y1,10）;magY1=abs（Y1（1:
1:
6））;
k1=0:
1:
5;w1=2*pi/10*k1;
subplot（2,2,2）;plot（w1/pi,magY1）;title（'DTFT幅度的样本）;
xlabel（'频率（单位pi）'）
n=[0:
1:
99];x=cos（0.48*pi*n）+cos（0.52*pi*n）;
n2=[0:
1:
99];y2=[x（1:
1:
10）,zeros（1,90）];
subplot（2,2,3）;stem（n2,y2）;title（'信号x（n），0<=n<=99'）;
Y2=dft（y2,100）;magY2=abs（Y2（1:
1:
51））;
k2=0:
1:
50;w2=2*pi/100*k2;
subplot（2,2,4）;plot（w2/pi,magY2）;title（'DTFT幅度的样本'）;
axis（[0,1,0,10]）
xlabel（'频率（单位pi）'）
clearall;closeall;clc;
n=[0:
1:
99];x=cos（0.48*pi*n）+cos（0.52*pi*n）;
subplot（2,1,1）;stem（n,x）;title（'信号x（n），0<=n<=99'）;xlabel（'n'）
X=dft（x,100）;magX=as（X（1:
1:
51））;
k=0:
1:
50;w=2*pi/100*k;
subplot（2,1,2）;plot（w/pi,magX）;title（'DTFT幅度'）;
xlabel（'频率（单位pi）'）
结果：
图一
图二
8、IIR数字滤波器设计及软件实现
clearall;closeall
Fs=10000;T=1/Fs;
st=mstg;
fp=280;fs=450;
wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60;
[N,wp]=ellipord（wp,ws,rp,rs）;
[B,A]=ellip（N,rp,rs,wp）;
y1t=filter（B,A,st）;
figure
（2）;subplot（3,1,1）;
myplot（B,A）;
yt='y_1（t）';
subplot（3,1,2）;tplot（y1t,T,yt）;
fpl=440;fpu=560;fsl=275;fsu=900;
wp=[2*fpl/Fs,2*fpu/Fs];ws=[2*fsl/Fs,2*fsu/Fs];rp=0.1;rs=60;
[N,wp]=ellipord（wp,ws,rp,rs）;
[B,A]=ellip（N,rp,rs,wp）;
y2t=filter（B,A,st）;
fp=890;fs=600;
wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60;
[N,wp]=ellipord（wp,ws,rp,rs）;
[B,A]=ellip（N,rp,rs,wp,'high'）;
y3t=filter（B,A,st）;
（a）低通滤波器损耗函数及其分离出的调幅信号y1（t）
（b）带通滤波器损耗函数及其分离出的调幅信号y2（t）
（c）高通滤波器损耗函数及其分离出的调幅信号y3（t）
9、FIR数字滤波器设计与软件实现
clearall;closeall;
N=1000;xt=xtg（N）;
fp=120;fs=150;Rp=0.2;As=60;Fs=1000;
%
（1）用窗函数法设计滤波器
wc=（fp+fs）/Fs;
B=2*pi*（fs-fp）/Fs;
Nb=ceil（11*pi/B）;
hn=fir1（Nb-1,wc,blackman（Nb））;
Hw=abs（fft（hn,1024））;
ywt=fftfilt（hn,xt,N）;
fb=[fp,fs];m=[1,0];
dev=[（10^（Rp/20）-1）/（10^（Rp/20）+1）,10^（-As/20）];
[Ne,fo,mo,W]=remezord（fb,m,dev,Fs）;
hn=remez（Ne,fo,mo,W）;
Hw=abs（fft（hn,1024））;
yet=fftfilt（hn,xt,N）;
升级会员 