如何用SciPy Signal无循环获取SISO系统对多输入信号的响应

嘿，刚好碰到过类似的需求！你现在的问题是scipy.signal.lsim默认只能处理单个输入序列，而你想直接把多个输入信号（就是你的x数组）传进去，不用一个个写x[0]、x[1]...对吧？下面给你两种可行的解决办法：

方法一：用numpy.apply_along_axis批量处理 #

这个方法最直观，相当于让numpy帮你自动遍历每一行输入信号，调用lsim计算响应，不用你手动写循环：

from __future__ import division
import numpy as np
from scipy import signal

nbr_inputs = 5
t_in = np.arange(0, 10, 0.2)
dim = (nbr_inputs, len(t_in))
x = np.cumsum(np.random.normal(0, 2e-3, dim), axis=1)

H = signal.TransferFunction([1, 3, 3], [1, 2, 1])

# 先写个小函数，处理单个输入序列的响应计算
def get_response(input_seq, time_vec, sys):
    _, output, _ = signal.lsim(sys, input_seq, time_vec)
    return output

# 对x的每一行（每个输入信号）应用这个函数
all_responses = np.apply_along_axis(get_response, axis=1, arr=x, time_vec=t_in, sys=H)

# 结果all_responses是(nbr_inputs, len(t_in))的数组，每一行对应一个输入的响应

方法二：利用LTI系统的线性叠加性（卷积法） #

因为你的系统是线性时不变的，我们可以先算出系统的脉冲响应，然后用卷积来计算每个输入的响应——这种方法有时候效率更高，尤其是输入信号很多的时候：

# 先获取系统的脉冲响应（和输入时间轴对齐）
_, impulse_response = signal.impulse(H, T=t_in)
dt = t_in[1] - t_in[0]  # 时间步长

# 对每个输入信号做卷积，乘以dt是为了近似连续卷积的积分
all_responses = np.array([np.convolve(input_seq, impulse_response, mode='same') * dt for input_seq in x])

这里要注意mode='same'能保证输出和输入的长度一致，刚好匹配你的时间轴。

为什么直接传x给lsim不行？ #

顺便说下，你直接传x给lsim会出错或者得到错误结果，是因为lsim把二维的x当成多输入系统的多个输入通道，而不是多个独立的SISO输入序列——这和你的需求完全不是一回事儿，所以必须用上面的方法来批量处理。

内容的提问来源于stack exchange，提问作者s_o_king