close all;clear all;clc
%内容1：调用filter解差分方程，由系统对u(n)的响应判断稳定性
A=[1,-0.9];B=[0.05,0.05]; %系统差分方程系数向量B和A x1n=[1 1 1 1 1 1 1 1 1 zeros(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=[1 1 1 1 1 1 1 1]; %产生信号x1n=R8n h1n=[ones(1,10) zeros(1,10)]; h2n=[1 2.5 2.5 1 zeros(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)')
function tstem(xn,yn) n=0:length(xn)-1; stem(n,xn,'.');box on 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);
box on;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);
box on;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);
box on;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);
box on;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,'.');box on
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),'.');box on
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,'.');box on 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),'.');box on
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,'.');box on 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),'.');box on
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),'.');box on
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),'.');box on
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))]) 实验三：
function st=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数字滤波器设计及软件实现 clear all;close all 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)]) 实验四：
function xt=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数字滤波器设计及软件实现 clear all;close all;
%===调用xtg产生信号xt,xt长度N=1000,并显示xt及其频谱===== N=1000;xt=xtg(N);
fp=120;fs=150;Rp=0.2;As=60;Fs=1000; %(1)用窗口函数法设计滤波器
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);
%以下为用窗函数法设计的绘图部分(滤波器损耗函数，滤波器输出信号波形)================= f=(0:1023)*Fs/1024; figure(2)
subplot(2,1,1)
plot(f,20*log10(Hw/max(Hw)));grid;title('(a)低通滤波器幅频特性') axis([0,Fs/2,-120,20]);
xlabel('f/Hz');ylabel('幅度'); N1=2000;
t=(0:N1-1)/Fs;Tp=N/Fs; subplot(2,1,2) plot(t,ywt);grid; length(t)
axis([0,Tp/2,-1,1]);xlabel('t/s');ylabel('y_w(t)'); title('(b)滤除噪声后的信号波形')
%(2)用等波纹最佳逼近法设计滤波器 fb=[fp,fs];m=[1,0];
dev=[(10^(Rp/20)-1)/(10^(Rp/20)+1),10^(-As/20)]; [Ne,f0,m0,W]=remezord(fb,m,dev,Fs); hn=remez(Ne,f0,m0,W); Hw=abs(fft(hn,1024)); yet=fftfilt(hn,xt,N);
%以下为用等波纹设计法的绘图部分========= figure(3);subplot(2,1,1); f=(0:1023)*Fs/1024;
plot(f,20*log10(Hw/max(Hw)));grid;title('(c)低通滤波器幅频特性') axis([0,Fs/2,-80,10]);
xlabel('f/Hz');ylabel('幅度'); subplot(2,1,2); plot(t,yet);grid;
axis([0,Tp/2,-1,1]);xlabel('t/s');ylabel('y_e(t)'); title('(d)滤除噪声后的信号波形') length(ywt) 实验五：
%DTMF双频拨号信号的生成和检测程序 %clear all;clc;
tm=[1,2,3,65;4,5,6,66;7,8,9,67;42,0,35,68]; %DTMF信号代表的16个数 N=205;K=[18,20,22,24,31,34,38,42]; f1=[697,770,852,941]; %行频率向量 f2=[1203,1336,1477,1633]; %列频率向量
TN=input('键入8位电话号码='); %输入六位数字
TNr=0; %接收端电话号码初值为0 for l=1:8;
d=fix(TN/10^(8-l)); TN=TN-d*10^(8-l); for p=1:4; for q=1:4;
if tm(p,q)==abs(d);break,end %检测码相符的列号q end
if tm(p,q)==abs(d);break,end %检测码相符的行号p end
n=0:1023; %为了发声，加长序列
x=sin(2*pi*n*f1(p)/8000)+sin(2*pi*n*f2(q)/8000); %构成双频信号 sound(x,8000); %发出声音 pause(0.1)
%接收检测端的程序
X=goertzel(x(1:205),K+1); %用Goertzel算法计算八点DFT样本 val=abs(X); %列出八点DFT向量
subplot(4,2,l);
stem(K,val,'.');grid;xlabel('k');ylabel('|X(k)|') %画出DFT(k)幅度 axis([10 50 0 120]) limit=80; for s=5:8;
if val(s)>limit,break,end %查找列号 end
for r=1:4;
if val(r)>limit,break,end %查找行号 end
TNr=TNr+tm(r,s-4)*10^(8-l); end
disp('接收端检测到的号码为：') %显示接收到的字符 disp(TNr) 实验一：
实验二：
实验三：
实验四：
实验五：
因篇幅问题不能全部显示，请点此查看更多更全内容