If you define an LTI system
scipy, you may conveniently feed an input
x to it to get its output
y as follows:
from scipy import signal import numpy as np import matplotlib.pyplot as plt # let's define as mass-damper-spring system: # x'' + beta x' + w0**2 x = F/m m = 10 w0 = 20*2*np.pi beta = 1.5 fs = 500 # sampling frequency (Hz) ts = np.arange(0,100,1./fs) # time interval of simulation sys = signal.lti(1./m,[1,beta, w0**2]) # simulate the response to an input x = 10+np.random.normal(0,1,size=ts.shape) # input _, y, _ = sys.output(U = x, T=ts, ) # output
However, it is also possible convolve the impulse response (which is also accessible via
lti module) with the input and find the output to the input:
_, kern = sys.impulse(T=ts) # impulse response y_conv = signal.convolve(kern/fs,x)[:len(ts)] #1/fs=dt of integration
I realized that the normalized
y_conv and normalized
y differ even though they should be the same theoretically (look below):
fig, axs= plt.subplots(3,1, figsize=(10,5), sharex=True) axs.plot(ts, y,label='output y') axs.scatter(ts, y-y_conv,marker='.', label='difference') axs.plot(ts, kern, label='kernel (impulse response)') plt.xlabel('time') for ax in axs: ax.legend()
Indeed, scipy team didn't implement the output of the linear system using convolution.
I did a similar experiment with an even simpler system (
sys = signal.lti(1,[1,1])) and realized that this error is indeed a function of the sampling frequency
fs, as the following images indicate (note the ranges of middle panels). The higher the
fs the lower the difference.
These pictures suggest that the difference vanishes asymptotically, as also commented below. However, my questions that concern finite sampling rate regimes are:
- How does the sampling frequency cause such differences mathematically? Is there any expression for the error bound?
- What is the role of bilinear mapping in this difference? Is it possible to tailor this mapping to mitigate such numerical discrepancy?
- Which method is more accurate to estimate the response to a general input, and why (implementation-wise)?