数字信号处理第三版高西泉上机实验完整版（含子程序）- 副本

实验一：系统响应及系统稳定性

一、 实验原理与方法

1、在时域求系统响应的方法有两种，第一种是通过解差分方程求得系统输出； 第二种是已知系统的单位脉冲响应，通过求输入信号和系统单位脉冲响应的线性卷积求得系统输出。

2、检验系统的稳定性，其方法是在输入端加入单位阶跃序列，观察输出波形，如果波形稳定在一个常数值上，系统稳定，否则不稳定。

二、 实验结果及程序清单（含分析及函数用法）

1、%-----（2）调用filter解差分方程以及单位脉冲响应------
close all;clear all
A=[1,-0.9];B=[0.05,0.05]; %系统差分方程系数向量B和A
x1n=[ones(1,8) zeros(1,25)]; %产生信号x1(n)=R8(n)，用zeros用来加点的个数
x2n=ones(1,30); %产生信号x2(n)=u(n)
hn=impz(B,A,25); %求系统单位脉冲响应h(n)
subplot(3,1,1);stem(hn); %调用函数stem绘图
title('(a) 系统单位脉冲响应h(n)');
y1n=filter(B,A,x1n); %求系统对x1(n)的响应y1(n)
subplot(3,1,2);stem(y1n);
title('(b) 系统对R8(n)的响应y1(n)');
y2n=filter(B,A,x2n); %求系统对x2(n)的响应y2(n)
subplot(3,1,3);stem(y2n);
title('(c) 系统对u(n)的响应y2(n)');

分析：

（a）25个点数和程序所写一致。

Filter函数实现线性常系数差分方程的递推求解。

调用格式如下：

Y=[filter(B,A,x)] 计算系统对输入信号x的零状态响应输出信号向量Y，

B,A是差分方程的系数向量。即

B=[a1,a2……am] A=[b1,b2……bn]

2、%-----（3）调用conv函数计算卷积-------
x1n=[ones(1,8)]; %产生信号x1(n)=R8(n)
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);
subplot(2,2,1);stem(h1n);
title('(d) 系统单位脉冲响应h1(n)');
subplot(2,2,3);stem(y21n);
title('(e) h1(n)与R8(n)的卷积y21(n)');
subplot(2,2,2);stem(h2n);
title('(f) 系统单位脉冲响应h2(n)');
subplot(2,2,4);stem(y22n);
title('(g) h2(n)与R8(n)的卷积y22(n)');

分析：

（d）（f）单位脉冲响应点数与程序要求一致

（e）（g）卷积点数满足M+N-1的要求，图形也满足要求。

Conv函数用于计算两个有限长序列的卷积

C=conv（A,B）计算两个有限长序列向量A和B的卷积

3、%-----（4）实验方法检查系统是否稳定-------
close all;clear all
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
y1n=filter(B,A,un); %谐振器对u(n)的响应y31(n)
y2n=filter(B,A,xsin); %谐振器对u(n)的响应y31(n)
subplot(2,1,1);stem(y1n);
title('(h) 谐振器对u(n)的响应y31(n)');
subplot(2,1,2);stem(y2n);
title('(i) 谐振器对正弦信号的响应y32(n)');

分析：

（h）中由：在系统的输入端加入单位阶跃序列，如果系统的输出趋近一个常数（包括零），就可以断定系统是稳定的

输出显然趋近于零，所以是稳定的

（i）中谐振器具有对某个频率进行谐振的性质，本实验中的谐振器的谐振频率是0.4 rad,因此稳定波形为sin(0.4n)。

三、 思考题简要分析

(1) 如果输入信号为无限长序列，系统的单位脉冲响应是有限长序列，可否用线性卷积法求系统的响应? 如何求？

答： 如果输入信号为无限长序列，系统的单位脉冲响应是有限长序列，可用分段线性卷积法求系统的响应

（2）如果信号经过低通滤波器，把信号的高频分量滤掉，时域信号会有何变化，用前面 第一个实验结果进行分析说明。

答：如果信号经过低通滤波器， 则信号的高频分量被滤掉， 时域信号的变化减缓，在有阶跃处附近产生过渡带。因此，当输入矩形序列时，输出序列的开始和终了都产生了明显的过渡带，见第一个实验结果的波形。

四、 总结

总结即在实验原理中说明的两点：

1、在时域求系统响应的方法有两种，第一种是通过解差分方程求得系统输出； 第二种是已知系统的单位脉冲响应，通过求输入信号和系统单位脉冲响应的线性卷积求得系统输出。

