《数字信号处理》上机全部源代码调试通过完整版.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