《数字信号处理第三版》第10章实验程序验证版西安电子科技大学出版 主编高西全 丁玉美.docx
《《数字信号处理第三版》第10章实验程序验证版西安电子科技大学出版 主编高西全 丁玉美.docx》由会员分享，可在线阅读，更多相关《《数字信号处理第三版》第10章实验程序验证版西安电子科技大学出版 主编高西全 丁玉美.docx（14页珍藏版）》请在冰豆网上搜索。
《数字信号处理第三版》第10章实验程序验证版西安电子科技大学出版主编高西全丁玉美
《数字信号处理（第三版）》第10章实验程序
注意：
每次实验前，请新建以下8个m文件并依实验需要（已注明）移动至matlab当前工作目录下，便于主程序调用。
各程序均已验证，请直接复制即可。
文件名：
tstem.m（实验一、二需要）
程序：
functiontstem（xn,yn）
%时域序列绘图函数
%xn:
被绘图的信号数据序列，yn:
绘图信号的纵坐标名称（字符串）
n=0:
length（xn）-1;
stem（n,xn,'.'）;
xlabel（'n'）;ylabel（'yn'）;
axis（[0,n（end）,min（xn）,1.2*max（xn）]）;
文件名：
tplot.m（实验一、四需要）
程序：
functiontplot（xn,T,yn）
%时域序列连续曲线绘图函数
%xn:
信号数据序列,yn:
绘图信号的纵坐标名称（字符串）
%T为采样间隔
n=0;length（xn）-1;t=n*T;
plot（t,xn）;
xlabel（'t/s'）;ylabel（yn）;
axis（[0,t（end）,min（xn）,1.2*max（xn）]）;
文件名：
myplot.m（实验一、四需要）
程序：
%
（1）myplot;计算时域离散系统损耗函数并绘制曲线图。
functionmyplot（B,A）
%B为系统函数分子多项式系数向量
%A为系统函数分母多项式系数向量
[H,W]=freqz（B,A,1000）
m=abs（H）;
plot（W/pi,20*log10（m/max（m）））;gridon;
xlabel（'\omega/\pi'）;ylabel（'幅度（dB）'）
axis（[0,1,-80,5]）;title（'损耗函数曲线'）;
文件名：
mstem.m（实验一、三需要）
程序：
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））]）;
文件名：
mpplot.m（实验一需要）
程序：
%
（2）mpplot;计算时域离散系统损耗函数和相频特性函数，并绘制曲线图。
functionmpplot（B,A,Rs）
%mpplot（B,A,Rs）
%时域离散系统损耗函数和相频特性绘图
%B为系统函数分子多项式系数向量
%A为系统函数分母多项式系数向量
%Rs为滤波器阻带最小衰减，省略则幅频曲线最小值取－80dB
ifnargin<3ymin=-80;elseymin=-Rs-20;end;%确定幅频曲线纵坐标最小值
[H,W]=[H,W]=freqz（B,A,1000）
m=abs（H）;
subplot（2,2,1）;
plot（W/pi,20*log10（m/max（m）））;gridon;
xlabel（'\omega/\pi'）;ylabel（'幅度（dB）'）
axis（[0,1,ymin,5]）;title（'损耗函数曲线'）;
subplot（2,2,3）;
plot（W/pi,p/pi）;
xlabel（'\omega/\pi'）;ylabel（'相位/\pi'）;gridon;
title（'（b）相频特性曲线'）;
文件名：
mfftplot.m（实验一需要）
程序：
functionmfftplot（xn,N）
%mfftplot（xn,N）计算序列向量xn的N点fft并绘制其幅频特性曲线
Xk=fft（xn,N）;%计算信号xn的频谱的N点采样
%===以下为绘图部分====
k=0:
N-1;wk=2*k/N;
m=abs（Xk）;mm=max（m）;
plot（wk,m/mm）;gridon;
xlabel（'\omega/\pi'）;ylabel（'幅度（dB）'）;
axis（[0,2,0,1.2]）;
title（'幅度特性曲线'）;
文件名：
mstg.m（实验四需要）
程序：
functionst=mstg
%产生信号序列向量st,并显示st的时域波形和频谱
%st=mstg返回三路调幅信号相加形成的混合信号，长度N=1600
N=1600;%N为信号st的长度
Fs=10000;T=1/Fs;Tp=N*T;%采样频率Fs=10kHz，Tp为采样时间
t=0:
T:
（N-1）*T;k=0:
N-1;f=k/Tp;
fc1=Fs/10;%第1路调幅信号的载波频率fc1=1000Hz
fm1=fc1/10;%第1路调幅信号的调制信号频率fm1=100Hz
fc2=Fs/20;%第2路调幅信号的载波频率fc2=500Hz
fm2=fc2/10;%第2路调幅信号的调制信号频率fm2=50Hz
fc3=Fs/40;%第3路调幅信号的载波频率fc3=250Hz
fm3=fc3/10;%第3路调幅信号的调制信号频率fm3=25Hz
xt1=cos（2*pi*fm1*t）.*cos（2*pi*fc1*t）;%产生第1路调幅信号
xt2=cos（2*pi*fm2*t）.*cos（2*pi*fc2*t）;%产生第2路调幅信号
xt3=cos（2*pi*fm3*t）.*cos（2*pi*fc3*t）;%产生第3路调幅信号
st=xt1+xt2+xt3;%三路调幅信号相加
fxt=fft（st,N）;%计算信号st的频谱
%====以下为绘图部分，绘制st的时域波形和幅频特性曲线====
subplot（3,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（3,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（'幅度'）;
文件名：
xtg.m（实验五需要）
程序：
functionxt=xtg（N）
%实验五信号x（t）产生,并显示信号的幅频特性曲线
%xt=xtg（N）产生一个长度为N,有加性高频噪声的单频调幅信号xt,采样频率Fs=1000Hz
%载波频率fc=Fs/10=100Hz,调制正弦波频率f0=fc/10=10Hz.
Fs=1000;T=1/Fs;Tp=N*T;
t=0:
T:
（N-1）*T;
fc=Fs/10;f0=fc/10;%载波频率fc=Fs/10，单频调制信号频率为f0=Fc/10;
mt=cos（2*pi*f0*t）; %产生单频正弦波调制信号mt，频率为f0
ct=cos（2*pi*fc*t）; %产生载波正弦波信号ct，频率为fc
xt=mt.*ct; %相乘产生单频调制信号xt
nt=2*rand（1,N）-1; %产生随机噪声nt
%设计高通滤波器hn,用于滤除噪声nt中的低频成分,生成高通噪声
%========================
fp=150;fs=200;Rp=0.1;As=70; %滤波器指标
fb=[fp,fs];m=[0,1]; %计算remezord函数所需参数f,m,dev
dev=[10^（-As/20）,（10^（Rp/20）-1）/（10^（Rp/20）+1）];
[n,fo,mo,W]=remezord（fb,m,dev,Fs）; %确定remez函数所需参数
hn=remez（n,fo,mo,W）; %调用remez函数进行设计,用于滤除噪声nt中的低频成分
yt=filter（hn,1,10*nt）; %滤除随机噪声中低频成分，生成高通噪声yt
xt=xt+yt; %噪声加信号
fst=fft（xt,N）;k=0:
N-1;f=k/Tp;
subplot（3,1,1）;plot（t,xt）;grid;xlabel（'t/s'）;ylabel（'x（t）'）;
axis（[0,Tp/5,min（xt）,max（xt）]）;title（'（a）信号加噪声波形'）
subplot（3,1,2）;plot（f,abs（fst）/max（abs（fst）））;grid;title（'（b）信号加噪声的频谱'）
axis（[0,Fs/2,0,1.2]）;xlabel（'f/Hz'）;ylabel（'幅度'）
10.1系统响应及系统稳定性
closeall;clearall;clc;
%内容1：
调用filter解差分方程，由系统对u（n）的响应判断稳定性
%========================
A=[1,-0.9];B=[0.05,0.05]; %系统差分方程系数向量B和A
x1n=[11111111zeros（1,50）]; %产生信号x1（n）=R8（n）
x2n=ones（1,128）; %产生信号x2（n）=u（n）
hn=impz（B,A,58）; %求系统单位脉冲响应h（n）
subplot（2,2,1）;y='h（n）';tstem（hn,y）; %调用函数tstem绘图
title（'（a）系统单位脉冲响应h（n）'）;boxon
y1n=filter（B,A,x1n）; %求系统对x1（n）的响应y1（n）
subplot（2,2,2）;y='y1（n）';tstem（y1n,y）;
title（'（b）系统对R8（n）的响应y1（n）'）;boxon
y2n=filter（B,A,x2n）; %求系统对x2（n）的响应y2（n）
subplot（2,2,4）;y='y2（n）';tstem（y2n,y）;
title（'（c）系统对u（n）的响应y2（n）'）;boxon
%内容2：
调用conv函数计算卷积
%========================
x1n=[11111111]; %产生信号x1（n）=R8（n）
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）'）;boxon
subplot（2,2,2）;y='y21（n）';tstem（y21n,y）;
title（'（e）h1（n）与R8（n）的卷积y21（n）'）;boxon
subplot（2,2,3）;y='h2（n）';tstem（h2n,y）; %调用函数tstem绘图
title（'（f）系统单位脉冲响应h2（n）'）;boxon
subplot（2,2,4）;y='y22（n）';tstem（y22n,y）;
title（'（g）h2（n）与R8（n）的卷积y22（n）'）;boxon
%内容3：
谐振器分析
%========================
un=ones（1,256）; %产生信号u（n）
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）; %谐振器对u（n）的响应y31（n）
y32n=filter（B,A,xsin）; %谐振器对u（n）的响应y31（n）
figure（3）
subplot（2,1,1）;y='y31（n）';tstem（y31n,y）;
title（'（h）谐振器对u（n）的响应y31（n）'）;boxon
subplot（2,1,2）;y='y32（n）';tstem（y32n,y）;
title（'（i）谐振器对正弦信号的响应y32（n）'）;boxon
10.2时域采样与频域采样
时域采样理论验证程序exp2a.m：
%========================
closeall;clearall;clc;
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完全相同。
频域采样理论验证程序exp2b.m：
%========================
closeall;clearall;clc;
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]）
10.3用FFT对信号做频谱分析
exp301.m：
%========================
clearall;closeall;clc;
%实验内容
（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））]）
exp302.m：
%实验内容
（2）周期序列谱分析
%========================
closeall;clearall;clc;
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））]）
exp303.m：
%实验内容（3）模拟周期信号谱分析
%========================
closeall;clearall;clc;
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点
升级会员 