2、检验系统的稳定性，其方法是在输入端加入单位阶跃序列，观察输出波形，如果波形稳定在一个常数值上，系统稳定，否则不稳定。

实验二：时域采样与频域采样

一、 实验原理与方法

1、 时域采样定理：

a) 对模拟信号以间隔T进行时域等间隔理想采样，形成的采样信号的频谱是原模拟信号频谱以采样角频率（）为周期进行周期延拓。公式为：

b) 采样频率必须大于等于模拟信号最高频率的两倍以上，才能使采样信号的

频谱不产生频谱混叠。

2、 频域采样定理：

公式为：

由公式可知，频域采样点数N必须大于等于时域离散信号的长度M(即N≥M)，才能使时域不产生混叠，则N点IDFT[]得到的序列就是原序列x(n),即=x(n)。

二、 实验结果及程序清单（含分析及函数用法）

1、% ---------（1）时域采样理论验证程序---------
Tp=64/1000; %观察时间Tp=64毫秒
%产生M长采样序列x(n)
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);
stem(xnt); %调用自编绘图函数stem绘制序列图
box on;title('(a) Fs=1000Hz');
k=0:M-1;fk=k/Tp;
subplot(3,2,2);stem(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完全相同。

Tp=64/1000; %观察时间Tp=64毫秒
%产生M长采样序列x(n)
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);
stem(xnt); %调用自编绘图函数stem绘制序列图
box on;title('(b) Fs=300Hz');
k=0:M-1;fk=k/Tp;
subplot(3,2,4);stem(fk,abs(Xk));title('(b) T*FT[xa(nT)],Fs=300Hz');
xlabel('f(Hz)');ylabel('幅度');axis([0,Fs,0,1.2*max(abs(Xk))]);

Tp=64/1000; %观察时间Tp=64毫秒
%产生M长采样序列x(n)
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);
stem(xnt); %调用自编绘图函数stem绘制序列图
box on;title('(c) Fs=200Hz');
k=0:M-1;fk=k/Tp;
subplot(3,2,6);stem(fk,abs(Xk));title('(c) T*FT[xa(nT)],Fs=200Hz');
xlabel('f(Hz)');ylabel('幅度');axis([0,Fs,0,1.2*max(abs(Xk))]);

分析：，当采样频率为1000 Hz时，频谱混叠很小：当采样频率为300 Hz时，频谱混叠很严重；当采样频率为200 Hz时，频谱混叠更很严重。

Fft函数的调用格式：

Xk=fft(xn,N)

调用参数xn为被交换的时域序列向量，N是DFT变换的区间长度，当N大于xn的长度时，fft函数自动在xn后面补零。当N小于xn的长度时，fft函数计算xn的前面N个元素构成的N长序列的N点DFT,忽略xn后面的元素。

2、%-------------(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[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);box on
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));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) 32点频域采样');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('(f)32点IDFT[X_3_2(k)]');

xlabel('n');ylabel('x_3_2(n)');axis([0,32,0,20])

分析;该图验证了频域采样理论和频域采样定理。 对信号x(n)的频谱函数X(ejω)在［0， 2π］上等间隔采样N=16时， N点IDFT［XN(k)］得到的序列正是原序列x(n)以16为周期进行周期延拓后的主值区序列

当N=16时，由于N，所以发生了时域混叠失真

当N=32时，由于N>M，频域采样定理，所以不存在时域混叠失真

Ifft函数用法同fft函数

三、 思考题简要分析

1、 如果序列x(n)的长度为M，希望得到其频谱在上的N点等间隔采样，当N时， 如何用一次最少点数的DFT得到该频谱采样？

答：先对原序列x(n)以N为周期进行周期延拓后取主值区序列，

再计算N点DFT则得到N点频域采样：

实验三：用FFT对信号作频谱分析

一、 实验原理与方法

1、周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。如果不知道信号周期，可以尽量选择信号的观察时间长一些。

2、对模拟信号进行谱分析时，首先要按照采样定理将其变成时域离散信号。如果是模拟周期信号，也应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的谱分析进行。

3、频谱分辨率直接和FFT的变换区间N有关，因为FFT能够实现的频率分辨率是，因此要求。

二、 实验结果及程序清单（含分析及函数用法）

