求解系统脉冲响应

前言

用伪随机信号求解系统的脉冲响应非常有意思，是我在反馈系统设计这门课听到的一种方法。最近终于有空可以学习一下相关的知识了。本人在系统相关分析方面还是小白，仅此记录下学习过程。

相关分析与系统

（一）脉冲输入

对于单输入单输出的线性系统，形如下图¹ ：

输入输出关系为：
$$y(t)={\int }_0^{\infty }x(t-\tau )g(\tau )d\tau$$

如果要得到$g(t)$,，那么只需要输入脉冲函数（狄拉克函数）$x(t)=\delta (t)$，则由狄拉克函数的筛选性：

$$y(t)={\int }_0^{\infty }\delta (t-\tau )g(\tau )d\tau =g(t)$$

这样貌似只要我们能得到理想的脉冲，我们就能得到脉冲响应。然而这只是理想情况，理想脉冲在物理上是难以实现的，而且理想脉冲的能量都集中在零点，因此这种方式求的脉冲响应很容易受到非零时刻的干扰。

（二）白噪声输入

结合上面的系统，如果不妨设$x(t)$的自相关函数为$R_{xx}(t)$, $x(t),y(t)$的互相关函数为$R_{xy}(t)$，那么，由维纳——霍普夫方程 可得：

$$R_{xy}(t)={\int }_0^{\infty }g(\tau )R_{xx}(t-\tau )d\tau$$

这个方程非常难以应用，但是如果输入是白噪声就不同了，白噪声的自相关函数为$R_{xx}(t)={\sigma }^2\delta (t)$，代入上式发现：
$$R_{xy}(t)={\sigma }^{2}g(t)$$

输入白噪声可以通过求输出的互相关函数得到脉冲响应。然而，理想的白噪声也是不容易实现的，但它在很宽的时域内都作用于系统，且具有一定的能量，因此它具有一定的抗干扰能力。

（三）m序列输入

m序列同白噪声都是随机信号，但m序列是可人工生成的周期性的伪随机信号。m序列是指最长线性移位寄存器序列，生成方式如下² :

n个移位寄存器生成的最长周期序列为$2^n-1$，但并不是所有的反馈系数$c_1,c_2,...,c_n$都可以产生最长的周期序列，而符合最长的周期序列的输出就是m序列。就不同的$n$所对应的反馈系数需要用到本原多项式，本原多项式各阶次的系数（除0阶外）作为各以为寄存器的反馈系数。对于不同的$n$，其对应的所有本原多项式可用如下代码求得：

>> n = 5;  %假如n为5
>> primpoly(n, 'all');

得到反馈系数后，可以使用如下代码生成一个周期的m序列³：

function [seq]=mseq(coef)
%***************************************************
% 此函数用来生成m序列
% coef为反馈系数向量
% Author: FastestSnail
% Date: 2017-10-03
%***************************************************
m=length(coef);
len=2^m-1; % 得到最终生成的m序列的长度     
backQ=0; % 对应寄存器运算后的值，放在第一个寄存器
seq=zeros(1,len); % 给生成的m序列预分配
registers = [1 zeros(1, m-2) 1]; % 给寄存器分配初始结果
for i=1:len
    seq(i)=registers(m);
    backQ = mod(sum(coef.*registers) , 2); %特定寄存器的值进行异或运算，即相加后模2
    registers(2:length(registers)) = registers(1:length(registers)-1); % 移位
    registers(1)=backQ; % 把异或的值放在第一个寄存器的位置
end
end

（四）矩阵法

对周期性的m序列的输入，设周期为$T$，其互相关函数为：
$$R_{xy}(t)=\frac{1}{T}{\int }_0^T{R}_{xx}(t-\tau )g(\tau )d\tau$$
因此，对于周期性伪随机信号，只需要在一个周期内计算即可。

