《数字信号处理》上机全部源代码调试通过完整版.docx
《《数字信号处理》上机全部源代码调试通过完整版.docx》由会员分享，可在线阅读，更多相关《《数字信号处理》上机全部源代码调试通过完整版.docx（24页珍藏版）》请在冰豆网上搜索。
《数字信号处理》上机全部源代码调试通过完整版
《数字信号处理》上机全部源代码调试通过，完整版
（高西全，第四版）
实验一
%实验1：
系统响应及系统稳定性
closeall;clearall
%调用fliter解差分方程，由系统对un的响应判断稳定性
%内容1：
调用filter解差分方程，由系统对u（n）的响应判断稳定性
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='y2（n）';tstem（y2n,y）;
title（'（c）系统对u（n）的响应y2（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='y2（n）';tstem（y2n,y）;
title（'（c）系统对u（n）的响应y2（n）'）
%内容2：
调用conv函数计算卷积
x1n=[11111111];%产生信号x1n=R8n
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）;
%调用函数tstem绘图
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）;%调用函数tstem绘图
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）;%产生信号un
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];
%系统差分方程系数向量B和A
y31n=filter（B,A,un）;%谐振器对un的响应y31n
y32n=filter（B,A,xsin）;
%谐振器对正弦信号的响应y32n
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）谐振器对正弦信号的响应y32（n）'）
functiontstem（xn,yn）
n=0:
length（xn）-1;
stem（n,xn,'.'）;
xlabel（'n'）;ylabel（'yn'）;
%xlabel（'n'）;ylabel（yn）;
axis（[0,n（end）,min（xn）,1.2*max（xn）]）;
实验二
%时域采样理论验证程序exp2a.m
Tp=64/1000;%观察时间Tp=64微秒
%产生M长采样序列x（n）
%Fs=1000;T=1/Fs;
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）;%M点FFT[xnt）]
yn='xa（nT）';subplot（3,2,1）;
tstem（xnt,yn）;%调用自编绘图函数tstem绘制序列图
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和Fs=200Hz的程序与上面Fs=1000Hz的程序完全相同。
%%%%%%%%%%%%fs=300Hz
Tp=64/1000;%观察时间Tp=64微秒
%产生M长采样序列x（n）
%Fs=1000;T=1/Fs;
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）;%M点FFT[xnt）]
yn='xa（nT）';figure
（2）;subplot（3,2,1）;
tstem（xnt,yn）;%调用自编绘图函数tstem绘制序列图
boxon;title（'（b）Fs=300Hz'）;
k=0:
M-1;fk=k/Tp;
subplot（3,2,2）;plot（fk,abs（Xk））;
title（'（b）T*FT[xa（nT）],Fs=1000Hz'）;
xlabel（'f（Hz）'）;ylabel（'幅度'）;
axis（[0,Fs,0,1.2*max（abs（Xk））]）
%%%%%%%%%%%%%%%%%%%fs=200Hz
Tp=64/1000;%观察时间Tp=64微秒
%产生M长采样序列x（n）
%Fs=1000;T=1/Fs;
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）;%M点FFT[xnt）]
yn='xa（nT）';figure（3）;subplot（3,2,1）;
tstem（xnt,yn）;%调用自编绘图函数tstem绘制序列图
boxon;title（'（c）Fs=200Hz'）;
k=0:
M-1;fk=k/Tp;
subplot（3,2,2）;plot（fk,abs（Xk））;
title（'（c）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;
%产生M长三角波序列x（n）
xa=0:
floor（M/2）;xb=ceil（M/2）-1:
-1:
0;xn=[xa,xb];
Xk=fft（xn,1024）;%1024点FFT[x（n）],用于近似序列x（n）的TF
X32k=fft（xn,32）;%32点FFT[x（n）]
x32n=ifft（X32k）;%32点IFFT[X32（k）]得到x32（n）
X16k=X32k（1:
2:
N）;%隔点抽取X32k得到X16（K）
x16n=ifft（X16k,N/2）;%16点IFFT[X16（k）]得到x16（n）
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,32,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]）
functiontstem（xn,yn）
n=0:
length（xn）-1;
stem（n,xn,'.'）;
xlabel（'n'）;ylabel（'yn'）;
%xlabel（'n'）;ylabel（yn）;
axis（[0,n（end）,min（xn）,1.2*max（xn）]）;
实验三
%实验三程序exp3.m
%用FFT对信号作频谱分析
%clearall;closeall
%实验内容
（1）================================
x1n=[ones（1,4）];%产生序列向量x1（n）=R4（n）
M=8;xa=1:
（M/2）;xb=（M/2）:
-1:
1;
x2n=[xa,xb];%产生长度为8的三角波序列x2（n）
x3n=[xb,xa];
X1k8=fft（x1n,8）;%计算x1n的8点DFT
X1k16=fft（x1n,16）;%计算x1n的16点DFT
X2k8=fft（x2n,8）;%计算x1n的8点DFT
X2k16=fft（x2n,16）;%计算x1n的16点DFT
X3k8=fft（x3n,8）;%计算x1n的8点DFT
X3k16=fft（x3n,16）;%计算x1n的16点DFT
%以下绘制幅频特性曲线
subplot（2,2,1）;mstem（X1k8）;
%绘制8点DFT的幅频特性图
title（'（1a）8点DFT[x_1（n）]'）;xlabel（'ω/π'）;
ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X1k8））]）
subplot（2,2,3）;mstem（X1k16）;
%绘制16点DFT的幅频特性图
title（'（1b）16点DFT[x_1（n）]'）;xlabel（'ω/π'）;
ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X1k16））]）
figure
（2）
subplot（2,2,1）;mstem（X2k8）;
%绘制8点DFT的幅频特性图
title（'（2a）8点DFT[x_2（n）]'）;xlabel（'ω/π'）;
ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X2k8））]）
subplot（2,2,2）;mstem（X2k16）;
%绘制16点DFT的幅频特性图
title（'（2b）16点DFT[x_2（n）]'）;xlabel（'ω/π'）;
ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X2k16））]）
subplot（2,2,3）;mstem（X3k8）;
%绘制8点DFT的幅频特性图
title（'（3a）8点DFT[x_3（n）]'）;xlabel（'ω/π'）;
ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X3k8））]）
subplot（2,2,4）;mstem（X3k16）;
%绘制16点DFT的幅频特性图
title（'（3b）16点DFT[x_3（n）]'）;xlabel（'ω/π'）;
ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X3k16））]）
%实验内容
（2）周期序列谱分析===================
N=8;n=0:
N-1;%FFT的变换区间N=8
x4n=cos（pi*n/4）;
x5n=cos（pi*n/4）+cos（pi*n/8）;
X4k8=fft（x4n）;%计算x4n的8点DFT
X5k8=fft（x5n）;%计算x5n的8点DFT
N=16;n=0:
N-1;%FFT的变换区间N=16
x4n=cos（pi*n/4）;
x5n=cos（pi*n/4）+cos（pi*n/8）;
X4k16=fft（x4n）;%计算x4n的16点DFT
X5k16=fft（x5n）;%计算x5n的16点DFT
figure（3）
subplot（2,2,1）;mstem（X4k8）;
%绘制8点DFT的幅频特性图
title（'（4a）8点DFT[x_4（n）]'）;
xlabel（'ω/π'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X4k8））]）
subplot（2,2,3）;mstem（X4k16）;
%绘制16点DFT的幅频特性图
title（'（4b）16点DFT[x_4（n）]'）;
xlabel（'ω/π'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X4k16））]）
subplot（2,2,2）;mstem（X5k8）;
%绘制8点DFT的幅频特性图
title（'（5a）8点DFT[x_5（n）]'）;xlabel（'ω/π'）;
ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X5k8））]）
subplot（2,2,4）;mstem（X5k16）;
%绘制16点DFT的幅频特性图
title（'（5b）16点DFT[x_5（n）]'）;xlabel（'ω/π'）;
ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X5k16））]）
%实验内容（3）模拟周期信号谱分析=================
figure（4）
Fs=64;T=1/Fs;
N=16;n=0:
N-1;%FFT的变换区间N=16
x6nT=cos（8*pi*n*T）+cos（16*pi*n*T）+cos（20*pi*n*T）;
%对x6（t）16点采样
X6k16=fft（x6nT）;%计算x6nT的16点DFT
X6k16=fftshift（X6k16）;%将零频率移到频谱中心
Tp=N*T;F=1/Tp;%频率分辨率F
k=-N/2:
N/2-1;fk=k*F;
%产生16点DFT对应的采样点频率（以零频率为中心）
subplot（3,1,1）;stem（fk,abs（X6k16）,'.'）;
boxon%绘制8点DFT的幅频特性图
title（'（6a）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;%FFT的变换区间N=16
x6nT=cos（8*pi*n*T）+cos（16*pi*n*T）+cos（20*pi*n*T）;
%对x6（t）32点采样
X6k32=fft（x6nT）;%计算x6nT的32点DFT
X6k32=fftshift（X6k32）;%将零频率移到频谱中心
Tp=N*T;F=1/Tp;%频率分辨率F
k=-N/2:
N/2-1;
fk=k*F;
%产生16点DFT对应的采样点频率（以零频率为中心）
subplot（3,1,2）;stem（fk,abs（X6k32）,'.'）;
boxon%绘制8点DFT的幅频特性图
title（'（6b）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;%FFT的变换区间N=16
x6nT=cos（8*pi*n*T）+cos（16*pi*n*T）+cos（20*pi*n*T）;
%对x6（t）64点采样
X6k64=fft（x6nT）;%计算x6nT的64点DFT
X6k64=fftshift（X6k64）;%将零频率移到频谱中心
Tp=N*T;F=1/Tp;%频率分辨率F
k=-N/2:
N/2-1;fk=k*F;
%产生16点DFT对应的采样点频率（以零频率为中心）
subplot（3,1,3）;stem（fk,abs（X6k64）,'.'）;
boxon%绘制8点DFT的幅频特性图
title（'（6a）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）
%mstem（Xk）绘制频域采样序列向量Xk的幅频特性图
M=length（Xk）;
k=0:
M-1;wk=2*k/M;%产生M点DFT对应的采样点频率（关于pi归一化值）
stem（wk,abs（Xk）,'.'）;boxon;%绘制M点DFT的幅频特性图?
xlabel（'w/\pi'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（Xk））]）;
实验四
%实验四程序exp4.m
%IIR数字滤波器设计及软件实现
%clearall;closeall
Fs=10000;T=1/Fs;%采样频率
%调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st
st=mstg;
%低通滤波器设计与实现=========================
fp=280;fs=450;
wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60;
%DF指标（低通滤波器的通、阻带边界频率）
[N,wp]=ellipord（wp,ws,rp,rs）;
%调用ellipord计算椭圆DF阶数N和通带截止频率wp
[B,A]=ellip（N,rp,rs,wp）;
%调用ellip计算椭圆带通DF系统函数系数向量B和A
y1t=filter（B,A,st）;%滤波器软件实现
%低通滤波器设计与实现绘图部分
figure
（2）;subplot（3,1,1）;
myplot（B,A）;
%调用绘图函数myplot绘制损耗函数曲线
yt='y_1（t）';
subplot（3,1,2）;tplot（y1t,T,yt）;
%调用绘图函数tplot绘制滤波器输出波形
%===========================================
%带通滤波器设计与实现=========================
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）;
%调用ellipord计算椭圆DF阶数N和通带截止频率wp
[B,A]=ellip（N,rp,rs,wp）;
%调用ellip计算椭圆带通DF系统函数系数向量B和A
y2t=filter（B,A,st）;%滤波器软件实现
%带通滤波器设计与实现绘图部分（省略）
%补充1
figure
（2）;subplot（3,1,3）;
myplot（B,A）;
%--
figure（3）;
yt='y_2（t）';
subplot（3,1,1）;tplot（y2t,T,yt）;
%===========================================
%高通滤波器设计与实现===========================
fp=890;fs=600;
wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60;
%DF指标（低通滤波器的通、阻带边界频率）
[N,wp]=ellipord（wp,ws,rp,rs）;
%调用ellipord计算椭圆DF阶数N和通带%截止频率wp
[B,A]=ellip（N,rp,rs,wp,'high'）;
%调用ellip计算椭圆带通DF系统函数系数向量B和A
y3t=filter（B,A,st）;%滤波器软件实现
%高低通滤波器设计与实现绘图部分（省略）
%补充2
figure（3）;subplot（3,1,2）;
myplot（B,A）;
%--
figure（3）;
yt='y_3（t）';
subplot（3,1,3）;tplot（y3t,T,yt）;
%===========================================
functionst=mstg
%产生信号序列向量st,并显示st的时域波形和频谱
%st=mstg返回三路调幅信号相加形成的混合信号,长度N=1600
N=160
升级会员 