2.差分方程及其求解---数字信号处理实验报告

计算机与信息工程学院验证性实验报告

一、实验目的

1．学习并掌握系统的差分方程表示方法以及差分方程的相关概念。 2．熟练使用filter 函数对差分方程进行数值求解。 3．掌握差分方程的求解及MATLAB 实现方法。

二、实验原理及方法

1．一LTI 系统可以用一个线性常系数差分方程表示：

()()N M

k m

k m a y n k b

x n m ==-=-∑∑，任意n

如果N a ≠0，那么这个差分方程就是N 阶的，已知系统的输入序列，用这个方程可以根据当前输入x(n)和以前M 点的输入x(n-m ),…,x(n-1),以及以前N 点的输出y(n-N),…,y(n -1)来计算当前输出y(n)。在实际中这个方程在时间上是从n =-∞到n =+∞朝前计算的，因此该方程的另一种形式是：

()()()M

N

m k m k y n b x n m a y n k ===---∑∑

方程的解能以下面形式求得：()()()H p y n y n y n =+分别为方程的齐次解跟特解部分。已知输入和差分方程的稀疏，可用filter 对差分方程进行数值求解。最简单形式为：y=filter(b,a,x) 其中b=[b0,b1,…,bM];a=[a0,a1,…,aN];

2．上面差分方程解的形式为齐次解和特解，另外还可以求零输入解和零状态解理论计算中要用到z 变换，请好好掌握z 变换的内容。用MATLAB 实现时，若已知初始条件，则应用y=filter(b,a, x, xic)来求完全响应。这里xic 是初始状态输入数组。MATLAB 还提供一种filtic 函数来得到xic 。

xic=filtic(b,a,Y,X )其中b 和a 是滤波器系数数组，Y 和X 是分别从y(n)和x(n)的初始条件得来的初始状态数组。 三.实验内容

1、已知下面差分方程：

()(1)0.9(2)()y n y n y n x n --+-=;任意n

要求：

a.计算并画出在n=-20,…,100的脉冲响应h(n);

b.计算并画出在n=-20,…,100的单位阶跃相应s(n). 2.解以下差分方程：

31

()(1)(2)(),022

y n y n y n x n n -

-+-=≥ 其中1

()()()4

n x n u n =，初始条件为y(-1)=4和y(-2)=10。

要求先用理论计算，再用MATLAB 编程实现，并对比两个结果。 参考流程图： 实验内容1．

实验内容2.

四、实验报告要求

1．总结差分方程的性质及应用．

2．写出实验程序．记录实验数据并与理论计算作比较，总结结果。

五、实验程序及结果

针对问题一：

问题a的程序如下：

先建立一个生成单位抽样序列的M文件

function [x,n]=impseq(n0,n1,n2)

if((n0n2)|(n1>n2))

error('参数必须满足n1<=n0<=n2')

end

n=n1:n2;

x=[(n-n0)==0];

主程序如下：

clc;clear all

b=1;a=[1 -1 0.9];

[x,n]=impseq(0,-20,100);

y=filter(b,a,x);

plot(n,y)

title('输出脉冲响应h(n)')

xlabel('n')

-20

20

4060

80

100

-1-0.8-0.6-0.4-0.200.20.40.60.8

1输出脉冲响应h(n)

n

问题b 的程序如下：

先建立一个生成单位阶跃序列的M 文件

function [x,n]=stepseq(n0,n1,n2) if((n0n2)|(n1>n2))

error('参数必须满足n1<=n0<=n2') end

n=[n1:n2];

x=[(n-n0)>=0]; 主程序如下： clc;clear all b=1;a=[1 -1 0.9];

[x,n]= stepseq(0,-20,100); y=filter(b,a,x); plot(n,y)

title('输出单位阶跃响应s(n)') xlabel('n')

-20

020

406080100

00.5

1

1.5

2

2.5

输出单位阶跃响应s(n)

n

针对问题二： 1、理论计算 系统函数为

12

()1

()()1 1.50.5Y z H z X z z z --=

=-+ 输入信号的Z 变换为

1

1

()10.25X z z

-=

- 对给出的差分方程求单边Z 变换，得：

112() 1.5[(1)()]0.5[y(1)(2)()]()Y z y z Y z z y z Y z X z ---=-+--+-++

1

1212

112()()1 1.50.51 1.50.5z Y z X z z z z z ------=+-+-+

111211(12)(10.25)()(1 1.50.5)(10.25)z z Y z z z z -----+--=

-+-

111

6117

110.510.25z z z ---=

-+

--- 故

()6()110.5()70.25()n n y n u n u n u n =-?+?

把n=1，2，3…11代入上式得：y(1)=2, y(2)=1.25, y(3)= 0.9375, y(4)= 0.7969,

y(5)= 0.7305, y(6)= 0.6982,y(7)= 0.6824, y(8)= 0.6745, y(9)= 0.6706, y(10)= 0.6686, y(11)= 0.6676

2、MATLAB编程实现

clc;clear all

b=[1];a=[1,-1.5,0.5]; %系统函数分子分母系数

ys=[4,10]; %初始条件

xic=filtic(b,a,ys); %初始条件形成的等效输入序列

bx1=[1];ax1=[1,-0.25]; %X(z）的系数

ay1=conv(a,ax1); %Y(z)的分母系数

by1=conv(b,bx1)+conv(xic,ax1); %Y(z)的分子系数

[r,p,c]=residuez(by1,ay1) %部分分式展开法

n=[0:10];x=0.25.^n;

y=filter(b,a,x,xic)

plot(n,y)

title('输出序列y(n)')

xlabel('n')

hold on

plot(n,y,'*')

n

由上结果可知：理论计算和MATLAB编程所算得的结果是完全一样的。

六、实验小节

经过本次实验掌握了系统的差分方程表示方法以及差分方程的相关概念。并熟练了使用filter函数对差分方程进行数值求解的方法。更进一步懂得了用z 变换来求零输入解和零状态解的理论知识。

教师签名：

年月日