1、clear all;close all
%---------------实验内容(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(3,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(3,2,2);mstem(X1k16); %绘制16点DFT的幅频特性图
title('(1b)16点DFT[x_1(n)]');xlabel('ω/π');ylabel('幅度');
axis([0,2,0,1.2*max(abs(X1k16))])
subplot(3,2,3);mstem(X2k8); %绘制8点DFT的幅频特性图
title('(2a) 8点DFT[x_2(n)]');xlabel('ω/π');ylabel('幅度');
axis([0,2,0,1.2*max(abs(X2k8))])
subplot(3,2,4);mstem(X2k16); %绘制16点DFT的幅频特性图
title('(2b)16点DFT[x_2(n)]');xlabel('ω/π');ylabel('幅度');
axis([0,2,0,1.2*max(abs(X2k16))])
subplot(3,2,5);mstem(X3k8); %绘制8点DFT的幅频特性图
title('(3a) 8点DFT[x_3(n)]');xlabel('ω/π');ylabel('幅度');
axis([0,2,0,1.2*max(abs(X3k8))])
subplot(3,2,6);mstem(X3k16); %绘制16点DFT的幅频特性图
title('(3b)16点DFT[x_3(n)]');xlabel('ω/π');ylabel('幅度');
axis([0,2,0,1.2*max(abs(X3k16))])

2、%-------------实验内容（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
subplot(2,2,1);mstem(X4k8); %绘制8点DFT的幅频特性图
title('(a) 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('(b)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('(a) 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('(b)16点DFT[x_5(n)]');xlabel('ω/π');ylabel('幅度');
axis([0,2,0,1.2*max(abs(X5k16))])

分析：

的周期为8，所以N=8和N=16均是其周期的整数倍，得到正确的单一频率正弦波的频谱，仅在0.25π处有1根单一谱线。如图（a）和（b）所示。

的周期为16，所以N=8不是其周期的整数倍，得到的频谱不正确，如图（a）所示。N=16是其一个周期，得到正确的频谱，仅在0.25π和0.125π处有2根单一谱线, 如图（b）所示。

3、%----------实验内容（3）---------
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),'.');box on %绘制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),'.');box on %绘制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),'.'); box on%绘制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))])

分析：

x6(t)有3个频率成分， f1=4 Hz, f2=8 Hz, f3=10 Hz。所以x6(t)的周期为0.5 s。 采样频率Fs=64 Hz=16f1=8f2=6.4f3。变换区间N=16时，观察时间Tp=16T=0.25 s，不是x6(t)的整数倍周期， 所以所得频谱不正确，如图8.3.1（6a）所示。变换区间N=32，64 时，观察时间Tp=0.5 s，1 s，是x6(t)的整数周期，所以所得频谱正确

补充：上述三个程序中用到的子程序：

function mstem(Xk)
M=length(Xk);
k=0:M-1; wk=2*k/M;
stem(wk,abs(Xk),'.'); box on;
xlabel('ω/π'); ylabel('幅度');
axis([0,2,0,1.2*max(abs(Xk))]);

三、 思考题简要分析

（1） 对于周期序列，如果周期不知道，如何用FFT进行谱分析？

答：可先截取M点进行DFT,再将截取长度扩大1倍，比较两次的结果。如果二者的主谱差别满足分析误差要求，则以两者中的一个近似表示周期序列的频谱，否则，继续把截取长度加倍，重复上述步骤。

（2） 如何选择FFT的变换区间？（包括非周期信号和周期信号）

答1、对于非周期信号：有频谱分辨率F，而频谱分辨率直接和FFT的变换区间有关，因为FFT能够实现的频率分辨率是2π/N...因此有最小的N>2π/F。就可以根据此式选择FFT的变换区间。

2、对于周期信号，周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。

（3）当N=8时，和的幅频特性会相同吗？为什么？N=16 呢？

答 ：不同，因为这样会影响是不是周期的整数倍的问题，即影响了频谱的正确性。

实验四：IIR数字滤波器设计及软件实现

一、 实验原理与方法

1、 IIR数字滤波器间接法基本设计过程：

将给定的数字滤波器的指标转换成过渡模拟滤波器的指标；

②设计过渡模拟滤波器；

③将过渡模拟滤波器系统函数转换成数字滤波器的系统函数

二、 实验结果及程序清单（含分析及函数用法）

1、----- --------三路调幅信号st------------------------

