常系数线性差分方程解法

常系数线性差分方程的一般形式为
$$\sum _{k=0}^N{a}_{k}y(n-k)=\sum _{r=0}^M{b}_{r}x(n-r)\dots \dots （1）$$

1. 解方程法

1.1 概念

离散系统的完全响应，可分解为自由响应与强迫响应之和，或分解为零输入响应与零状态响应之和。
零输入响应y_{zi}：当x(n)=0（即方程右端为零）时，系统的输出函数，对应方程（1）的齐次解。
零状态响应y_{zs}：当接入激励前（n<0），系统的输出状态均为0时，系统的输出函数，对应y(-1)，y(-2)，……，y(-N)均为0时（系统满足零状态条件），方程（1）的通解。
自由响应：自由响应对应于方程（1）的齐次解，意为没有外部输入时系统的输出。
强迫响应：对应于方程（1）的特解，对应于加上外部输入后系统的输出。
当系统满足零状态条件时，零状态响应为自由响应和强迫响应的叠加即方程（1）的通解。
一般情况下，零输入响应是自由响应的一部分，零状态响应包括自由响应的另一部分和强迫响应
（那为啥自由响应不等于零输入响应？？？）

1.2 求解步骤

求解分三步进行
①求齐次方程的通解
方程（1）对应的齐次方程为
$$\sum _{k=0}^N{a}_{k}y(n-k)=0\dots \dots （2）$$
设y(n)=${\alpha }^n$，则$\alpha$为一元N次方程$\sum _{k=0}^N{a}_k{\alpha }^{n-k}=0$的根。
若方程没有重根，则N个根分别为${\alpha }_1，{\alpha }_2，\dots \dots ，{\alpha }_N$，则齐次方程的通解为
$$y_1(n)=\sum _{i=1}^N{c}_i{\alpha }_i^n=c_1{\alpha }_1^n+c_2{\alpha }_2^n+\dots \dots +c_N{\alpha }_N^n\dots \dots (3)$$
若方程有1个m重根${\alpha }_1$，其余为非重根${\alpha }_2，\dots \dots ，{\alpha }_{N-m+1}$，则齐次方程通解为
$$y_1(n)=c_{11}{\alpha }_1^n+c_{12}n{\alpha }_2^n+\dots \dots +c_{1m}{n}^{m-1}{\alpha }_1^n+c_2{\alpha }_2^n+c_3{\alpha }_3^n+\dots \dots +c_{N-m+1}{\alpha }_{N-m+1}^n\dots \dots (4)$$
其中任意常数c_i、c_{ir}由系统的初始状态（方程的边界条件）确定。

②确定差分方程（1）的一个特解D(n)
求出右端具体表达式 -> 根据结果选择含待定系数的特解 -> 带入方程确定系数

输入x[n]	特解D(n)
常数	常数
$cos(\Omega n+\varphi )$或$sin(\Omega n+\varphi )$	$c_{1}cos(\Omega n)+c_{2}sin(\Omega n)$
$a^n$	$\sum _{i=0}^k{c}_i{n}^i{a}^n$（a是方程的k重特征根）
$n^k$	$\sum _{i=0}^k{c}_i{n}^i$
$n^k{a}^{n}cos(\Omega n)$或$n^k{a}^{n}sin(\Omega n)$	$(B_1{n}^k+\dots \dots +B_{k}n+B_{k+1})a^{n}cos(\Omega n)+(D_1{n}^k+\dots \dots +D_{k}n+D_{k+1})a^{n}sin(\Omega n)$

③确定完全解中的待定系数

2. 迭代法

将方程（1）进行移项
$$y(n)=-\sum _{k=0}^N\frac{a_k}{a_0}y(n-k)+\sum _{r=0}^M\frac{b_r}{a_0}x(n-r)\dots \dots (5)$$

import numpy as np

def solveCCLDE(n,x,y0,a,b):
    """迭代法解常系数线性差分方程
    N+1=len(a),M+1=len(b)
    :param n:返回输出序列的长度
    :param x:输入x(n)，从x(0)到x(n-1)
    :param y0:初始状态，从y(-N)到y(-1)
    :param a:输出项系数
    :param b:输入项系数
    :return y:长度为n的输出序列,从y(0)到y(n-1)
    """
    N=len(a)-1
    M=len(b)-1
    #y[-N]~y[-1]对应x[-M]~x[0]，补充x[-M]~x[-1]
    x=np.array([0]*M+x)
    a=np.mat(a).transpose()
    b=np.mat(b).transpose()
    b=b/a[-1,0]
    a=a/a[-1,0]
    
    y=[]
    for _ in range(n):
        xx=np.mat(x[_:_+M+1])
        yy=np.mat(y0)
        y.append((-yy*a[:-1]+xx*b).tolist()[0][0])
        y0.append(y[-1])
        y0=y0[1:]
        
    return y

3. z变换法

利用单边z变换的移位性质
$$\mathcal{Z}[x(n-m)u(n)]=z^{-m}[X(z)+\sum _{k=-m}^{-1}x(k)z^{-k}]\dots \dots （6）$$
对方程(1）两边取单边z变换得
$$\sum _{k=0}^N{a}_k{z}^{-k}[Y(z)+\sum _{l=-k}^{-1}y(l)z^{-l}]=\sum _{r=0}^M{b}_r{z}^{-r}[X(z)+\sum _{m=-r}^{-1}x(m)z^{-m}]\dots \dots （7）$$
因为输入序列x(n)为因果序列，则$\sum _{m=-r}^{-1}x(m)z^{-m}=0$，方程变为
$$\sum _{k=0}^N{a}_k{z}^{-k}[Y(z)+\sum _{l=-k}^{-1}y(l)z^{-l}]=\sum _{r=0}^M{b}_r{z}^{-r}X(z)\dots \dots （8）$$
解出Y(z)得
$$Y(z)=\frac{\sum _{r=0}^M{b}_r{z}^{-r}}{\sum _{k=0}^N{a}_k{z}^{-k}}X(z)-\frac{\sum _{k=0}^N[a_k{z}^{-k}\sum _{l=-k}^{-1}y(l)z^{-l}]}{\sum _{k=0}^N{a}_k{z}^{-k}}\dots \dots （9）$$
（9）式中第一项仅与输入有关，与起始状态无关，是零状态响应的z变换，记为$Y_{zs}(z)$，第二项仅与起始状态有关，与输入无关，是零输入响应的z变换，记为$Y_{zi}(z)$，令
$$H(z)=\frac{\sum _{r=0}^M{b}_r{z}^{-r}}{\sum _{k=0}^N{a}_k{z}^{-k}}\dots \dots （10）$$
H(z)完全又系统特性所决定，称为离散时间系统的“系统函数”。