数字信号处理实验报告一系统响应及系统稳定性
数字信号实验报告
实验名称 ： 数字信号处理第一次实验 指导老师 ： 学院 ： 信息工程学院 班级 ： 12电子信息工程2班 姓名 ：
学号 ：
实验一：系统响应及系统稳定性
实验目的
1. 掌握求系统响应的方法
2. 掌握时域离散系统的时域特性
3. 分析观察及检验系统的稳定性
实验原理
在时域中，描写系统特性的方法是差分方程和单位脉冲响应，在频域可以用系统函数描述系统特性。已知输入信号可以由差分方程/单位脉冲响应或系统函数求出系统对于该输入信号的响应。本实验仅在时域求解。在计算机上适合用递推法求差分方程的解，最简单的方法是采用MA TLAB 语言工具箱函数filter 函数。也可以用MA TLAB 语言的工具箱函数conv 函数计算输入信号和系统的单位脉冲响应的线性卷积，求出系统的响应。系统的时域特性指的是系统的线性时不变性质/因果性和稳定性。重点分析实验系统的稳定性，包括观察系统的暂态响应和稳定响应。系统的稳定性是指对任意有界的输入信号，系统都能得到有界的系统响应。或者系统的单位脉冲响应满足绝对可和的条件。系统的稳定性由其差分方程的系数决定。实际中检查系统是否稳定，不可能检查系统对所有有界的输入信号，输出是否都是有界输出，或者检查系统的单位脉冲响应满足绝对可和的条件。可行的方法是在系统的输入端加入单位阶跃序列，如果系统的输出趋近于一个常数（包括零），就可以断定系统是稳定的。系统的稳态输出是指当n 趋向于无穷时，系统的输出。如果系统稳定，则信号加入系统后，系统输出的开始一段称为暂态效应，随着n 的加大，幅度趋于稳定，达到稳态输出。
实验内容及步骤
1. 编制程序，包括产生输入信号、单位脉冲响应序列的子程序，用filter 函数或conv 函数求解系统输出响应的主程序。程序中要有绘制信号波形的功能。 2. 给定一个地通滤波器的差分方程为y (n )=0. 05x (n )+0. 05x (n -1) +0. 9y (n -1) 输入信号：
x 1(n ) =R 8(n ) x 2(n ) =μ(n )
①分别求出x 1(n ) =R 8(n ) 和x 2(n ) =μ(n ) 的系统响应y 1(n ) 和y 2(n ) ，并画出其
波形
②求出系统的单位脉冲响应，画出其波形。
（1） 给定系统的单位脉冲响应为 h 1(n ) =R 10(n )
h 2(n ) =δ(n ) +0. 25δ(n -1) +2. 5δ(n -2) +δ(n -3)
用线性卷积法求x 1(n ) =R 8(n ) 分别对系统h 1(n ) 和h 2(n ) 的输出响应y 21(n ) 和y 22(n ) 并画出波形。
（2） 给定一谐振器的差分方程为
谐振器y (n ) =1. 8237y (n -1) -0. 9801y (n -2) +b 0x (n ) -b 0x (n -2) 令b 0=1/100. 49，
的谐振频率为0.4rad
①用实验方法检查系统是否稳定，输入信号为μ(n ) 时，画出系统输出波形y 31(n ) 。 ②给定输入信号为x (n ) =sin(0. 014n ) +sin(0. 4n ) ，求出系统的输出响应y 32(n ) ，并画出其波形。
实验程序
实验1—1程序：
close all ; clear all ; A=[1,-0.9]; B=[0.05,0.05];
x1n=[1 1 1 1 1 1 1 1 zeros(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)'); box on
y1n=filter(B,A,x1n); subplot(2,2,2);y='y1(n)'; tstem(y1n,y);
title('(b) 系统对R8(n)的响应y1(n)'); box on
y2n=filter(B,A,x2n); subplot(2,2,4);y='y2(n)'; tstem(y2n,y);
title('(c) 系统对u(n)的响应y2(n)'); box on
波形图：
(a) 系统单位脉冲响应
h(n)(b) 系统对R8(n)的响应y1(n)
y n
n
y n
n
(c) 系统对u(n)的响应y2(n)y n
n
实验1—2程序：
close all ; clear all ;
x1n=[1 1 1 1 1 1 1 1 ]; 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);
title('(d) 系统单位脉冲响应h1(n)'); box on
subplot(2,2,2);y='y21(n)'; tstem(y21n,y);
title('(e) h1(n)与R8(n)的卷积y21(n)'); box on
subplot(2,2,3);y='h2(n)'; tstem(h2n,y);
title('(f) 系统单位脉冲响应h2(n)'); box on
subplot(2,2,4);
y='y22(n)'; tstem(y22n,y);
title('(g) h2(n)与R8(n)的卷积y22(n)'); box on
波形图：
(d) 系统单位脉冲响应h1(n)
(e) h1(n)与R8(n)的卷积y21(n)
y n
y n
n
(f) 系统单位脉冲响应
h2(n)
n
(g) h2(n)与R8(n)的卷积y22(n)y n
y n
n
n
实验1—3程序：
close all; clear all; un=ones(1,256); 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]; y31n=filter(B,A,un); y32n=filter(B,A,xsin); figure(3) subplot(2,1,1); y='y31(n)'; tstem(y31n,y);
title('(h) 谐振器对u(n)的响应y31(n)'); box on
subplot(2,1,2);
y='y32(n)'; tstem(y32n,y);
title('(i) 谐振器对正弦信号的响应y32(n)'); box on
波形图：
(h) 谐振器对u(n)的响应
y31(n)
y n
n
(i) 谐振器对正弦信号的响应y32(n)
y n
n
分析讨论：
（1） 综合起来，在时域求系统响应的方法有两种，第一是通过解差分方程求得系统输出，
注意合理地选择初始条件；第二种是已知系统的单位脉冲响应，通过求输入信号和系统单位脉冲响应的线性卷积求得系统输出。用计算机求解时最好使用MATLAB 语言进行。
（2） 实际中要检验系统的稳定性，其方法是在输入端加入单位阶跃序列，观察输出波形，
如果波形稳定在一个常数值上，系统稳定，否则不稳定。如第三个实验就是稳定的。
（3） 谢正奇具有对某个频率进行谐振的性质，本实验中的谐振器的谐振频率是0.4rad ，
因此稳定波形为sin(0.4n)。
思考题
思考题1-1
如果输入信号为无限长序列，系统的单位脉冲响应是有限长序列，可用分段线性卷积法求系统的响应。
思考题1-2
如果信号经过经过低通滤波器，则信号的高频分量被滤掉，时域信号的变化减缓，在有阶跃处附近产生过渡变化时间。因此，当输入矩形序列时，输出序列的开始和终了都产生了明显的上升和下降时间，详细可见第一个实验结果的波形。
实验报告要求
（1） （2） （3） （4） （5）
简述在时域求系统响应的方法
简述通过实验判断系统稳定性的方法。分析上面第三个实验的稳定输出的波形。 对各实验所得的结果进行简单的分析和解释。 简要回答思考题
打印程序清单和要求的各信号波形。
实验二：时域采样与频域采样
实验目的
时域采样理论与频域采样理论是数字信号处理中的重要理论。要求掌握模拟信号采样前后的频谱变化，以及如何选择采样频率才能使采样后的信号不丢失信息；要求掌握频域采样会引起时域周期化的概念，以及频率域采样定理及其对频域采样点数选择的指导作用。
实验原理
1. 时域采样定理的要点
(1)对模拟信号x a (t ) 以T 进行时域等间隔理想采样，形成的采样信号的频谱X (j Ω) 会以采样角频率Ωs (Ωs =2π/T ) 为周期进行周期延拓。公式为
1∞
X a (j Ω) =FT [x a (t ) ]=∑X a (j Ω-jn Ωs )
T n =-∞
∧
∧
(2)采样频率Ωs 必须大于等于模拟信号最高频率的两倍以上，才能使采样信号的频谱不产生频谱混叠。
利用计算机计算X a (j Ω) 并不方便，下面我们导出另外一个公式，以便在计算机上进行实验。
理想采样信号x a (t ) 和模拟信号x a (t ) 之间的关系为
∧
∧
x a (t ) =x a (t ) ∑δ(t -nT ) 对上式进行傅立叶变换，得到
n =-∞
∧
∞
X a (j Ω) =⎰[x a (t ) ∑δ(t -nT ) ]e
-∞
n =-∞
∧
∞
∞
-j Ωt
dt =
n =-∞
∑⎰
∞
∞
-∞
x a (t ) δ(t -nT ) e -j Ωt dt
在上式的积分号内只有当t=nT时，才有非零值，因此
X a (j Ω) =
∧
n =-∞
∑x (nT ) e
a
∞
-j ΩnT
上式中，在数值上x a (nT ) =x (n ) ，再将ω=Ω代入，得到
X a (j Ω) =
∧
n =-∞
∑x (n ) e
∞
-jwn
上式的右边就是序列的傅立叶变换X (e ) ，即
jw
X a (j Ω) =X (e jw ) |w =ΩT
上式说明理想 的傅立叶变换可用相应的采样序列的傅立叶变换得到，只要将自变量ω用ΩT 代替即可。 2. 频域采样定理的要点
（1）对信号x(n)的频谱函数X (e ) 在[0,2π]上等间隔采样N 点，得到
jw
∧
X N (k ) =X (e jw ) |
w =
2k πN
（ k=0,1,2…N-1），则N 点IDFT [X N (k )]得到的序列就是原序列x(n)以N 为周期进行周期延拓后的主值区序列，公式x N (n ) =IDFT [X N (k )]N =[
i =-∞
∑x (n +iN ) ]R
∞
N
(n )
（2）由上式可知，频率采样点数N 必须大于等于时域离散信号的长度M （即N ≥M ），才能使时域不产生混叠，这时N 点
[IDFT (X N k )]得到的序列x N (n ) 就是原序列x (n ) ，即
x N (n ) =x (n ) 。如果N>M，则x N (n ) 比原序列尾部多了N-M 个零点；如果N
而且N
因此
x N (n ) 与x (n ) 不相同。
实验程序
实验2—1程序： close all; clear all; Tp=64/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); 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(Xk));
title('(a) T*FT[xa(nT)],Fs=1000Hz'); xlabel('f(Hz)');ylabel('幅度'); axis([0,Fs,0,1.2*max(abs(Xk))]) Tp=64/300; 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); 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(Xk));
title('(a) T*FT[xa(nT)],Fs=300Hz'); xlabel('f(Hz)');ylabel('幅度'); axis([0,Fs,0,1.2*max(abs(Xk))]) Tp=64/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); 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(Xk));
title('(a) T*FT[xa(nT)],Fs=200Hz'); xlabel('f(Hz)');
ylabel('幅度');
axis([0,Fs,0,1.2*max(abs(Xk))]) 波形图：
(a) Fs=1000Hz
n
(a) Fs=300Hz
n
(a) Fs=200Hz
11
(a) T*FT[xa(nT)],Fs=1000Hz
幅度
y n
0.500
5001000f(Hz)
(a) T*FT[xa(nT)],Fs=300Hz
幅度
y n
0.500
100
200
300
f(Hz)
(a) T*FT[xa(nT)],Fs=200Hz
幅度
n
y n
0.500
50
100f(Hz)
150
200
实验2—2程序：
M=27; N=32; n=0:M;
xa=0:floor(M/2); xb= ceil(M/2)-1:-1:0; xn=[xa,xb]; Xk=fft(xn,1024); X32k=fft(xn,32); x32n=ifft(X32k); X16k=X32k(1:2:N); x16n=ifft(X16k,N/2); 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])
波形图：
(a)FT[x(n)]
200
(b) 三角波序列x(n)
|X (e j ω) |
10000
0.5ω/π
(c) 16点频
域采样
1
x (n )
n
(d) 16点
IDFT[X16(k)]
|X 16(k ) |
k
x 16(n )
10
20
30
(e) 32点频域采样
n
(f) 32点IDFT[X32(k)]
|X 32(k ) |
k
x 32(n )
n
分析讨论：
（1） 时域采样定理的验证程序的运行结果如图实验2-1的波形图所示。由图可见，当采
样频率为1000Hz 时，频谱混叠很小；当采样频率为300Hz 时，频谱混叠很严重；当采样频率为200Hz 时，频谱混叠更严重。
（2） 频域采样定理的验证程序的运行结果如图实验2-2的波形图所示。该图验证了频域
采样定理和频域采样理论。对信号x(n)的频谱函数X (e ) 在［0,2π］上等间隔采
jw
样N=16时，N 点IDFT ［X N (k)］得到的序列正是原序列x(n)以16为周期进行周期延拓后的主值区序列，由于N
思考题
先对原序列x (n ) 以N 为周期进行周期延拓后去主值区序列，
x N (n ) =[∑x (n +Ni ) ]R N (n )
i =-∞
∞
再计算N 点的DFT ，则得到N 点的频域采样：
X N (k ) =DFT [x N (n )]N =X (e jw ) |
w =
2πk N
k =0, 1, 2,..., N -1
实验报告要求
（1） （2） （3） （4）
运行程序，打印要求显示的图形。
对各实验所得的结果进行简单的分析和解释。 简要回答思考题。
打印程序清单和要求的各信号波形。
实验三：用FFT 对信号作频谱分析
实验目的
学习用FFT 对连续信号的时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便正确应用FFT 。
实验原理
用FFT 对信号作频谱分析是学习数字信号处理的重要内容。经常需要进行谱分析的信号是模拟信号和时域离散信号。对信号进行谱分析的重要问题是频谱分辨率D 和分析误差。频谱分辨率直接和FFT 的变换区间N 有关，因为FFT 能够实现的频率分辨率是2∏/N，因此要求2∏/N≤D 。可以根据此式选择FFT 的变换区间N 。误差主要来自于用FFT 频谱分析时，得到的是离散谱，而信号（周期信号除外）是连续谱的，只有当N 较大时，离散谱的包络才能逼近于连续谱，因此N 要适当选择大一些。
周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT ，得到离散谱才能代表周期信号的频谱。如果不知道信号周期，可以尽量选择信号的观察时间长一些。
对模拟信号进行谱分析时，首先要按照采样定理将其变成时域离散信号。如果是模拟周期信号，也应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的谱分
析进行。
实验步骤及内容
（1） 对以下序列进行谱分析：
x 1(n ) =R 4(n )
⎧n +1(0≤n ≤3) ⎪
x 2(n ) =⎨8-n (4≤n ≤7)
⎪0(其他n ) ⎩⎧4-n (0≤n ≤3) ⎪
x 3(n ) =⎨n -3(4≤n ≤7)
⎪0(其他n ) ⎩
选择FFT 的变换区间N 为8和16两种情况进行频谱分析。分别打印其幅频特性曲线，并进行对比，分析和讨论。
（2） 对以下周期序列进行谱分析：
x 4(n ) =cos x 5(n ) =cos
π
4
n
πn
4
+cos
πn
8
选择FFT 的变换区间N 为8和16两种情况分别对以上序列进行频谱分析。分别打印其幅频特性曲线，并进行对比，分析和讨论。 （3） 对模拟周期信号进行谱分析：
x 6(t ) =cos 8πt +cos 16πt +cos 20πt
选择采样频率F s =64Hz ，对变换区间N=16，32，64三种情况进行谱分析。分别打印其幅频特性曲线，并进行分析和讨论。
实验程序
实验3—1程序： clear all; close all;
x1n=[ones(1,4)]; M=8;xa=1:(M/2); xb=(M/2):-1:1; x2n=[xa,xb]; x3n=[xb,xa];
X1k8=fft(x1n,8); X1k16=fft(x1n,16);
X2k8=fft(x2n,8); X2k16=fft(x2n,16); X3k8=fft(x3n,8); X3k16=fft(x3n,16); subplot(2,2,1); mstem(X1k8);
title('(1a)8点DFT[x_1(n)]'); xlabel('ω/π');
ylabel('幅度');
axis([0,2,0,1.2*max(abs(X1k8))]) subplot(2,2,2); mstem(X1k16);
title('(1b)16点DFT[x_1(n)]'); xlabel('ω/π');
ylabel('幅度');
axis([0,2,0,1.2*max(abs(X1k16))]) subplot(2,2,3); mstem(X2k8);
title('(2a)8点DFT[x_2(n)]'); xlabel('ω/π');
ylabel('幅度');
axis([0,2,0,1.2*max(abs(X2k8))]) subplot(2,2,4); mstem(X2k16);
title('(2b)16点DFT[x_2(n)]'); xlabel('ω/π');
ylabel('幅度'); figure(2)
axis([0,2,0,1.2*max(abs(X2k16))]) subplot(2,2,1); mstem(X3k8);
title('(3a) 8点DFT[x_3(n)]'); xlabel('ω/π');
ylabel('幅度');
axis([0,2,0,1.2*max(abs(X3k8))]) subplot(2,2,2); mstem(X3k16);
title('(3b)16点DFT[x_3(n)]'); xlabel('ω/π'); ylabel('幅度');
axis([0,2,0,1.2*max(abs(X3k16))])