脉冲响应函数python代码 脉冲响应函数结果分析

求解系统脉冲响应

• 前言
• 相关分析与系统

• （一）脉冲输入
• （二）白噪声输入
• （三）m序列输入
• （四）矩阵法
• （五）迭代法
• （六）参数选择

• 举例说明

• （一）问题
• （二）参数选择
• （三）代码
• （四）结果对比

• 引用

前言

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

相关分析与系统

（一）脉冲输入

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

输入输出关系为：
$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^2g(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{cases} a^2,\quad &t = 0 \\ -\frac{a^2}{N}, &t = t_s,2t_s,...(N-1)t_s \end{cases}$

设$\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) \triangle$

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

写成矩阵形式：$R_{xy} = \triangle Rg$
其中：
$R_{xy}= \begin{bmatrix} R_{xy}(0) \\ R_{xy}(1) \\ \vdots \\ R_{xy}(N-1) \end{bmatrix} g = \begin{bmatrix} g(0) \\ g(1) \\ \vdots \\ g(N-1) \end{bmatrix}$

$\begin{aligned} R &= \begin{bmatrix} R_{xx}(0) &R_{xx}(1) &\dots &R_{xx}(N-1) \\ R_{xx}(-1) &R_{xx}(0) &\dots &R_{xx}(N-2) \\ \vdots &\vdots &\ddots &\vdots \\ R_{xx}(-N+1) &R_{xx}(-N+2) &\dots &R_{xx}(0) \end{bmatrix} \\ &= a^2 \begin{bmatrix} 1 &-1/N &\dots &-1/N \\ -1/N &1 &\dots &-1/N \\ \vdots &\vdots &\ddots &\vdots \\ -1/N &-1/N &\dots &1 \end{bmatrix} \end{aligned}$

这时得到$g = \frac{1}{\triangle} R^{-1}R_{xy}$。
而根据互相关函数$R_{xy}$的定义可求：
$R_{xy}(\tau) = \frac{1}{N}\sum_{k=0}^{N-1}x(k-\tau)y(k)$
写成矩阵形式：
$R_{xy} = \begin{bmatrix} x(0) &x(1) &\dots &x(N-1) \\ x(-1) &x(0) &\dots &x(N-2) \\ \vdots &\vdots &\ddots &\vdots \\ x(-N+1) &x(-N+2) &\dots &x(0) \end{bmatrix} \begin{bmatrix} y(0) \\ y(1) \\ \vdots \\ y(N-1) \end{bmatrix} =XY$

所以最终得到¹：
$g = \frac{1}{a^2(N+1)t_s} \begin{bmatrix} 2 &1 &\dots &1 \\ 1 &2 &\dots &1 \\ \vdots &\vdots &\ddots &\vdots \\ 1 &1 &\dots &2 \end{bmatrix} XY$

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

（五）迭代法

迭代法在引用一¹

（六）参数选择

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

举例说明

（一）问题

尝试绘制典型二阶系统的脉冲响应的波形：
$G(s) = \frac{w_n^2}{s^2+2\zeta w_ns+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$进行试验。
• 根据$Nt_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.1s \quad N=15$
• $t_s = 0.01 \quad N=127$

引用

---

1. 佚名.辨识课件——相关分析法[EB/OL].https://wenku.baidu.com/view/3103ae1d0622192e453610661ed9ad51f01d541e.html,2018-05-06. ↩︎ ↩︎ ↩︎ ↩︎ ↩︎
2. Angelo_pj.m序列产生原理及其性质[EB/OL].. ↩︎
3. FastestSnail.基于MATLAB的m序列产生函数及其调用方法[EB/OL].. ↩︎
4. 佚名.系统辨识 第4章 经典辨识方法[EB/OL].https://wenku.baidu.com/view/d5b1e2f47e192279168884868762caaedd33ba89.html,2018-02-05. ↩︎