数字信号处理实验参考程序.docx
《数字信号处理实验参考程序.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验参考程序.docx（19页珍藏版）》请在冰豆网上搜索。
数字信号处理实验参考程序
实验一
一、程序清单
%内容1：
调用filter解差分方程，由系统对u（n）的响应判断系统稳定性
%=========================================================
closeall
clearall
A=[1,-0.9];B=[0.05,0.05];
x1n=[11111111zeros（1,50）];
x2n=ones（1,128）;
hn=impz（B,A,58）;
subplot（2,2,1）;y='h（n）';tstem（hn,y）;
title（'（a）系统单位脉冲响应h（n）'）
y1n=filter（B,A,x1n）;
subplot（2,2,2）;y='y1（n）';tstem（y1n,y）;
title（'（b）系统对R8（n）的响应y1（n）'）
y2n=filter（B,A,x2n）;
subplot（2,2,4）;y='y1（n）';tstem（y1n,y）;
title（'（c）系统对u（n）的响应y2（n）'）
%==========================================================
%内容2：
调用conv函数计算卷积
x1n=[11111111];
h1n=[ones（1,10）zeros（1,10）];
h2n=[12.52.51zeros（1,10）];
y21n=conv（h1n,x1n）;
y22n=conv（h2n,x1n）;
figure
（2）
subplot（2,2,1）;y='h1（n）';tstem（h1n,y）;
title（'（d）系统单位脉冲响应h1（n）'）
subplot（2,2,2）;y='y21（n）';tstem（y21n,y）;
title（'（e）h1（n）与R8（n）的卷积y21（n）'）
subplot（2,2,3）;y='h2（n）';tstem（h2n,y）;
title（'（f）系统对单位脉冲响应h2（n）'）
subplot（2,2,4）;y='y22（n）';tstem（y22n,y）;
title（'（g）h2（n）与R8（n）的卷积y22（n）'）
%=========================================================
%内容3：
谐振器分析
un=ones（1,256）;
n=0:
255;
xsin=sin（0.014*n）+sin（0.4*n）;
A=[1,-1.8237,0.9801];B=[1/100.49,0,-1/100.49];
y31n=filter（B,A,un）;
y32n=filter（B,A,xsin）;
figure（3）
subplot（2,1,1）;y='y31（n）';tstem（y31n,y）;
title（'（h）谐振器对u（n）的响应y31（n）'）
subplot（2,1,2）;y='y32（n）';tstem（y32n,y）;
title（'（i）谐振器对xsin的响应y32（n）'）
functiontstem（xn,yn）
n=0:
length（xn）-1;
stem（n,xn,’.’）;
xlabel（’n’）;ylabel（’yn’）;
axis（[0,n（end）,min（xn）,1.2*max（xn）]）;
二、实验结果与波形
三、分析与讨论
1、综合来说，在时域求系统响应的方法有两种：
一是通过解差分方程求得系统的输出，注意要合理选择初始条件；二是已知系统的单位脉冲响应，通过求输入信号和系统单位脉冲响应的线性卷积求得系统输出。
用计算机求解时最好使用MATLAB语言进行。
2、实际中要检验系统的稳定性，方法是在输入端加入单位阶跃序列，观察输出波形，如果波形稳定在一个常数值上，则系统稳定，否则不稳定。
上面第三个实验是稳定的。
3、谐振器具有对某个频率进行谐振的特性，本实验中的谐振器的谐振频率是0.4rad，因此稳定波形为sin（0.4n）。
4、如果输入信号为无限长序列，系统的单位脉冲响应是有限长序列，可以用分段线性卷积法求系统的响应，具体参考DFT一章内容。
如果信号经过低通滤波器，则信号的高频分量被滤掉，时域信号的变化减缓，在有阶跃处附近产生过渡带。
因此，当输入矩形序列时，输出序列的开始和终了都产生了明显的过渡带，见第一个实验结果的波形。
实验二
一、程序清单
%时域采样理论验证程序
%Fs=1000Hz
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））]）;
%====================================================================
%Fs=300Hz
Tp=64/1000;
Fs=300;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,3）;
tstem（xnt,yn）;
boxon;title（'（a）Fs=1000Hz'）;
k=0:
M-1;fk=k/Tp;
subplot（3,2,4）;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））]）;
%====================================================================
%Fs=200Hz
Tp=64/1000;
Fs=200;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,5）;
tstem（xnt,yn）;
boxon;title（'（a）Fs=1000Hz'）;
k=0:
M-1;fk=k/Tp;
subplot（3,2,6）;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））]）;
%频域采样理论验证程序
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）;
figure
（2）;
subplot（3,2,2）;stem（n,xn,'.'）;boxon
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）,'.'）;boxon
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,'.'）;boxon
title（'（d）16点IDFT[X_1_6（k）]'）;xlabel（'n'）;ylabel（'x_1_6（n）'）;axis（[0,16,0,20]）;
k=0:
N-1;
subplot（3,2,5）;stem（k,abs（X32k）,'.'）;boxon
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,'.'）;boxon
title（'（f）32点IDFT[X_3_2（k）]'）;xlabel（'n'）;ylabel（'x_3_2（n）'）;axis（[0,32,0,20]）;
二、实验结果与波形
三、分析讨论
1、时域采样理论验证，由图可见，当采样频率为1000Hz时，频谱混叠很小；当采样频率为30Hz时，频谱混叠很严重；当采样频率为200Hz时，频谱混叠更严重。
2、频域采样理论验证，对信号x（n）的频谱函数X（ejω）在【0，2Pi】上等间隔采样N=16时，N点IDFT[XN（k）]得到的序列正是原序列x（n）以16为周期进行周期延拓后的主值区序列：
xN（n）=IDFT[XN（k）]N=[∑x（n+iN）]RN（n），由于N实验三
一、程序清单及运行结果
%用FFT对信号进行谱分析
%内容1======================================================
x1n=[ones（1,4）];
M=8;xa=1:
（M/2）;xb=（M/2）:
-1:
1;x2n=[xa,xb];
x3n=[xb,xa];
X1k8=fft（x1n,8）;
X1k16=fft（x1n,16）;
X2k8=fft（x2n,8）;
X2k16=fft（x2n,16）;
X3k8=fft（x3n,8）;
X3k16=fft（x3n,16）;
subplot（3,2,1）;mstem（X1k8）;title（'（1a）8点DFT[x_1（n）]'）;
subplot（3,2,2）;mstem（X1k16）;title（'（1b）16点DFT[x_1（n）]'）
subplot（3,2,3）;mstem（X2k8）;title（'（2a）8点DFT[x_2（n）]'）;
subplot（3,2,4）;mstem（X2k16）;title（'（2b）16点DFT[x_2（n）]'）;
subplot（3,2,5）;mstem（X3k8）;title（'（3a）8点DFT[x_3（n）]'）;
subplot（3,2,6）;mstem（X3k16）;title（'（3b）16点DFT[x_3（n）]'）;
%用FFT对周期信号进行谱分析
%内容2======================================================
N=8;n=0:
N-1;
x4n=cos（pi*n/4）;
x5n=cos（pi*n/4）+cos（pi*n/8）;
X4k8=fft（x4n）;
X5k8=fft（x5n）;
N=16;n=0:
N-1;
x4n=cos（pi*n/4）;
x5n=cos（pi*n/4）+cos（pi*n/8）;
X4k16=fft（x4n）;
X5k16=fft（x5n）;
figure
（2）
subplot（2,2,1）;mstem（X4k8）;
title（'（1a）8点DFT[x_4（n）]'）;
subplot（2,2,2）;mstem（X4k16）;
title（'（1b）16点DFT[x_4（n）]'）
subplot（2,2,3）;mstem（X5k8）;
title（'（2a）8点DFT[x_5（n）]'）;
subplot（2,2,4）;mstem（X5k16）;
title（'（2b）16点DFT[x_5（n）]'）;
%模拟周期信号谱分析
%内容3======================================================
Fs=64;T=1/Fs;
N=16;n=0:
N-1;
x6nT=cos（8*pi*n*T）+cos（16*pi*n*T）+cos（20*pi*n*T）;
X6k16=fft（x6nT）;
X6k16=fftshift（X6k16）;
Tp=N*T;F=1/Tp;
k=-N/2:
N/2-1;fk=k*F;
figure（3）
subplot（3,1,1）;stem（fk,abs（X6k16）,'.'）;boxon
title（'（a）16点|DFT[x_6（nT）]|'）;xlabel（'f（Hz）'）;ylabel（'幅度'）;
axis（[-N*F/2-1,N*F/2-1,0,1.2*max（abs（X6k16））]）
N=32;n=0:
N-1;
x6nT=cos（8*pi*n*T）+cos（16*pi*n*T）+cos（20*pi*n*T）;
X6k32=fft（x6nT）;
X6k32=fftshift（X6k32）;
Tp=N*T;F=1/Tp;
k=-N/2:
N/2-1;fk=k*F;
subplot（3,1,2）;stem（fk,abs（X6k32）,'.'）;boxon
title（'（b）32点|DFT[x_6（nT）]|'）;xlabel（'f（Hz）'）;ylabel（'幅度'）;
axis（[-N*F/2-1,N*F/2-1,0,1.2*max（abs（X6k32））]）
N=64;n=0:
N-1;
x6nT=cos（8*pi*n*T）+cos（16*pi*n*T）+cos（20*pi*n*T）;
X6k64=fft（x6nT）;
X6k64=fftshift（X6k64）;
Tp=N*T;F=1/Tp;
k=-N/2:
N/2-1;fk=k*F;
subplot（3,1,3）;stem（fk,abs（X6k64）,'.'）;boxon
title（'（c）64点|DFT[x_6（nT）]|'）;xlabel（'f（Hz）'）;ylabel（'幅度'）;
axis（[-N*F/2-1,N*F/2-1,0,1.2*max（abs（X6k64））]）
functionmstem（Xk）
M=length（Xk）;
k=0:
M-1;wk=2*k/M;
stem（wk,abs（Xk）,'.'）;boxon
xlabel（'\omega/\pi'）;ylabel（'幅度'）;axis（[0,2,0,1.2*max（abs（Xk））]）
实验四
一、程序清单及运行结果
functionst=mstg
N=1600;
Fs=10000;T=1/Fs;Tp=N*T;
t=0:
T:
（N-1）*T;k=0:
N-1;f=k/Tp;
fc1=Fs/10;
fm1=fc1/10;
fc2=Fs/20;
fm2=fc2/10;
fc3=Fs/40;
fm3=fc3/10;
xt1=cos（2*pi*fm1*t）.*cos（2*pi*fc1*t）;
xt2=cos（2*pi*fm2*t）.*cos（2*pi*fc2*t）;
xt3=cos（2*pi*fm3*t）.*cos（2*pi*fc3*t）;
st=xt1+xt2+xt3;
fxt=fft（st,N）;
subplot（2,1,1）;plot（t,st）;grid;xlabel（'t/s'）;ylabel（'s（t）'）;
axis（[0,Tp/8,min（st）,max（st）]）;title（'（a）s（t）的波形'）
subplot（2,1,2）;stem（f,abs（fxt）/max（abs（fxt））,'.'）;grid;title（'（b）s（t）的频谱'）
axis（[0,Fs/5,0,1.2]）;xlabel（'f（Hz）'）;ylabel（'幅度'）;
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）;subplot（2,1,1）;yt='y_1t';myplot（B,A）;subplot（2,1,2）;tplot（y1t,T,yt）;
%带通滤波器设计与实现============================================
fp1=440;fpu=560;fs1=275;fsu=900;
wp=[2*fp1/Fs,2*fpu/Fs];ws=[2*fs1/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）;figure
（2）
subplot（2,1,1）;myplot（B,A）;yt='y_2t';subplot（2,1,2）;tplot（y2t,T,yt）;
%高通滤波器设计与实现============================================
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）;figure（3）
subplot（2,1,1）;myplot（B,A）;yt='y_3t';subplot（2,1,2）;tplot（y3t,T,yt）;
实验五
一、程序清单及运行结果
functionxt=xtg（N）
%采样频率Fs=1000Hz
%载波频率fc=Fs/10=100Hz,调制正弦波频率f0=fc/10=10Hz.
N=2000;Fs=1000;T=1/Fs;Tp=N*T;t=0:
T:
（N-1）*T;
fc=Fs/10;f0=fc/10;
mt=cos（2*pi*f0*t）;ct=cos（2*pi*fc*t）;xt=mt.*ct;
nt=2*rand（1,N）-1;
%=====设计高通滤波器hn,用于滤除噪声nt中的低频成分,生成高通噪声=====
fp=150;fs=200;Rp=0.1;As=70;
fb=[fp,fs];m=[0,1];
dev=[10^（-As/20）,（10^（Rp/20）-1）/（10^（Rp/20）+1）];
[n,fo,mo,W]=remezord（fb,m,dev,Fs）;
hn=remez（n,fo,mo,W）;
yt=filter（hn,1,10*nt）;
%===========================
xt=xt+yt;
fst=fft（xt,N）;k=0:
N-1;f=k/Tp;
subplot（2,1,1）;plot（t,xt）;grid;xlabel（'t/s'）;ylabel（'x（t）'）;
axis（[0,Tp/5,min（xt）,max（xt）]）;title（'（a）信号加噪声波形'）
subplot（2,1,2）;plot（f,abs（fst）/max（abs（fst）））;grid;title（'（b）信号加噪声的频谱'）
axis（[0,Fs/2,0,1.2]）;xlabel（'f/Hz'）;ylabel（'幅度'）
clearall;closeall;
N=1000;xt=xtg（N）;
fp=120;fs=150;Rp=0.2;As=60;Fs=1000;
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）;
figure
（2）
subplot（2,1,1）;
T=1/Fs;N=1024;k=0:
N-1;Tp=N*T;f=k/Tp;plot（f,20*log10（Hw））;
axis（[0,500,-110,10]）;gridon;
subplot（2,1,2）;tplot（ywt,1/Fs,ywt）;gridon;
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）;
figure（3）
subplot（2,1,1）;plot（f,20*log10（Hw））;axis（[0,500,-80,10]）;gridon;
subplot（2,1,2）;tplot（yet,1/Fs,yet）;gridon;
functiontplot（xn,T,yn）
n=0:
length（xn）-1;t=n*T;
plot（t,xn）;
xlabel（'t/s'）;ylabel（'yn'）;
axis（[0,0.5,min（xn）,1.2*max（xn）]）;
functionmyplot（B,A）
[H.W]=freqz（B,A,1000）;
m=abs（H）;
plot（W/pi,20
升级会员 