Python3 脉冲响应和方差分解 VAR模型 脉冲响应分析stata

背景介绍

在实际工程中，由于系统的测量都是载噪的，而且噪声对观测数据的影响常常达到不可忽略的地步，因此当噪声影响足以使得要求的精度不足时，就必须考虑噪声的影响。实际中，系统噪声存在各种难以精确描述的因素：数学模型中未加以考虑的各种干扰作用、模型线性化及其他类似的假设所引起的误差等（系统误差）、输入输出的量测误差等（量测误差）等等，这些误差都难以用精准的数学模型来描述，且难以使用物理工具来测量。这些难以用确定性模型来精确描述的随机因素对过程的影响不容忽视，因此处理它们的方式应尽量简化。一般我们在确定性模型上以迭加输入的方式来考虑噪声的影响，并将多个噪声源综合成一个等效的噪声源来描述。
由于相关分析法可以提高对噪声的抗干扰性，伪随机二位式序列（PRBS，Pseudo Random Binary Signal）可以保证系统所有模态，且不会干扰系统正常运行，能够进行在线辨识，因此我们使用伪随机二位式信号与相关分析法相结合进行辨识系统的脉冲响应函数。
现需对三阶系统$G(s)=(b_0 s+b_1)/(a_0 s^3+a_1 s^2+a_2 s^1+a_3 )$使用M序列作为信号输入产生脉冲响应，利用Hankel矩阵法进行辨识，并将辨识结果与原系统进行比较分析，其中$b_0=1，b_1=15，a_0=1，a_1=5，a_2=4，a_3=15。$

相关分析法介绍