%-----------------信号产生函数mstg---------------
function st=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(2,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(2,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('幅度')

2、------- 低通滤波器设计与实现--------------------------

clear all;close all
Fs=10000;T=1/Fs; %采样频率
%调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st
st=mstg;

fp=280;fs=450; %下面wp,ws,为fp,fs的归一化值范围为0-1
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
[h,w]= freqz(B,A);
y1t=filter(B,A,st); %滤波器软件实现
figure(2);subplot(2,1,1);
plot(w,20*log10(abs(h)));
axis([0,1,-80,0])
subplot(2,1,2);
t=0:T:(length(y1t)-1)*T;
plot(t,y1t);
%axis([0,1,-80,0])

3、-------------带通滤波器设计与实现----------------------

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
[h,w]= freqz(B,A);
y2t=filter(B,A,st);
figure(3);subplot(2,1,1);
plot(w,20*log10(abs(h)));
axis([0,1,-80,0])
subplot(2,1,2);
t=0:T:(length(y2t)-1)*T;
plot(t,y2t);

4、-----------------高通滤波器设计与实现---------------

fp=900;fs=550;
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
[h,w]= freqz(B,A);
y3t=filter(B,A,st);
figure(4);subplot(2,1,1);
plot(w,20*log10(abs(h)));
axis([0,1,-80,0])
subplot(2,1,2);
t=0:T:(length(y3t)-1)*T;
plot(t,y3t);

分析：

由图可见，三个分离滤波器指标参数选取正确，损耗函数曲线达到所给指标。分离出的三路信号y1(n)、 y2(n)、 y3(n)的波形是抑制载波的单频调幅波。

Ellipord函数的调用:

[n,Wp]=ellipord(Wp,Ws,Rp,Rs)

用于计算满足指标的椭圆数字滤波器的最低阶数N和通带边界频率wp,其中wp,ws,rp,rs为四个基本指标，需要注意的是wp,ws应为归一化之后的值。

三、 思考题简要分析

（1） 请阅读信号产生函数mstg，确定三路调幅信号的载波频率和调制信号频率。

答：第1路调幅信号的载波频率fc1=1000Hz,调制信号频率fm1=100Hz

第2路调幅信号的载波频率fc2=500Hz，调制信号频率fm2=50Hz

第3路调幅信号的载波频率fc3=250Hz,调制信号频率fm3=25Hz

（2）信号产生函数mstg中采样点数N=800，对st进行N点FFT可以得到6根理想谱线。如果取N=1000，可否得到6根理想谱线？为什么？N=2000呢？请改变函数mstg中采样点数N的值，观察频谱图验证您的判断是否正确。

答：， s(t)的每个频率成分都是25 Hz的整数倍。 采样频率Fs=10 kHz=25×400 Hz，即在25 Hz的正弦波的1个周期中采样400点。 所以， 当N为400的整数倍时一定为s(t)的整数个周期。 因此，采样点数N=1800时， 对s(t)进行N点FFT不能得到6根理想谱线。 如果取N=2000，能得到6根理想谱线。

实验五：FIR数字滤波器设计与软件实现

一、 实验原理与方法

1、窗函数法设计线性相位低通滤波器设计步骤：

1、选择窗函数的类型，并估计窗口长度N

2、构造希望逼近的频率响应函数()

3、计算(n)

4、加窗得到设计结果：h(n)=(n)w(n)

二、 实验结果及程序清单（含分析及函数用法）

1、%--------------------信号产生函数------------------------
function xt=xtg
%信号x(t)产生,并显示信号的幅频特性曲线
%xt=xtg(N) 产生一个长度为N,有加性高频噪声的单频调幅信号xt,采样频率Fs=1000Hz
%载波频率fc=Fs/10=100Hz,调制正弦波频率f0=fc/10=10Hz.
N=1000;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(2,1,1);plot(t,xt);grid;xlabel('t/s');ylabel('x(t)');
axis([0,Tp/5,min(xt),max(xt)]);title('(a) 信号加噪声波形')
subplot(2,1,2);plot(f,abs(fst)/max(abs(fst)));grid;title('(b) 信号加噪声的频谱')
axis([0,Fs/2,0,1.2]);xlabel('f/Hz');ylabel('幅度')

2、clear all
xt=xtg;
N=1000;Fs=1000;T=1/Fs;Tp=N*T;
k=0:N-1;f=k/Tp;
t=0:T:(N-1)*T;
fp=120;fs=150;Rp=0.1;As=60;Fs=1000;
wc=(fp+fs)/Fs;
B=2*pi*(fs-fp)/Fs;
M=ceil(11*pi/B);
hn=fir1(M-1,wc,blackman(M));
Hw=abs(fft(hn,N));
ywt=fftfilt(hn,xt,N);
figure;
subplot(2,1,1);
plot(f,20*log10(Hw)/max(Hw));grid on
xlabel('f/Hz');ylabel('幅度（dB）');
title('(a)低通滤波器的幅频特性')
axis([0,500,-160,5]);
subplot(2,1,2);
plot(t,ywt);grid on
xlabel('t/s');ylabel('y_1(t)');
title('(b)滤除噪声后的信号波形')

3、clear all
xt=xtg;
N=1024;Fs=1000;T=1/Fs;Tp=N*T;
k=0:N-1;f=k/Tp;
K=1000;
t=0:T:(K-1)*T;
fp=120;fs=150;Rp=0.2;As=60;
fb=[fp,fs];
m=[1,0];
dev=[(10^(Rp/20)-1)/(10^(Rp/20)+1),10^(-As/20)];
[Ne,fo,mo,W]=remezord(fb,m,dev,Fs);
hn=remez(Ne,fo,mo,W);
Hw=abs(fft(hn,1024));
yet=fftfilt(hn,xt,N);
figure;
subplot(2,1,1);
plot(f,20*log10(Hw)/max(Hw));grid on
xlabel('f/Hz');ylabel('幅度（dB）');
title('(a)低通滤波器的幅频特性')
axis([0,500,-80,5]);
subplot(2,1,2);
plot(t,yet);grid on
xlabel('t/s');ylabel('y_2(t)');
title('(b)滤除噪声后的信号波形')

三、 思考题简要分析

(1)如果给定通带截止频率和阻带截止频率以及阻带最小衰减,如何用窗函数法设计线性相位低通滤波器?请写出设计步骤.

答：1、选择窗函数的类型，并估计窗口长度N

2、构造希望逼近的频率响应函数()

3、计算(n)

4、加窗得到设计结果：h(n)=(n)w(n)

(2)如果要求用窗函数法设计带通滤波器,且给定通带上、下截止频率为和，阻带上、下截止频率为和，试求理想带通滤波器的截止频率。

答：

（3）解释为什么对同样的技术指标，用等波纹最佳逼近法设计的滤波器阶数低？

答：①用窗函数法设计的滤波器，如果在阻带截止频率附近刚好满足，则离开阻带截止频率越远，阻带衰减富裕量越大，即存在资源浪费；

② 几种常用的典型窗函数的通带最大衰减和阻带最小衰减固定，且差别较大，又不能分别控制。所以设计的滤波器的通带最大衰减和阻带最小衰减通常都存在较大富裕。如本实验所选的blackman窗函数，其阻带最小衰减为74dB,而指标仅为60dB。

③ 用等波纹最佳逼近法设计的滤波器，其通带和阻带均为等波纹特性，且通带最大衰减和阻带最小衰减可以分别控制，所以其指标均匀分布，没有资源浪费，所以期阶数低得多。

实验六 数字信号处理在双音多频拨号系统中的应用

          实验内容

    ①  运行仿真程序exp6.m，任意送入6位电话号码，打印出相应的幅度谱。观察程序运行结果，对照表10.10.1和表10.10.2，判断程序谱分析的正确性。

    ②  分析该仿真程序，将产生、检测和识别6位电话号码的程序改为能产生、检测和识别8位电话号码的程序，并运行一次，打印出相应的幅度谱和8位电话号码。

==================

实验6程序：exp6.m

% 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=[1209,1336,1477,1633];               % 列频率向量

TN=input('键入6位电话号码= ');          % 输入6位数字

TNr=0;                                  %接收端电话号码初值为零

for l=1:6;

    d=fix(TN/10^(6-l));

    TN=TN-d*10^(6-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(3,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^(6-l);

end

disp('接收端检测到的号码为：')         % 显示接收到的字符

disp(TNr)

运行程序，根据提示键入6位电话号码123456，回车后可以听见6位电话号码对应的DTMF信号的声音，并输出相应的6幅频谱图如图10.10.1所示，左上角的第一个图在k=18和k=31两点出现峰值，所以对应第一位号码数字1。最后显示检测到的电话号码123456。

 2、实验内容②  只要对6位电话号码检测程序exp6.m作如下修改，即可产生、检测和识别8位电话号码。

（1）将第8行改为TN=input('键入8位电话号码= ');

（2）将第10~12行改为

for l=1:8;

  d=fix(TN/10^(1-8)); 

  TN=TN-d*10^(1-8); 

（3）将第26行改为 subplot(4,2,l);

（4）将第36行改为TNr=TNr+tm(r,s-4)*10^(1-8);

修改后的程序为exp6_8.m，程序清单见程序集。运行程序exp6_8.m，输入输入8位电话号码87654321，则输出相应的8幅频谱图如图10.10.2所示。最后显示检测到的电话号码12345678。

本文来源：https://www.2haoxitong.net/k/doc/ddfb02e5d0d233d4b14e69d8.html