数字信号处理课程设计
题目： 试实现线性常系数差分方程的求解 学院： 专业：
班级：
学号：
组员：
指导教师：
题目：用Matlab 实现线性常系数差分方程求解
一． 设计要求
1． 2． 3．
掌握线性常系数差分方程的求解
熟练掌握Matlab 基本操作和各类函数调用 结合Matlab 实现线性常系数差分方程的求解
二．设计原理
1．差分与差分方程
与连续时间信号的微分及积分运算相对应，离散时间信号有差分及序列求和运算。设有序列f(k),则称…,f(k+2),f(k+1),…,f(k－1),f(k－2), …为f(k)的移位序列。序列的差分可以分为前向差分和后向差分。一阶前向差分定义为
∆f (k ) =f (k +1) -f (k ) （3.1—1）
一阶后向差分定义为
∆f (k ) =f (k ) -f (k -1) （3.1—2）
式中Δ和Δ称为差分算子。由式（3.1—1）和式（3.1—2）可见，前向差分与后向差分的关系为
∆f (k ) =∆f (k -1) （3.1—3）
二者仅移位不同，没有原则上的差别，因而它们的性质也相同。此处主要采用后向差分，并简称其为差分。
由查分的定义，若有序列f 1(k ) 、f 2(k ) 和常数a 1，a 2则
∆[a 1f 1(k ) +a 2f 2(k )]=[a 1f 1(k ) +a 2f 2(k )]-[a 1f 1(k -1) +a 2f 2(k -1)]=a 1[f 1(k ) -f 1(k -1)]+a 2[f 2(k ) -f 2(k -1)]=a 1∆f 1(k ) +a 2∆f 2(k )
这表明差分运算具有线性性质。 二阶差分可定义为
（3.1—4）
∆2f (k ) =∆[∆f (k )]=∆[f (k ) -f (k -1)]=∆f (k ) -∆f (k -1) =f (k ) -2f (k -1) +f (k -2)
（3.1—5）
类似的，可定义三阶、四阶、…、n 阶差分。一般地，n 阶差分
∆f (k ) =∆[∆
式中
n n -1
⎛n ⎫
f (k )]=∑(-1) j ⎪f (k -j ) （3.1—6）
j =0⎝j ⎭
n
⎛n ⎫n !
=, j =0,1,2, ⎪
⎝j ⎭(n -j )! j !
为二项式系数 序列f(k)的求和运算为
, n （3.1—7）
i =-∝
∑f (i ) （3.1—8）
k
差分方程是包含关于变量k 的未知序列y(k)及其各阶差分的方程式，它的一般形式可写为
F ⎡⎣k , y (k ), ∆y (k ), , ∆n y (k ) ⎤⎦=0 （3.1—9a ）
式中差分的最高阶为n 阶，称为n 阶差分方程。由式（3.1—6）可知，各阶差分均可写为y(k)及其各移位序列的线性组合，故上式常写为
G [k , y (k ), y (k -1), , y (k -n ) ]=0 （3.1—9b ）
通常所说的差分方程是指式（3.1—9b ）形式的方程。
若式（3.1—9b ）中，y(k)及其各移位序列均为常数，就称其为常系数差分方程；如果某些系数是变量k 的函数，就称其为变系数差分方程。描述LTI 离散系统的是常系数线性差分方程。
差分方程是具有递推关系的代数方程，若一直初始条件和激励，利用迭代法渴求的差分方程的数值解。
2. 差分方程的经典解
一般而言，如果但输入—单输出的LTI 系统的激励f(k)，其全响应为y(k)，那么，描述该系统激励f(k)与响应y(k)之间关系的数学模型式n 阶常系数线性差分方程，它可写为
y (k ) +a n -1y (k -1) ++a 0y (k -n )
+b 0f (k -m )
=b m f (k ) +b m -1f (k -1) +
式中a j (j =0,1,
（3.1—10a ）
, n -1) 、b i (i =0,1,
m
, m ) 都是常数。上式可缩写为
∑a
j =0
n
n -j
y (k -j ) =∑a m -i f (k -i )(式中a n =1) （3.1—10b ）
i =0
与微分方程的经典解类似，上述差分方程的解由齐次解和特解两部分组成。齐次解用y h (k ) 表示，特解用y p (k ) 表示，即
y (k )=y h (k ) +y p (k ) （3.1—11）
a. 齐次解
当式（3.1—10) 中的f(k)及其各移位项均为零时，齐次方程
y (k ) +a n -1y (k -1) +的解称为齐次解。
+a 0y (k -n ) =0 （3.1—12）
首先分析最简单的一阶差分方程。若一阶差分方程的齐次方程为
y (k ) +ay (k -1) =0 （3.1—13）
它可改写为
y (k )
=-a
y (k -1)
y(k)与y(k－1) 之比等于－a 表明，序列y(k)是一个公比为－a 的等比级数，因此y(k)应有如下形式
y (k ) =C (-a ) k （3.1—14）
式中C 式常数，有初始条件确定。
k k
C λC λ对于n 阶齐次差分方程，它的齐次解由形式为的序列组合而成，将
代入到式（3.1—12），得
C λk +a n -1C λk -1+
+a 1C λk -n -1+a 0C λk -n =0
k -n
λ由于C ≠0，消去C ；且λ≠0，以除上式，得
λn +a n -1λn -1+
+a 1λ+a 0=0（3.1—15）
上式称为差分方程式（3.1—10）和式（3.1—12）的特征方程，它有n 个根
λj (j =0,1, , n )
，称为差分方程的特征根。显然，形式为
C j λk j
的序列都满足式
（3.1—12），因而它们是式（3.1—10）方程的齐次解。依特征根取值的不同，
C D A θ
差分方程齐次解的形式见表3—1, 其中j 、j 、j 、j 等为待定常数
表3—1 不同特征根所对应的齐次解
-θ0)]
b. 特解
特解的函数形式与激励的函数形式有关，表3—2列出了集中典型的激励f(k)所对应的特解
y p (k )
。选定特解后代入原差分方程，求出其待定系数
P j (A 或θ)
等，
就得出方程的特解。
表3—2 不同激励所对应的特解
c. 全解
式（3.1—10）的线性差分方程的全解是齐次解与特解之和。如果方程的特征根均为单根，则差分方程的全解为
y (k )=y h (k ) +y p (k ) =∑C j λj k +y p (k )
j =1
n
（3.1—16）
如果特征根λ1为r 重根，而其余n －r 个特征根为单根时，差分方程的全解为
y (k )=∑C j k r -j λj k +
j =1
r
j =r +1
∑C λ
j
n
k j
+y p (k )
（3.1—17）
式中各系数
C j
由初始条件确定。
如果激励信号是在k=0时接入的，差分方程的解适合于k ≥0。对于n 阶差分方程，用给定的n 个初始条件y(0),y(1),…,y(n－1) 就可确定全部待定系数件y(0),y(1),…,y(n－1) 分别代入到式（3.1—16），可得
。如
果差分方程的特解都是单根，则方程的全解为式（3.1—16），将给定的初始条
y (0)=C 1+C 2++C n +y p (0)+λn C n +y p (1)
+λn n -1C n +
y p (n -1)
（3.1—18）
y (1)=λ1C 1+λ2C 2+
y (n -1) =λ1n -1C 1+λ2n -1C 2+
由以上方程可求得全部待定系数C j (j
=0,1, , n ) 。
2.1零输入响应
系统的激励为零，仅由系统的初始状态引起的响应，称为零输入响应，用
y zi (k ) 表示。在零输入条件下，式（3.1—10）等号右端为零，化为齐次方程，
即
∑a
j =0
n
n -j
y zi (k -j ) =0 （3.1—25）
一般设定激励是在k=0时接入系统的，在k ＜0时，激励尚未接入，故式（3.1—25）的几个初始状态满足
y zi (-1) =y (-1) y zi (-2) =y (-2)
（3.1—26）
y zi (-n ) =
y (-n )
式（3.1—26）中的y(－1),y(－2), …,y(－n) 为系数的初始状态，由式（3.1—25）和式（3.1—26）可求得零输入响应
y zi (k ) 。
2.2零状态响应
当系统的初始状态为零，仅由激励f(k)所产生的响应，称为零状态响应，用 表示。在零状态情况下，式（3.1—10）仍是非齐次方程，其初始状态为零，即零状态响应满足
∑a
j =0
n
n -j
y zs (k -j ) =∑a m -i f (k -i )
i =0
m
y zs (-1) =y zs (-2) =
=y zs (-n ) =0（3.1—30）
的解。若其特征根均为单根，则其零状态响应为
y zs (k ) =∑C zsj λj k +y p (k )
j =1
n
（3.1—31）
C y (k ) 式中zsj 为待定常数，p 为特解。需要指出，零状态响应的初始状态
y zs (-1), y zs (-2), , y zs (-n ) 为零为零，但其初始值y zs (0),y zs (1),, y zs (n -1) 不一定等于零。
3. 线性常系数差分方程
3.1一个N 阶线性常系数差分方程可用下式表示：
y (n ) =∑b i x (n -i ) -∑a i y (n -i )
i =0
i =1
M N
（1.4.1）
或者
∑a y (n -i ) =∑b x (n -i ), a
i
i
i =0
i =0
N M
=1
（1.4.2）
式中，x （n ）和y(n)分别是系统的输入序列和输出序列，ai 和bi 均为常系数，式中y(n-i)和x(n-i)项只有一次幂，也没有相互交叉相乘项，故称为线性常系数差分方程。差分方程的阶数是用方程y(n-i)项中i 的最大取值与最小取值之差确定的。在（1.4.2）式中，y(n-i)项i 最大的取值N ，i 的最小取值为零，因此称为N 阶差分方程。
4． 线性常系数差分方程的求解
已知系统的输入序列，通过求解差分方程可以求出输出序列。求解差分方程的基本方法有以下三种：
（1） 经典解法。这种方法类似于模拟系统中求解微分方程的方法，它包括齐次解与
特解，由边界条件求待定系数，上节已作简单介绍，这里不作介绍。
（2） 递推解法。这种方法简单，且适合用计算机求解，但只能得到数值解，对于阶
次较高的线性常系数差分方程不容易得到封闭式（公式）解答。
（3） 变换域方法。这种方法是将差分方程变换到z 域进行求解，方法简便有效。
当然还可以不直接求解差分方程，而是先由差分方程求出系统的单位脉冲响应，再与已知的输入序列进行卷积运算，得到系统输出。但是系统的单位脉冲响应如果不是预
先知道，仍然需要求解差分方程，求其零状态响应解。
（4） 卷积法：由差分方程求出系统的h(n)，再与已知的x(n) 进行卷积，得到
y(n)。
观察（1.4.1）式，求n 时刻的输出，要知道n 时刻以及n 时刻以前的输入序列值，还要知道n 时刻以前的N 个输出序列值。因此求解差分方程在给定输入序列的条件下，还需要确定N 个初始条件。如果求n0时刻以后的输出，n0时刻以前N 个输出值y （n0-1）、y （n0-2）、∙∙∙∙∙∙、y （n0-N ）就构成了初始条件。
（1.4.1）式表明，已知输入序列和N 个初始条件，则可以求出n 时刻的输出；如果将该公式中的n 用n+1代替，可以求出n+1时刻的输出，因此（1.4.1） 式表示的差分方程本身就是一个适合递推法求解的方程。
三．设计过程
1. 用MATLAB 求解差分方程
MATLAB 信号处理工具箱提供的filter 函数实现线性常系数差分方程的递推
求解，调用格式如下：
yn=filter(B,A.xn) 计算系统对输入信号向量xn 的零状态响应输出信号向量yn ，yn 与xn 长度相等，其中，B 和A 是（1.4.2）式所给差分方程的系数向量，即
B=[b0,b1,∙∙∙,bM], A=[a0,a1,∙∙∙,aN]
其中a0=1，如果a0≠1，则filter 用a0对系数向量B 和A 归一化。
yn=filter(B,A.xn,xi) 计算系统对输入信号向量xn 的全响应输出信号yn 。所谓全响应，就是由初始状态引起的零输入响应和由输入信号xn 引起的零状态响应之和。其中，xi 是等效初始条件的输入序列，所以xi 是由初始条件确定的。
MATLAB 信号处理工具箱提供的filtic 就是由初始条件计算xi 的函数，其调用格式如下：
xi=filtic（B,A,ys,xs ）
其中，ys 和xs 是初始条件向量：ys=[y(-1),y(-2),y(-3),∙∙
∙,y(-N)],xs=[x(-1),x(-2),x(-3),∙∙∙,x(-M)]。如果xn 是因果序列，则xs=0.调用时可缺省xs 。
例1.4.1的MATLAB 求解程序ep141.m 如下：
%1.4.1.m:调用MATLAB 解差分方程y(n)-0.8y(n-1)=x(n)
a=0.8；ys=1; %设差分方程系数a=0.8，初始状态：y(-1)=1
xn=[1,zeros(1,30)]; %x(n)=单位脉冲序列，长度N=31 B=1;A=[1,-0.8]; %差分方程系数
xi=filtic(B,A,ys); %由初始条件计算等效初始条件的输入序列xi
yn=filter(B,A,xn,xi); %调用fiter 解差分方程，求系统输出信号y(n)
n=0:length(yn)-1;
stem(n,yn,'.' )
title(' 时域波形图)' );xlabel('n' );ylabel('y(n)')
程序中取查分方程系数a=0.8时，得到系统输出y （n ）如图1.4.1（a ）所示，与例1.4.1的解析递推结果完全相同。如果令初始条件y （-1）=0（仅修改程序中ys=0），则得到系统输出y （n ）=h（n ），如图1.4.1（b ）所示。
(a) (b)
图（a ）为a=0.8，y （-1）=1时，系统输出时域波形图，图（b ）为a=0.8，y(-1)=0时，系统输出时域波形图。
四．设计代码及结果
MATLAB 源程序
源程序如下
%1.m:调用MATLAB 解差分方程y(n)-0.8y(n-1)=x(n)
ys=1; %初始状态：y(-1)=1
xn=[1,zeros(1,30)]; %x(n)=单位脉冲序列，长度N=31
B=1;A=[1,-0.8]; %差分方程系数
xi=filtic(B,A,ys); %由初始条件计算等效初始条件的输入序列xi yn=filter(B,A,xn,xi); %调用filter 解差分方程，求系统输出信号y(n) n=0:length(yn)-1; %n的取值范围
stem(n,yn,'.' ) %画出时域波形图
title(' 时域波形图)' );xlabel('n' );ylabel('y(n)') %x轴、y 轴分别代表n ，x(n)
n=-5:5; %n的取值范围
xn=0.5.^n; %xn=0.5.^n
stem(n,xn,'fill' ),grid on %画出时域波形图
xlabel('n' ),ylabel('x(n)'), title('时域波形图') %x轴、y 轴分别代表n ，x(n)
n=-5:5; %n的取值范围
a=0.5; %设a=0.5
xen=a*[a.^n+a.^(-n)]; %xen=a*[a.^n+a.^(-n)]
xon=a*[a.^n-a.^(-n)]; %xon=a*[a.^n-a.^(-n)]
figure(1); stem(n,xen,'filled' ),grid on %画出xe(n)时域波形图
title(' 时域波形图' );xlabel('n' );ylabel('xe(n)') %x轴、y 轴分别代表n ，xe(n) figure(2); stem(n,xon,'filled' ),grid on %画出xo(n)时域波形图
title(' 时域波形图' );xlabel('n' );ylabel('xo(n)') %x轴、y 轴分别代表n ，xo(n)
b0=2;b2=0;b2=-1;a1=-0.7;a2=0.1;ys=0; %设差分方程系数，初始状态：y(-1)=1
B=[2,0,-1];A=[1,-0.7,0.1]; %差分方程系数
n=-5:5; %n的取值范围
xn=0.5.^n; %x(n)=0.5.^n
xi=filtic(B,A,ys); %由初始条件计算等效初始条件的输入序列xi yn=filter(B,A,xn,xi); %调用fiter 解差分方程，求系统输出信号y(n)
stem(n,yn,'.' ) %画出时域波形图
title('(a)');xlabel('n' );ylabel('y(n)') %x轴、y 轴分别代表n ，y(n)
五．程序运行结果如下：
差分序列时域波形图
输入参数a1=-0.7,a2=0.1,b0=2,b1=0,b2=-1 得到二阶线性常系数差分方程为
y(n)-0.7y(n-1)+0.1y(n-2)=2x(n)-x(n-2)
结果如下：
x （n ）的时域波形图
共轭对称分量xe （n ）的时域波形图
共轭对称分量xo （n ）的时域波形图
输入初始条件ys=y(-1)=1,得到结果如下：
输入初始条件ys=y(-1)=1时，输出y(n)波形图
改变初始条件ys=y(-1)=0,得到结果如下：
输入初始条件ys=y(-1)=0时，输出y(n)波形图
改变初始条件ys=y(-1)=100,得到结果如下:
输入初始条件ys=y(-1)=100时，输出y(n)波形图
六．比较结果总结
由上实验可知，通过改变初始条件ys 的值，得到的输出波形大小并不一致，即输出信号y(n)是不相同的；从而我们可以得出，对于同一个差分方程和同一个输入信号，因为初始条件不同，得到的输出信号是不相同的。
七．收获与体会
本次MATLAB 课程设计让我熟悉了该软件的一些功能，但是对于灵活应用MATLAB ，以及掌握各方面的设计思维以及技巧，还需要投入更多的时间。在熟悉MATLAB 程序和操作的同时培养了我的独立思考能力，专研精神，解决问题能力和动手能力。
在此之前了解到MATLAB 是一个很重要很有用的工具，但我并没有完全理解，本课程设计中，通过查阅资料，阅读网上程序并读写程序，对于MATLAB 的应用有了更深的了解，同时也认识到MATLAB 功能非常的强大，有着很多方
面的应用，如绘制函数，处理音频，图像数据，创建用户界面等功能，实为一个功能强大的软件。
本次课程设计我完成了基于MATLAB 的线性常系数差分方程求解的题目，通过实际操作回顾所学的内容，强化基础，实践理论知识。相信在以后的学习中，还会更加深入的了解MATLAB ，应用它。
随着课程设计报告的基本完成，本次课程设计终于接近了尾声。本次课程设计要求我们利用上学期所学的信号与线性系统分析的知识结合MATLAB 编程工具，完成差分方程求解设计的题目，通过实际操作，回顾所学内容，务实基础，强化理论知识，并体验理论与实际相结合的过程。
设计过程中遇到的第一个问题便是对于MATLAB 语言的不熟悉，其实现在想想这个问题不应该成为问题。毕竟本专业曾开设过MATLAB 程序设计这门课，而且老师还特别提醒过课程设计会用到MATLAB 语言。可惜的是，老师的话并没有引起我的足够重视，才导致要真正设计的时候一头雾水，无从下手。通过这件事，我明白了对于一件事情，想要做得很好，提前做功课和准备是必不可少的，机会只垂青有准备的人。
当然，通过本次课程设计，我还是基本熟悉了一些MATLAB 模块以及与本课程有关的一些函数的用法。但是由于上学期信号与线性系统分析课程学的不是很好，所以也就是在懂得同学的指导下按部就班的写了一些代码，也没有什么创新，但还是很好的回顾了差分方程求解的内容。老师在任务书中提到的目的要求，我自认为多多少少完成了一些。但不可否认的是还有很多没有达到的。总之，很多东西还是需要自己在不断摸索中找到答案。
同时，通过本次数字信号处理课程设计，使我全面系统的了解了差分方程求解的原理。通过结合课本的知识去完成课程设计也让我明白了理论与实践的相结合的重要性。只有理论知识是远远不够的，只有把所学的理论知识与实际相结合起来，从理论中得出结论，才是真正的知识，才能提高自己的实际动手能力和独立思考能力。从而加深了对课本知识的理解，同时也让我认识到了知识点方面自己的薄弱点。弄清楚了差分方程求解的设计及软件实现，总结了课程设计中的错误，提高了自己在实验过程中的效率和准确性。也把以前所学过的知识重新温故，巩固了所学的知识。通过这次实验，我对用MATLAB 求解差分方程的线性常系
数有了更深刻的理解，熟悉了MATLAB 的运行环境。初步掌握了MATLAB 语言在数字信号处理实验中一些基本库函数的调用和编写基本程序等应用，熟悉了差分方程求解的设计的一般步骤及处理后的结果，并能很好的理解数字信号处理中的基本概念，基本分析方法，使理论联系了实际，进一步认识到将所学的东西运用于实践的重要性。
另外，通过本次课程设计让我深深的了解到MATLAB 语言的重要性，它在以后的学习中还会用到，相信在以后学习的过程中，我会比较好的掌握并运用它。
八．参考文献
[1] 高西全，丁玉美. 数字信号处理[M].西安：西安电子科技大学出版社，2001
[2] 楼顺天，李博涵. 基于MATLAB 的系统分析与设计-信号处理[M].西安：西安电子科技大学出版社，1998
[3] 薛山.MATLAB 基础教程[M].北京：清华大学出版社，2011
[4] 吴大正. 信号与线性系统分析[M].北京：高等教育出版社，2005
[5] 刘敏，魏玲.MATLAB 通信仿真与应用[M]．北京：国防工业出版社，2001
[6] 赵红怡，张常年. 数字信号处理及其MATLAB 实现[M].北京：化学工业出版社，2002
[7] 何强，何英.MATLAB 扩展编程[M].北京：清华大学出版社，2002：293-296