我们不妨使用离散方式逼近连续的系统。假设采样时间为$t_s$，m序列电平幅值为$\pm a$，每个电平持续时间和采样周期相同，那么m序列的自相关序列为¹：
$$R_{xx}(t)=\{\begin{array}{ll}{a}^2,& t=0\\ -\frac{a^2}{N},& t=t_s,2t_s,...(N-1)t_s\end{array}$$

设$\tau =\mu t_s$为个信号周期内各个采样点时刻，$\mu =0,1,2,...,N-1$。并且直接把$t_s$用单位1代替，带入离散形式的维纳——霍普夫方程：

$$R_{xy}(\tau )=\frac{1}{N}\sum _{k=0}^{N-1}g(k)R_{xx}(k-\tau )△$$

为了能计算上式，必须至少观察2个信号周期，然后用后一个周期进行计算。

写成矩阵形式：$R_{xy}=△Rg$
其中：
$$R_{xy}=\left[\begin{array}{c}{R}_{xy}(0)\\ R_{xy}(1)\\ ⋮\\ R_{xy}(N-1)\end{array}\right]g=\left[\begin{array}{c}g(0)\\ g(1)\\ ⋮\\ g(N-1)\end{array}\right]$$

$$\begin{array}{rl}R& =\left[\begin{array}{cccc}{R}_{xx}(0)& R_{xx}(1)& \dots & R_{xx}(N-1)\\ R_{xx}(-1)& R_{xx}(0)& \dots & R_{xx}(N-2)\\ ⋮& ⋮& \ddots & ⋮\\ R_{xx}(-N+1)& R_{xx}(-N+2)& \dots & R_{xx}(0)\end{array}\right]\\ & =a^2\left[\begin{array}{cccc}1& -1/N& \dots & -1/N\\ -1/N& 1& \dots & -1/N\\ ⋮& ⋮& \ddots & ⋮\\ -1/N& -1/N& \dots & 1\end{array}\right]\end{array}$$

这时得到$g=\frac{1}{△}{R}^{-1}{R}_{xy}$。
而根据互相关函数$R_{xy}$的定义可求：
$$R_{xy}(\tau )=\frac{1}{N}\sum _{k=0}^{N-1}x(k-\tau )y(k)$$
写成矩阵形式：
$$R_{xy}=\left[\begin{array}{cccc}x(0)& x(1)& \dots & x(N-1)\\ x(-1)& x(0)& \dots & x(N-2)\\ ⋮& ⋮& \ddots & ⋮\\ x(-N+1)& x(-N+2)& \dots & x(0)\end{array}\right]\left[\begin{array}{c}y(0)\\ y(1)\\ ⋮\\ y(N-1)\end{array}\right]=XY$$

所以最终得到¹：
$$g=\frac{1}{a^2(N+1)t_s}\left[\begin{array}{cccc}2& 1& \dots & 1\\ 1& 2& \dots & 1\\ ⋮& ⋮& \ddots & ⋮\\ 1& 1& \dots & 2\end{array}\right]XY$$

注意这里$g$的系数，分子是1，而不是ppt中的$N$

（五）迭代法

迭代法在引用一¹ 中有详细记述(公式过多)，优点是内存消耗小，但经尝试原理上和精度上未和矩阵法有较大不同。下面只给出尝试的代码。

（六）参数选择

• m序列电平幅值，经过尝试，可能不受影响，准确的说对采样周期和电平持续时间相等的情况不受影响。但引用一¹ 上提到$a$ 取大一点为好。
• 经尝试，$t_s$越小越好，但可能实际情况并不允许很小的$t_s$。选择方式参考了另一资料⁴ ，$1/(3t_s)\ge f_{max}$，$f_{max}$为最大系统工作频率，可以考虑为截止频率。
• m序列周期长度$N$，最好选择为满足$N{t}_s\ge T_s$，$T_s为过渡过程时间，在这里就是脉冲响应衰减为0所需时间。

举例说明

（一）问题

尝试绘制典型二阶系统的脉冲响应的波形：
$$G(s)=\frac{w_n^2}{s^2+2\zeta w_{n}s+w_n^2}$$
其中，$w_n=10,\zeta =0.5$

（二）参数选择

1.计算传递函数的截止频率
一般认为，传递函数的幅频曲线降低到-3dB时对应的频率为截止频率$w_c$，经计算得：$w_c=12.7rad/s\approx 2Hz$。
2.确定脉冲响应过渡过程时间
可对系统输入一方波，观察输出的过渡时间，但这里为了方便直接使用Matlab施加理想脉冲：

观察得过渡过程时间为1s。

• 根据$1/(3t_s)\ge f_{max}$， 可推算$t_s\le 0.167s$，下面取$t_s=0.1s,0.01s$进行试验。
• 根据$N{t}_s\ge T_s$ ，可推算$N>10,100$，下面取$N=15,127$进行试验。

（三）代码

• 矩阵法

clc
clear
% 构造已知传递函数
s = tf('s');
wn = 10;
damp = 0.5;
Gs = wn^2/(s^2 + 2*damp*wn*s + wn^2);
%离散化
Ts = 0.1;
Gz = c2d(Gs,Ts);

% 差分方程系数
a = Gz.Numerator{1}(2);
b = Gz.Numerator{1}(3);
c = Gz.Denominator{1}(2);
d = Gz.Denominator{1}(3);

% 闭环系统，最大工作频率如果取闭环截止频率，2Hz
% 那么不妨取dt=0.1s，N=15。
%构造m序列
coef = [0 0 1 1];
seq = mseq(coef);
seq(seq==0) = -1;
%幅值a
A = 10;
seq = A*seq;

% 构造X阵
r = 2;
n = 4;
N = 2^n - 1;
xtemp = repmat(seq, 1, 1+r);
X = zeros(N, r*N);
for i = 1:N
	X(N-i+1,:) = xtemp(1,i:i+r*N-1);
end

%构造输入脉冲,用r+1个周期的m序列输入，得到同样多个响应，再取后面r个周期的
x1 = [0 xtemp];	%rN+1个
y1 = zeros((1+r)*N+1,1);	%列向量

for i=3:(r+1)*N+1
	y1(i) = -c*y1(i-1)-d*y1(i-2) + a*x1(i-1) +b*x1(i-2);
end
Y = y1(end-r*N+1:end);	%r*N x 1

% 矩阵法
Rxy = 1/r/N*X*Y;
g = N/A^2/(N+1)/Ts*(ones(N)+eye(N))* Rxy;

%作图
t1 = [0:N]*Ts;
figure(1);
scatter(t1, [0;g]);
hold on
impulse(Gs,[0:1e-2:1.5])；
legend('m序列','理想脉冲');

• 迭代法

clc
clear

%构造已知传递函数
s = tf('s');
wn = 10;
damp = 0.5;
Gs = wn^2/(s^2 + 2*damp*wn*s + wn^2);
%离散化
Ts = 0.1;
Gz = c2d(Gs,Ts);

% 差分方程系数
a = Gz.Numerator{1}(2);
b = Gz.Numerator{1}(3);
c = Gz.Denominator{1}(2);
d = Gz.Denominator{1}(3);

% 闭环系统，最大工作频率如果取闭环截止频率，2Hz
% 那么不妨取dt=0.1s，N=15。
%构造m序列
coef = [0 0 1 1];
seq = mseq(coef);
seq(seq==0) = -1;
%幅值a
A = 10;
seq = A*seq;
r = 2;
n = 4;
N = 2^n - 1;

%仿真时间
Tsim = 100;
simlen = round(Tsim/Ts);
remainder = mod(simlen-1, N);
nPeriod = round((simlen-1-remainder)/N);

u1 = [0 repmat(seq, 1, nPeriod) seq(1:remainder)];
y1 = zeros(1,simlen);
gm = zeros(N,1);
m = 0;
for i = 3:simlen	% 0 Ts 2Ts ... 
	y1(i) = -c*y1(i-1)-d*y1(i-2) + a*u1(i-1) +b*u1(i-2);	
	if i == 1+2*N	%利用前N个计算Rxy
		m = N-1;
		Rxy = zeros(N,1);
		for j = 0:N-1
			for k = 0:m
				Rxy(j+1) = Rxy(j+1) + y1(i-k)*u1(i-k-j);
			end
		end
		Rxy = Rxy/(m+1);
		gm = N/A^2/(N+1)/Ts*(ones(N) + eye(N))*Rxy;
    end
    if i > 1+2*N
        m=m+1;
        gm = m/(m+1)*gm + N/A^2/(N+1)/Ts/(m+1)*[ones(N)+eye(N)] *y1(i)* u1(i:-1:i-N+1)'; %递推计算
    end
end

tgm = [0:N-1]*Ts;
t1 = [0:simlen-2]*Ts;
figure(1);
scatter(tgm, gm);
hold on
impulse(Gs,[0:1e-2:3])
legend('m序列','理想脉冲');

（四）结果对比

• $t_s=0.1sN=15$
  
• $t_s=0.01N=127$
  

引用

---

1. 佚名.辨识课件——相关分析法[EB/OL].https://wenku.baidu.com/view/3103ae1d0622192e453610661ed9ad51f01d541e.html,2018-05-06. ↩︎ ↩︎ ↩︎ ↩︎ ↩︎

2. Angelo_pj.m序列产生原理及其性质[EB/OL].https://blog.csdn.net/weixin_45015947/article/details/89891757,2019-05-06. ↩︎

3. FastestSnail.基于MATLAB的m序列产生函数及其调用方法[EB/OL].https://blog.csdn.net/cjbct/article/details/78153616,2017-10-03. ↩︎

4. 佚名.系统辨识 第4章 经典辨识方法[EB/OL].https://wenku.baidu.com/view/d5b1e2f47e192279168884868762caaedd33ba89.html,2018-02-05. ↩︎