数字信号.docx
《数字信号.docx》由会员分享，可在线阅读，更多相关《数字信号.docx（23页珍藏版）》请在冰豆网上搜索。
数字信号
实验一：
closeall;clearall;clc
%内容1：
调用filter解差分方程，由系统对u（n）的响应判断稳定性
A=[1,-0.9];B=[0.05,0.05];%系统差分方程系数向量B和A
x1n=[111111111zeros（1,50）];%产生信号x1n=R8n
x2n=ones（1,128）;%产生信号x2n=un
hn=impz（B,A,58）;%求系统单位脉冲响应h（n）
y1n=filter（B,A,x1n）;%求系统对x1n的响应y1n
subplot（2,2,1）;y='h（n）';tstem（hn,y）;%调用函数tstem绘图
title（'（a）系统单位脉冲响应h（n）'）
subplot（2,2,2）;y='y1（n）';tstem（y1n,y）;
title（'（b）系统对R8（n）的响应y1（n）'）
y2n=filter（B,A,x2n）;%求系统对x2n的响应y2n
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,'.'）;boxon
xlabel（'n'）;ylabel（yn）;
axis（[0,n（end）,min（xn）,1.2*max（xn）]）axis（[0,n（end）,min（xn）,1.2*max（xn）]）
实验二：
（1）时域采样定理的验证程序
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）;
boxon;title（'（a）Fs=1000Hz'）;
k=0:
M-1;fk=k/Tp;
subplot（3,2,2）;plot（fk,abs（'（b）T*FT[xa（nT）],Fs=1000Hz'）;
xlabel（'f（Hz）'）;ylabel（'幅度'）；axis（[0,Fs,0,1.2*max（abs（Xk））]）
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）';subplot（3,2,3）;
tstem（xnt,yn）;
boxon;title（'（a）Fs=300Hz'）;
k=0:
M-1;fk=k/Tp;
subplot（3,2,4）;plot（fk,abs（'（d）T*FT[xa（nT）],Fs=300Hz'）;
xlabel（'f（Hz）'）;ylabel（'幅度'）；axis（[0,Fs,0,1.2*max（abs（Xk））]）
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）';subplot（3,2,3）;
tstem（xnt,yn）;
boxon;title（'（a）Fs=300Hz'）;
k=0:
M-1;fk=k/Tp;
subplot（3,2,4）;plot（fk,abs（'（d）T*FT[xa（nT）],Fs=300Hz'）;
xlabel（'f（Hz）'）;ylabel（'幅度'）；axis（[0,Fs,0,1.2*max（abs（Xk））]）
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）';subplot（3,2,5）;
tstem（xnt,yn）;
boxon;title（'（a）Fs=200Hz'）;
k=0:
M-1;fk=k/Tp;
subplot（3,2,6）;plot（fk,abs（'（f）T*FT[xa（nT）],Fs=200Hz'）;
xlabel（'f（Hz）'）;ylabel（'幅度'）；axis（[0,Fs,0,1.2*max（abs（Xk））]）
（2）频域采样定理的验证程序
%频域采样定理的验证程序
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[xn],用于近似序列xn的TF
X32k=fft（xn,32）;
x32n=ifft（X32k）;
X16k=X32k（1:
2:
N）;
x16n=ifft（X16k,N/2）;
k=0:
1023;wk=2*k/1024;
subplot（3,2,1）;plot（wk,abs（Xk））;title（'（a）FT[xn]'）;
xlabel（'\omega^pi'）;ylabel（'|X（e^j^\omega）|'）;axis（[0,1,0,200]）
subplot（3,2,2）;stem（n,xn,'.'）;boxon
title（'（b）三角波序列xn'）;xlabel（'n'）;ylabel（'xn'）;axis（[0,32,0,20]）
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）16点频域采样'）;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（'（d）32点IDFT[X_3_2（k）]'）;xlabel（'n'）;ylabel（'x_3_2（n）'）;axis（[0,32,0,20]）
（3）用FFT对信号作频谱分析程序
%用FFT对信号作频谱分析程序
%实验内容
（1）==============
x1n=ones（1,4）;
M=8;xa=1:
（M/2）;xb=（M/2）:
-1:
1;x2n=[xa,xb];%产生长度为8的三角波序列
x3n=[xb,xa];
N1=8;m=0:
N1-1;
x1=0:
7;w1=2*x1/N1;%FFT的变换区间N=8
N2=16;n=0:
N2-1;
x2=0:
15;w2=2*x2/N2;
X1k8=fft（x1n,8）;%计算x1n的8点DFT
X1k16=fft（x1n,16）;
X2k8=fft（x2n,8）;
X2k16=fft（x2n,16）;
X3k8=fft（x3n,8）;
X3k16=fft（x3n,16）;
%以下绘制幅频特性曲线
subplot（2,2,1）;stem（w1,abs（X1k8）,'.'）;%绘制8点DFT的幅频特性图
title（'（1a）8点DTF[x_1（n）]'）;xlabel（'w/pi'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X1k8））]）
subplot（2,2,3）;stem（w2,abs（X1k16）,'.'）;%绘制16点DFT的幅频特性图
title（'（1b）8点DTF[x_1（n）]'）;xlabel（'w/pi'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X1k16））]）
figure
（2）
subplot（2,2,1）;stem（w1,abs（X2k8）,'.'）;%绘制8点DFT的幅频特性图
title（'（2a）8点DTF[x_2（n）]'）;xlabel（'w/pi'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X2k8））]）
subplot（2,2,2）;stem（w2,abs（X2k16）,'.'）;%绘制16点DFT的幅频特性图
title（'（2b）16点DTF[x_2（n）]'）;xlabel（'w/pi'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X2k16））]）
subplot（2,2,3）;stem（w1,abs（X3k8）,'.'）;%绘制8点DFT的幅频特性图
title（'（3a）8点DTF[x_3（n）]'）;xlabel（'w/pi'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X3k8））]）
subplot（2,2,4）;stem（w2,abs（X3k16）,'.'）;%绘制16点DFT的幅频特性图
title（'（3b）16点DTF[x_3（n）]'）;xlabel（'w/pi'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X3k16））]）
%实验内容
（2）周期序列频谱分析==============
x4n=cos（pi*m/4）;
x5n=cos（pi*m/4）+cos（pi*m/8）;
X4k8=fft（x4n;%计算x4n的8点DTF
X5k8=fft（x5n）;
x4n=cos（pi*n/4）;
x5n=cos（pi*n/4）+cos（pi*n/8）;
X4k16=fft（x4n）;
X5k16=fft（x5n）;
figure（3）
subplot（2,2,1）;stem（w1,abs（X4k8）,'.'）;%绘制8点DFT的幅频特性图
title（'（4a）8点DTF[x_4（n）]'）;xlabel（'w/pi'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X4k8））]）
subplot（2,2,3）;stem（w2,abs（X4k16）,'.'）;%绘制16点DFT的幅频特性图
title（'（4b）16点DTF[x_4（n）]'）;xlabel（'w/pi'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X4k16））]）
subplot（2,2,2）;stem（w1,abs（X5k8）,'.'）;%绘制8点DFT的幅频特性图
title（'（5a）8点DTF[x_5（n）]'）;xlabel（'w/pi'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X5k8））]）
subplot（2,2,4）;stem（w2,abs（X5k16）,'.'）;%绘制16点DFT的幅频特性图
title（'（5b）8点DTF[x_5（n）]'）;xlabel（'w/pi'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X5k16））]）
figure（4）
%实验内容（3）模拟周期信号谱分析======================
Fs=64;T=1/Fs;
N=16;n=0:
N-1;%FFT的变换区间N=16
x6nT=cos（8*pi*n*T）+cos（116*pi*n*T）+cos（20*pi*n*T）;%对x6（t）16点采样
X6k16=fft（x6nT）;计算x6nT的16点DTF
X6k16=fftshift（X6k16）;将零频率移到频谱中心
Tp=N*T;F=1/Tp;
k=-N/2:
N/2-1;fk=k*F;%产生16点DTF对应的采样点频率（以零频率为中心）
subplot（3,1,1）;stem（fk,abs（X6k16）,'.'）;boxon
title（'（6a）16点|DTF[x_6（n）]|'）;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（116*pi*n*T）+cos（20*pi*n*T）;
X6k32=fft（x6nT）;计算x6nT的16点DTF
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（'（6b）32点|DTF[x_6（n）]|'）;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=64
x6nT=cos（8*pi*n*T）+cos（116*pi*n*T）+cos（20*pi*n*T）;
X6k64=fft（x6nT）;计算x6nT的16点DTF
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（'（6c）32点|DTF[x_6（n）]|'）;xlabel（'f（HZ）'）;ylabel（'幅度'）
axis（[-N*F/2-1,N*F/2-1,0,1.2*max（abs（X6k32））]）
实验三：
functionst=mstg
%st产生信号序列向量s（t）,并显示s（t）的时域波形和频谱
%st=mstg返回三路调幅信号相加形成的混合信号,长度N=1600
N=1600%N为信号st的长度
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）;
%===========以下为绘图部分，绘制st的时域波形和幅频特性曲线=========
subplot（3,1,1）
plot（t,st）;grid;xlabel（'t/s'）;ylabel（'s（t）'）;
axis（[0,Tp/2,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（'幅度'）
%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;
[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）;
%带通滤波器设计与实现=============
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（3）;subplot（3,1,1）;
myplot（B,A）;
yt='y_2（t）';
subplot（3,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（4）;subplot（3,1,1）;
myplot（B,A）;
yt='y_3（t）';
subplot（3,1,2）;tplot（y3t,T,yt）;
绘图函数tplot（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）]）
实验四：
functionxt=xtg（N）
%实验信号x（t）产生函数,并显示信号的幅频特性曲线
%xt=xtg产生一个长度为N,有加性高频噪声的单频调幅信号xt,采样频率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,f0,m0,W]=remezord（fb,m,dev,Fs）;
hn=remez（n,f0,m0,W）;
yt=filter（hn,1,10*nt）;
%===========================================
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）信号加噪声波形'）
axis（[0,Fs/2,0,1.2]）;xlabel（'f/Hz'）;ylabel（'幅度'）;
%FIR数字滤波器设计及软件实现
clearall;closeall;
%===调用xtg产生信号xt,xt长度N=1000,并显示xt及其频谱=====
N=1000;x
升级会员 