对于单输入单输出（SISO，Single Input Single Output）线性系统，若输入u(t)为一个平稳随机过程，则输出z(t)也为一个平稳随机过程，而平稳随机过程有两大数字特征：均值为常数与自相关函数为单变量函数（即不随时间变化）。
由相关理论可知，若输入u(t)为一个各态历经的平稳随机过程，则输出z(t)也为一个各态历经的平稳随机过程，而对于各态历经的平稳随机过程，它的集平均(均值、自相关函数等)，可以用一个样本的时间平均来代替。因此我们就可以将系统的自相关函数用过程的样本时间平均来代替：
$\left\{\begin{aligned} R_{u y}(\tau) &=R_{y u}(-\tau)=\lim _{T \rightarrow \infty} \frac{1}{T} \int_{0}^{T} y(t) u(t-\tau) d t \\ R_{u u}(\tau) &=R_{u u}(-\tau)=\lim _{T \rightarrow \infty} \frac{1}{T} \int_{0}^{T} u(t) u(t+\tau) d t \\ R_{u z}(\tau) &=R_{z u}(-\tau)=\lim _{T \rightarrow \infty} \frac{1}{T} \int_{0}^{T} z(t) u(t-\tau) d t \\ &=\lim _{T \rightarrow \infty} \frac{1}{T} \int_{0}^{T}[y(t)+n(t)] u(t-\tau) d t \\ &=R_{u y}(\tau)+\lim _{T \rightarrow \infty} \frac{1}{T} \int_{0}^{T} n(t) u(t-\tau) d t \end{aligned}\right.$

设干扰噪声n(t)与u(t)不相关，且n(t)为零均值，则：
$\begin{aligned} R_{u z}(\tau) &=R_{u y}(\tau)+\lim _{T \rightarrow \infty} \frac{1}{T} \int_{0}^{T} n(t) u(t-\tau) d t \\ &=R_{u y}(\tau)+R_{u n}(\tau)=R_{u y}(\tau) \end{aligned}$
上式表明，通过相关分析得到实测输出$z(t)$与输入$u(t)$之间的互相关函$R_{uz} (τ)$，在$n(t)$与$u(t)$不相关的条件下（此条件常常能满足)，等价于真实输出$y(t)$与输入$u(t)$间的互相关函数$R_{uy} (τ)$，因此，相关分析法有较强的抗干扰能力。下面我们就来求解$R_{uy}(τ)$。$\begin{aligned} R_{u y}(\tau) & \stackrel{\Delta}{=} E\{u(t) \cdot y(t+\tau)\}=E\left\{u(t) \int_{0}^{\infty} g(\lambda) u(t+\tau-\lambda) d \lambda\right\} \\ &=\int_{0}^{\infty} E\{u(t) g(\lambda) u(t+\tau-\lambda)\} d \lambda \\ &=\int_{0}^{\infty} g(\lambda) E\{u(t) u(t+\tau-\lambda)\} d \lambda \end{aligned}$
所以进一步可以写为：
$R_{u y}(\tau)=\int_{0}^{\infty} g(\lambda) R_{u u}(\tau-\lambda) d \lambda$
上式也被成为维纳霍普方程，它是相关分析法的理论基础。与古典法相比，相关分析法解决了抗干扰问题，但是又引出了积分难的问题，对此我们的解决方法是将采用白噪声作为输入，因为其均值为常数，即$E[u(t)]=0$，其自相关函数$R_{uu}(t)=Kδ(t)$，将其带入上式中可得：$R_{u z}(\tau)=\int_{0}^{\infty} g(\lambda) K \delta(-t) d t=K g(\tau)$所以：$\hat{g}(\tau)=\frac{1}{K} R_{u z}(\tau)$由此相关分析法也解决了古典法另一不能解决的问题：干扰系统正常工作，但是由于$\hat{g}(\tau)=\frac{1}{K} R_{u z}(\tau)$要求积分时间$T$无穷大，实际中不能做到，我们就采用具有周期性的、近似白噪声的伪随机信号作为输入信号，以解决积分时间长的问题。因为$u(t)=u(t+T)$，所以：$R_{u u}(\tau+T)=E[u(t) u(t+\tau+T)]=E[u(t) u(t+\tau)]=R_{u u}(\tau)$又因为：$\begin{aligned} R_{u u}(\tau) &=\lim _{T^{\prime} \rightarrow \infty} \frac{1}{T^{\prime}} \int_{0}^{T^{\prime}} u(t) u(t+\tau) d t \\ &=\lim _{k \rightarrow \infty} \frac{1}{k T} \int_{0}^{k T} u(t) u(t+\tau) d t \\ &=\lim _{k \rightarrow \infty} \frac{k}{k T} \int_{0}^{T} u(t) u(t+\tau) d t \\ &=\frac{1}{T} \int_{0}^{T} u(t) u(t+\tau) d t \end{aligned}$所以：$\begin{aligned} R_{u y}(\tau) &=\int_{0}^{\infty} g(\lambda) R_{u u}(\tau-\lambda) d \lambda \\ &=\int_{0}^{\infty} g(\lambda)\left[\frac{1}{T} \int_{0}^{T} u(t) u(t+\tau-\lambda) d t\right] d \lambda \\ &=\frac{1}{T} \int_{0}^{T} u(t)\left[\int_{0}^{\infty} g(\lambda) u(t+\tau-\lambda) d \lambda\right] d t \\ &=\frac{1}{T} \int_{0}^{T} u(t) y(t+\tau) d t \end{aligned}$进一步：$\begin{aligned} \mathrm{R}_{\mathrm{uu}}(\tau) &=\mathrm{K} \cdot \delta(\tau)=\mathrm{K} \cdot \delta(\tau+\mathrm{T}) \\ \mathrm{R}_{\mathrm{uy}}(\tau) &=\int_{0}^{\infty} \mathrm{g}(\lambda) \mathrm{R}_{\mathrm{uu}}(\tau-\lambda) \mathrm{d} \lambda \\ &=\int_{0}^{\mathrm{T}} \mathrm{g}(\lambda) \mathrm{R}_{\mathrm{uu}}(\tau-\lambda) \mathrm{d} \lambda+\int_{\mathrm{T}}^{2 \mathrm{~T}}+\int_{2 \mathrm{~T}}^{3 \mathrm{~T}}+\cdots \\ &=\int_{0}^{\mathrm{T}} \mathrm{g}(\lambda) \cdot \mathrm{K} \cdot \delta(\tau-\lambda) \mathrm{d} \lambda+\int_{\mathrm{T}}^{2 \mathrm{~T}} \mathrm{~g}(\lambda) \cdot \mathrm{K} \cdot \delta(\tau-\lambda+\mathrm{T}) \mathrm{d} \lambda+\int_{2 \mathrm{~T}}^{3 \mathrm{~T}}+\cdots \\ &=\mathrm{K}[\mathrm{g}(\tau)+\mathrm{g}(\tau+\mathrm{T})+\mathrm{g}(\tau+2 \mathrm{~T})+\cdots] \end{aligned}$

伪随机二位式序列

我们把接近于白噪声、序列中只含有+1、-1的一个周期序列叫做伪随机二位式序列（PRBS，Pseudo Random Binary Signal），它具有规律性、且可通过移位寄存器人为产生。PRBS移动若干位后，再和原来的PRBS对应位相加（模二加法），所得结果仍然是PRBS，只相当于原来的PRBS移动几位的效果，利用该性质人们可产生互不相关的伪随机信号。
对于一个$n$级的移位寄存器，它能形成序列的最大可能长度为$2n-1$个状态的组合，若一个$n$级的移位寄存器生成序列的周期长度$T$达到$T= 2n-1$，则称该序列为二位式最大长度序列，或称M序列，其具有如下性质：

1. 是一个确定的周期性序列，它的周期长度$N_p = 2^n -1$；
2. 一个周期内，“0”状态比“1”状态少1个，“1”状态有$2^{n-1}=(N_p+1)/2$ 个，“0”状态有 $2^{n-1}-1=(N_p-1)/2$
3. 若将序列中相邻状态不变的那一部分长度称“游程”(或“段”），则在一个周期内的游程总数为$m$，游程长度为$i$的有$m/2^i$ 个，游程长度为$n$的有1个；
4. 若将一个M序列与将其延迟了$r$个码以后的序列，按模2加法原则相加，所得的新序列还是M序列，不过延迟了$q$个码，$r、q$均为整数，且$r≥1，q≤ N_p-1$，即$u(k)⊕u(k+r)=u(k+q)$，我们也将这种性质称为移位相加性;
5. M序列具有近似离散的白噪声性质；

将$2N_p$个码的M序列{u(k)}与$2N_p$个码的方波信号{m(k)}，按模2加法规则相加，即可得到逆重复M序列{l(k)}，即$u(k)⊕m(k)=l(k)$，相较于M序列，逆重复M序列更接近白噪声，它具有以下性质：

1. 逆重复M序列{l(k)}的周期等于$2T_p$（即为原来M序列{u(k)}周期的两倍）为偶数；
2. 逆重复M序列的前、后半个周期是逆重复的，即$l(t)=-l(k+T_p)$；
3. 一个周期内，“0”状态与“1”状态的个数相等，各为$N_p$；
4. 逆重复M序列{l(k)}与M序列{u(k)}不相关，即：$R_{ul}(τ)=0$；
5. 逆重复M序列自相关函数R(t)与M序列自相关函数$R_uu (T)$的关系为: $R_{11}(μΔ)=(-1)^μ R_{uu}(μΔ)，μ=0,1,2,⋯,2N_P-1$；

M序列辨识线性系统的脉冲响应函数

我们采用多个周期进行计算$\hat{g}(u)$，首先我们定义$\hat{g}=\left[\begin{array}{c} \hat{g}(0) \\ \hat{g}(1) \\ \vdots \\ \bar{g}\left(N_{p}-1\right) \end{array}\right] \quad R_{v y}=\left[\begin{array}{c} R_{u y}(0) \\ R_{u y}(1) \\ \vdots \\ R_{u y}\left(N_{p}-1\right) \end{array}\right]$$\mathrm{R}_{\mathrm{uy}}(\mu)=\frac{1}{\mathrm{r} \cdot \mathrm{N}_{\mathrm{p}}} \sum_{\mathrm{k}=0}^{\mathrm{N} \mathrm{N}_{\mathrm{p}}-1} \mathrm{y}(\mathrm{k}) \cdot \mathrm{u}(\mathrm{k}-\mu) \quad \mu=0,1, \cdots, \mathrm{N}_{\mathrm{p}}-1$因此就可以得出：$\mathrm{R}_{\mathrm{uy}}=\frac{1}{\mathrm{r} \cdot \mathrm{N}_{\mathrm{p}}}\left[\begin{array}{cccc} \mathrm{u}(0) & \mathrm{u}(1) & \cdots & \mathrm{u}\left(\mathrm{rN}_{\mathrm{p}}-1\right) \\ \mathrm{u}(-1) & \mathrm{u}(0) & \cdots & \mathrm{u}\left(\mathrm{rN}_{\mathrm{p}}-2\right) \\ \vdots & \vdots & \cdots & \vdots \\ \mathrm{u}\left(-\mathrm{N}_{\mathrm{p}}+1\right) & \mathrm{u}\left(-\mathrm{N}_{\mathrm{p}}+2\right) & \cdots & \mathrm{u}\left(\mathrm{rN}_{\mathrm{p}}-\mathrm{N}_{\mathrm{p}}\right) \end{array}\right] \cdot\left[\begin{array}{c} \mathrm{y}(0) \\ \mathrm{y}(1) \\ \vdots \\ \mathrm{y}\left(\mathrm{rN}_{\mathrm{p}}-1\right) \end{array}\right]=\frac{1}{\mathrm{r} \cdot \mathrm{N}_{\mathrm{p}}} \mathrm{U} \cdot \mathrm{Y}$在结合相关分析法中介绍的内容可得：$\begin{aligned} \mathrm{g} &=\frac{\mathrm{N}_{\mathrm{p}}}{\mathrm{a}^{2} \Delta\left(\mathrm{N}_{\mathrm{p}}+1\right)}\left[\begin{array}{cccc} 2 & 1 & \cdots & 1 \\ 1 & 2 & \cdots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 1 & 1 & \cdots & 2 \end{array}\right] \mathrm{R}_{\mathrm{uy}} \\ =& \frac{1}{\mathrm{r} \cdot \mathrm{a}^{2} \Delta\left(\mathrm{N}_{\mathrm{p}}+1\right)}\left[\begin{array}{cccc} 2 & 1 & \cdots & 1 \\ 1 & 2 & \cdots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 1 & 1 & \cdots & 2 \end{array}\right] \mathrm{U} \cdot \mathrm{Y} \end{aligned}$这种基于一次完成法的计算$\hat{g}(u)$具备以下特点：

1. —次离线求出$\hat{g}(u)，μ=0,1,⋯,N_p-1$；
2. 需要输入$r+1$个周期的$u(k)：u(-N_p+1)~u(rN_p-1)$；
3. 精度要求较高时，$R_{uy}(u)$的计算精度要高，$r$的数目要大，所以数据存储量大；
4. 不是递推公式，无法在线辨识。
   同理，我们可以得出使用逆重复$M$序列辨识系统的脉冲响应函数，公式不做具体推导，如下式：
   $\mathrm{R}_{\mathrm{uy}}=\frac{1}{\mathrm{r} \cdot \mathrm{N}_{\mathrm{p}}}\left[\begin{array}{cccc} \mathrm{u}(0) & \mathrm{u}(1) & \cdots & \mathrm{u}\left(\mathrm{rN}_{\mathrm{p}}-1\right) \\ \mathrm{u}(-1) & \mathrm{u}(0) & \cdots & \mathrm{u}\left(\mathrm{rN}_{\mathrm{p}}-2\right) \\ \vdots & \vdots & \cdots & \vdots \\ \mathrm{u}\left(-\mathrm{N}_{\mathrm{p}}+1\right) & \mathrm{u}\left(-\mathrm{N}_{\mathrm{p}}+2\right) & \cdots & \mathrm{u}\left(\mathrm{rN}_{\mathrm{p}}-\mathrm{N}_{\mathrm{p}}\right) \end{array}\right] \cdot\left[\begin{array}{c} \mathrm{y}(0) \\ \mathrm{y}(1) \\ \vdots \\ \mathrm{y}\left(\mathrm{rN}_{\mathrm{p}}-1\right) \end{array}\right]=\frac{1}{\mathrm{r} \cdot \mathrm{N}_{\mathrm{p}}} \mathrm{U} \cdot \mathrm{Y}$$\hat{g}=\frac{N_{p}}{a^{2} \Delta\left(N_{p}+1\right)} R_{1 y}=\frac{1}{r \cdot a^{2} \Delta\left(N_{p}+1\right)} L \cdot Y$综上所述，我们可以列出PRBS辨识系统的步骤：
5. 预备试验，粗略了解系统的过渡过程时间$T_s$和系统工作频带: $f_M$（截止频率）、$f_m$
6. 选择PRBS信号的参数（$N_p，Δ，a$），要求信号的有效频带必须能完全覆盖系统的频带，即$Δ≤0.3 1/f_M$ 与$N_p≥1/(Δ⋅f_m )$，为了加宽频带，当$T_p$一定时，$Δ$下降，$N_p$上升；要求足够的信噪比，为了提高激励功率，当$T_p$一定时，$N_p$下降，$Δ$上升；要求$(N_p-1)Δ>T_s$；在达到最大信噪比的要求下，一般取既保证系统的线性，又不超出设备允许公差的最大幅值$a$；
7. 用计算机或专用仪器输出PRBS信号；
8. 给系统以预激励，在推导维纳—霍甫方程时，假设$u(t)$是平稳随机过程由于系统本身存在惯性，开始输入PRBS激励时，受非零初始条件的影响，使得$y(t)$为非平稳的；
9. 计算$R_{uy}(τ)$，在生产现场做试验，一般是在系统的正常工作状态$u_0$上再附加一个PRBS输入$u(t)$。

程序编写思路

根据第一部分对使用M序列辨识系统的脉冲响应函数的介绍，我们使用$M$序列与逆重复$M$序列作为信号输入，结合Hankel矩阵法对原始系统进行辨识，求出其传递函数，整体算法流程图如下图所示：

首先我们对系统进行初始化，分别设置$a、∆、n、r$，由此可以计算得到$N_p$，接着调用自定义函数产生$M$序列，利用公式1-11计算$R_{uy}$，利用前文公式计算系统脉冲响应估计量$\hat{g}$，由此构造Hankel矩阵$H(3,1)$，根据Hankel矩阵求解辨识系统的离散域的传递函数系数，由双线性变换法求解辨识系统连续域的传递函数，最后对原系统与辨识系统进行比较分析。

M序列及逆重复M序列的产生

我们需要使用M序列及逆重复M序列作为输入信号，由于多次需要使用，我们对其单独进行了自定义函数编写，然后在后续的程序中调用即可，产生M序列及逆重复M序列的算法流程图如下图所示：

首先我们需要输入$a、total、n、model$，其中$a$为PSBS序列的幅值，total等于PRBS长度的2倍，$n$代表$n$级移位寄存器，model代表选择输出的模式，若其为‘M’即输出M序列，否则输出逆M重复序列；接着采用ones函数进行寄存器初始化状态设置，将其状态都置为1；产生方波与$R(n)$进行异或运算，$R(n)$与$R(n-i)$进行异或运算（$n-i$为第$n-i$级状态），寄存器状态右移，最后当$k$迭代至total进行Model判断选择输出序列，否则将不断地迭代至满足条件$k=total$。我们最后选择Model判断虽然节省了代码空间，但是增大了内存运算量，实际中我们可以根据实际需求选择判断模式的位置。

这里需要特别注意，产生M序列时，不同$n$对应的不同寄存器，需要将不同的状态反馈进行异或运算，我们结合下图来说明：

由上图可知，我们需要确定的是如何选取第$n-i$级状态与第$n$级状态反馈进行异或运算反馈，这也是产生M序列的难点所在，所幸根据查阅资料可以得如下表所示的反馈系数状态表：

级数n	反馈系数i	第n-i状态
3	1	n-1
4	1	n-1
5	2	n-2
6	1	n-1
7	1	n-1
8	2	n-2
9	3	n-3
10	4	n-4
11	2	n-2

这里的$i$代表从第$n$级开始，从后往前数第一个可以能够与第$n$级状态进行异或运算的状态，后续产生不同级的寄存器需要根据此表不断改变第$n-i$级状态。
接下来我们对系统进行仿真实验。

仿真实验

预备实验

我们首先需要选择M序列的参数$n$与$∆$，而调节时间$T_s$和截止频率$f_M$是选择M序列参数的依据，利用step函数绘制系统的单位阶跃响应，利用margin函数绘制系统的幅频特性曲线，如下图所示：

原系统阶跃响应及幅频特性图

可以求得调节时间$T_s=30.7s，f_M=2.121rad/s=0.3376Hz$，这里需要注意Matlab默认的调节时间$T_s$为系统达到终值的2%误差带内所需的最短时间，我们需要手动调为达到系统终值5%误差带内的最短时间。

根据$\left\{\begin{array}{c} \frac{1}{3 \Delta} \geq f_{M} \\ \left(N_{p}-1\right) \Delta>T_{s} \end{array}\right.$,我们可以求得$\left\{\begin{array}{c} \Delta \leq 0.98730 s \\ N_{p}>32.1 \end{array}\right.$，所以$n>5.05$。至此我们大致确定了M序列的参数的范围，后续实验中需要通过不断地调试得到最优的辨识效果。

PRBS序列的产生

调用自定义的函数产生一个M序列与逆M序列，如下图所示：

M重复序列的产生图

逆M重复序列的产生图

上图是选择$a=2、total=126、n=6$的情况下产生的PRBS，从图中可以看出逆M重复序列的周期为M序列的两倍，并且逆M重复序列的前后半个周期是逆重复的，符合我们课本所学的知识。后续我们需要将产生的M序列或者逆M序列作为信号输入。

脉冲响应估计量

在产生M序列或逆M重复序列之后就可以将其作为信号输入，利用lsim函数产生在不同时刻下系统的输出值，如图3-3所示：

M序列作为输入的系统输出图

逆M序列作为输入的系统输出图上图是选择$a=2、total=126、n=6$的情况下产生的PRBS作为输入信号，我们选择0~12s的时间内作为输出显示。在计算出系统的输出值之后，我们就可以算系统脉冲响应的估计值，由于在6级寄存器产生的M序列辨识效果不明显，在经过不断地调试之后，我们得到了11级寄存器产生M序列作为输入信号辨识效果最佳，而∆的选择也会关系到系统辨识的效果，如下图是M序列作为输入信号不同的$∆$的选择产生的脉冲响应：

∆=0.5脉冲响应估计量图

∆=0.2脉冲响应估计量图

可以从上图中明显看出，$∆=0.5$时，脉冲响应的真实值与估计值不完全重合，$∆=0.2$时，脉冲响应的真实值与估计值几乎完全重合。对脉冲响应的真实值与脉冲响应估计值进行均方误差的计算，如下所示：$\operatorname{MSE}(\hat{g})=E(\hat{g}-g)^{2}$在不同的∆下，取不同时刻下脉冲响应的真实值与脉冲响应的估计值，计算其均方误差如下表所示：

∆	$a=2、total=126、n=6$
0.8	0.0643399
0.65	0.0446066
0.5	0.0270118
0.35	0.0136712
0.2	0.0045212

由此可见∆越小，脉冲响应的估计值与真实值之间的均方误差越小。

与M序列相同，我们只需改变系数矩阵即可得到逆M序列的辨识脉冲响应结果，如下图所示：

∆=0.5脉冲响应估计量图

∆=0.2脉冲响应估计量图

计算在不同的∆下，取不同时刻下脉冲响应的真实值与脉冲响应的估计值的均方误差，如下表所示：

∆
0.8	0.0643387
0.65	0.0446061
0.5	0.0270116
0.35	0.0136714
0.2	0.0045212

由上表可知∆越小，脉冲响应的估计值与真实值之间的均方误差越小，得到的结论与M序列作为输入是一致的。

Matlab代码

1.	function Out = Ssequence(n,a,total,model)
2.	% 选择产生逆M序列或者逆M序列
3.	% n级寄存器，a为M序列的幅值，del为时钟脉冲
4.	% total为M序列的长度,model为M序列或者逆M序列
5.	% 初始化n级寄存器
6.	R = ones(1,n);%初始n级寄存器各状态均为1
7.	m = 1;
8.	
9.	for k = 1 : total
10.	    temp1 = xor(m,R(n));%进行异或运算,产生逆M序列
11.	    if temp1 == 0
12.	        l(k) = -a;
13.	    else
14.	        l(k) = a;
15.	    end
16.	    m = not(m); %产生方波
17.	    temp2 = xor(R(n-2),R(n));
18.	    R = circshift(R,1);
19.	    R(1) = double(temp2);
20.	    if temp2 == 0
21.	        u(k) = -a;
22.	    else
23.	        u(k) = a;
24.	    end
25.	end
26.	
27.	if model == 'M'
28.	    Out = u;
29.	else
30.	    Out = l;
31.	end
32.	
33.	End
34.	
35.	% M序列作为输入辨识系统脉冲响应
36.	clc;
37.	close all;
38.	clear all;
39.	
40.	%% 产生M序列
41.	n = 11;%n级寄存器
42.	a = 2;%M序列的幅值
43.	del = 0.8;%del为时钟脉冲
44.	Np = 2^n - 1;%M序列的循环长度
45.	Total = 2 * Np;%逆M序列的长度
46.	r = 1;
47.	u = Ssequence(n,a,Total,'L');% 产生M序列或逆M序列
48.	% 产生输出响应y
49.	G = tf([1,15],[1,5,4,15]);
50.	TF = Total * del;
51.	tim = 0 : del : TF - del;
52.	y = lsim(G,u,tim);% G系统对Out输入的时间响应
53.	
54.	figure(1);
55.	stairs(tim,u);%绘制阶梯状图
56.	axis([0,12,-2.5,2.5]);
57.	hold on
58.	plot(tim,y,'r');
59.	hold off
60.	grid on
61.	
62.	%% 计算Ruy并得到系统脉冲响应的估计值
63.	C = 1 / (r * a^2 * del) / (Np + 1);
64.	C_mat = ones(Np,Np) + diag(ones(1,Np));
65.	
66.	j = 1;
67.	U = zeros(Np,Np*r);
68.	for i = 1 : Np
69.	    U(i,:) = u(1,Np+j:2*Np+j-1);
70.	    j = j - 1;
71.	end
72.	
73.	% ghat = C .* C_mat * U * y(Np+1:2*Np);% M序列作为信号输入系统脉冲响应估计值
74.	ghat = C * U * y(Np+1:2*Np);% 逆M序列作为信号输入系统脉冲响应估计值
75.	Result = [ghat y(Np+1:2*Np)];
76.	
77.	figure(2)
78.	impulse(G,50);
79.	grid on
80.	hold on
81.	t = 0 : del : 50;
82.	plot(t',ghat(1:length(t)),'r','LineWidth',1);
83.	legend('脉冲响应真实值','脉冲响应估计值')
84.	
85.	%% 计算均方误差
86.	% t = 0 : T0 : 60;
87.	y1 = impulse(G,t);
88.	y2 = ghat(1:length(t));
89.	delta_y1 = mean((y2 - y1).^2);
90.	
91.	%% 构造Hankel矩阵并计算辨识系统的传递函数
92.	H = [y2(2),y2(3),y2(4);y2(3),y2(4),y2(5);y2(4),y2(5),y2(6)];
93.	T0 = 0.8;%采样周期
94.	if det(H)<=0.0001
95.	    disp('Hankel矩阵奇异，无法求逆');
96.	else
97.	    A = inv(H)*[-y2(5);-y2(6);-y2(7)];
98.	    B = [1 0 0;A(3) 1 0;A(2) A(3) 1]*[y2(2);y2(3);y2(4)];
99.	    numd = B';
100.	    dend = [1 A(3) A(2) A(1)];
101.	    bssysd = tf(numd,dend,T0); % 创建1个采样时间为T0的离散时间传递函数
102.	end
103.	% bianshisys1 = d2c(bssysd,'tustin');%辨识出的传递函数
104.	bianshisys1 = d2c(bssysd);
105.	%% 绘制原系统与辨识系统的单位阶跃响应
106.	figure(3)
107.	step(T0*bianshisys1,'r')
108.	hold on 
109.	step(G,'g')
110.	grid on
111.	legend('辨识传递函数','实际传递函数');
112.	%% 计算均方误差
113.	t1 = 0 : T0 : 60;
114.	y3 = step(G,t1);
115.	y4 = step(T0*bianshisys1,t1);
116.	delta_y2 = mean((y3 - y4).^2)
117.	
118.	%% 乘同余法产生0-1均匀分布的随机数
119.	A = 5^13; 
120.	x0 = 1; 
121.	M = 10^36;
122.	N = 100;
123.	for k = 1 : N
124.	    x2 = A * x0;
125.	    x1 = mod(x2,M);
126.	    v1 = x1 / M;
127.	    v(:,k) = v1;
128.	    x0 = x1;
129.	end
130.	plot(1:N,v,'r');
131.	xlabel('k');
132.	ylabel('噪声幅值');
133.	title('白噪声序列');
134.	
135.	%% 混合同余法产生0-1均匀分布的随机数
136.	N = 7;
137.	M = 2^35;
138.	c = (0.5 + sqrt(3)/6) * M;
139.	v = [];
140.	for k = 1 : N
141.	    x2 = A * x0 + c;
142.	    x1 = mod(x2,M);
143.	    v1 = x1 / M;
144.	    v(:,k) = v1;
145.	    x0 = x1;
146.	end
147.	plot(1:N,v,'r');
148.	xlabel('k');
149.	ylabel('噪声幅值');
150.	title('白噪声序列');
151.	
152.	%% 白噪声的randn的产生及有色噪声序列产生
153.	L = 500;%仿真长度
154.	d = [1 -1.5 0.7 0.1];%有色噪声分母系数 
155.	c = [1 0.5 0.2];%有色噪声分子系数(可用roots命令求其根)
156.	nd = length(d)-1;%有色噪声分母阶次 
157.	nc = length(c)-1;%有色噪声分子阶次
158.	xik = zeros(nc,1); %白噪声初值，相当于4(k-1)(k-nc)
159.	ek = zeros(nd,1);%有色噪声初值
160.	xi = randn(L,1); %randn产生均值为0，方差为1的高斯随机序列(白噪声序列)
161.	for k = 1 : L
162.	    e(k) = -d(2:nd+1) * ek + c * [xi(k);xik]; %产生有色噪声%数据更新
163.	    for i = nd : -1 : 2
164.	        ek(i) = ek(i-1);
165.	    end
166.	    ek(1) = e(k);
167.	    for i = nc : -1 : 2 
168.	        xik(i) = xik(i-1);
169.	    end
170.	    xik(1) = xi(k);
171.	end
172.	subplot(2,1,1);
173.	plot(xi);
174.	xlabel('k');
175.	ylabel('噪声幅值');
176.	title('白噪声序列');
177.	subplot(2,1,2);
178.	plot(e);
179.	xlabel('k');
180.	ylabel('噪声幅值');
181.	title('有色噪声序列');
182.	
183.	G = tf([1,15],[1,5,4,15]);
184.	% bianshisys1 = tf([1,15],[1,5,4,15]);
185.	% bianshisys2 = tf([0.4,2.4],[1,2,0.64,0.96]);
186.	% bianshisys3 = tf([2,60],[1,10,16,120]);
187.	% bianshisys4 = tf([1,15],[1,5,4,15]);
188.	
189.	bianshisys1 = tf([0.0002232,0.996,15.09],[1,5.039,3.974,15.14]);
190.	bianshisys2 = tf([0.0002232,0.3984,2.414],[1,2.016,0.6359,0.969]);
191.	bianshisys3 = tf([0.0002232,1.992,60.34],[1,10.08,15.9,121.1]);
192.	bianshisys4 = tf([0.01851,0.9056,15.75],[1,5.241,4.045,15.75]);
193.	
194.	subplot(2,2,1)
195.	step(bianshisys1,'r')
196.	hold on 
197.	step(G,'g')
198.	grid on
199.	title('\Delta=0.2,T_{0}=0.2');
200.	legend('辨识传递函数','实际传递函数');
201.	
202.	subplot(2,2,2)
203.	step(bianshisys2,'r')
204.	hold on 
205.	step(G,'g')
206.	grid on
207.	title('\Delta=0.2,T_{0}=0.5');
208.	legend('辨识传递函数','实际传递函数');
209.	
210.	subplot(2,2,3)
211.	step(bianshisys3,'r')
212.	hold on 
213.	step(G,'g')
214.	grid on
215.	title('\Delta=0.2,T_{0}=0.1');
216.	legend('辨识传递函数','实际传递函数');
217.	
218.	subplot(2,2,4)
219.	step(bianshisys4,'r')
220.	hold on 
221.	step(G,'g')
222.	grid on
223.	title('\Delta=0.5,T_{0}=0.5');
224.	legend('辨识传递函数','实际传递函数');