了解 Scipy 卷积

Question

我试图了解 Scipy 提供的离散卷积与人们将获得的分析结果之间的差异。我的问题是输入信号的时间轴和响应函数如何与离散卷积输出的时间轴相关联？

为了尝试回答这个问题，我考虑了一个带有分析结果的示例。我的输入信号是高斯信号，我的响应函数是具有阶跃函数的指数衰减。这两个信号卷积的解析结果是修改后的高斯（https://en.wikipedia.org/wiki/Exponentially_modified_Gaussian_distribution）。

Scipy 给出了三种卷积模式，“same”、“full”、“valid”。我应用了“相同”和“完整”卷积，并根据下面的分析解决方案绘制了结果。

您可以看到它们非常匹配。

需要注意的一个重要特征是，对于“全”离散卷积，生成的时域比输入信号时域大得多（参见。https://www.researchgate.net/post/How_can_I_get_the_convolution_of_two_signals_in_time_domain_by_just_having_the_values_of_amplitude_and_time_using_Matlab），但对于“相同”离散卷积，时域是相同的，输入和响应域（我在这个例子中选择了相同的域）。

当我观察到更改填充响应函数的域会改变卷积函数的结果时，我产生了困惑。这意味着当我设置 t_response = np.linspace(-5,10,1000) 而不是 t_response = np.linspace(-10,10,1000) 时，我得到了不同的结果，如下所示

如您所见，解决方案略有变化。我的问题是输入信号的时间轴和响应函数如何与输出的时间轴相关联？我在下面附上了我正在使用的代码，我们将不胜感激。

import numpy as np
import matplotlib as mpl
from scipy.special import erf
import matplotlib.pyplot as plt
from scipy.signal import convolve as convolve
params = {'axes.labelsize': 30,'axes.titlesize':30, 'font.size': 30, 'legend.fontsize': 30, 'xtick.labelsize': 30, 'ytick.labelsize': 30}
mpl.rcParams.update(params)

def Gaussian(t,A,mu,sigma):
    return A/np.sqrt(2*np.pi*sigma**2)*np.exp(-(t-mu)**2/(2.*sigma**2))
def Decay(t,tau,t0):
    ''' Decay expoential and step function '''
    return 1./tau*np.exp(-t/tau) * 0.5*(np.sign(t-t0)+1.0)
def ModifiedGaussian(t,A,mu,sigma,tau):
        ''' Modified Gaussian function, meaning the result of convolving a gaussian with an expoential decay that starts at t=0'''
        x = 1./(2.*tau) * np.exp(.5*(sigma/tau)**2) * np.exp(- (t-mu)/tau)
        s = A*x*( 1. + erf(   (t-mu-sigma**2/tau)/np.sqrt(2*sigma**2)   ) )
        return s

### Input signal, response function, analytic solution
A,mu,sigma,tau,t0 = 1,0,2/2.344,2,0  # Choose some parameters for decay and gaussian
t = np.linspace(-10,10,1000)         # Time domain of signal
t_response = np.linspace(-5,10,1000)# Time domain of response function

### Populate input, response, and analyitic results
s = Gaussian(t,A,mu,sigma)
r = Decay(t_response,tau,t0)
m = ModifiedGaussian(t,A,mu,sigma,tau)

### Convolve
m_full = convolve(s,r,mode='full')
m_same = convolve(s,r,mode='same')
# m_valid = convolve(s,r,mode='valid')
### Define time of convolved data
t_full = np.linspace(t[0]+t_response[0],t[-1]+t_response[-1],len(m_full))
t_same = t
# t_valid = t

### Normalize the discrete convolutions
m_full /= np.trapz(m_full,x=t_full)
m_same /= np.trapz(m_same,x=t_same)
# m_valid /= np.trapz(m_valid,x=t_valid)

### Plot the input, repsonse function, and analytic result
f1,(ax1,ax2,ax3) = plt.subplots(nrows=3,ncols=1,num='Analytic')
ax1.plot(t,s,label='Input'),ax1.set_xlabel('Time'),ax1.set_ylabel('Signal'),ax1.legend()
ax2.plot(t_response,r,label='Response'),ax2.set_xlabel('Time'),ax2.set_ylabel('Signal'),ax2.legend()
ax3.plot(t_response,m,label='Output'),ax3.set_xlabel('Time'),ax3.set_ylabel('Signal'),ax3.legend()

### Plot the discrete  convolution agains analytic
f2,ax4 = plt.subplots(nrows=1)
ax4.scatter(t_same[::2],m_same[::2],label='Discrete Convolution (Same)')
ax4.scatter(t_full[::2],m_full[::2],label='Discrete Convolution (Full)',facecolors='none',edgecolors='k')
# ax4.scatter(t_valid[::2],m_valid[::2],label='Discrete Convolution (Valid)',facecolors='none',edgecolors='r')
ax4.plot(t,m,label='Analytic Solution'),ax4.set_xlabel('Time'),ax4.set_ylabel('Signal'),ax4.legend()
plt.show()

Answer 1

问题的症结在于，在第一种情况下，您的信号具有相同的采样率，而在第二种情况下则不然。

我觉得在频域中考虑这一点更容易，其中卷积是乘法。当您创建具有相同时间轴np.linspace(-10, 10, 1000) 的信号和滤波器时，它们现在具有相同的采样率。对每个应用傅里叶变换，得到的阵列为信号和滤波器提供相同频率的功率。所以你可以直接将这些数组的对应元素相乘。

但是，当您创建一个时间轴为 np.linspace(-10, 10, 1000) 的信号和一个时间轴为 np.linspace(-5, 10, 1000) 的滤波器时，情况就不再适用了。应用傅里叶变换并乘以相应元素不再正确，因为相应元素处的频率不一样。

让我们用你的例子来具体说明一下。 信号 (s) 的变换的第一个元素（即np.fft.fftfreq(1000, np.diff(t).mean())[1]）的频率约为0.05 Hz。但是对于 filter (r)，第一个元素的频率大约是0.066 Hz。因此，将这两个元素相乘就是将两个不同频率的功率相乘。 （这就是为什么您经常看到信号处理示例首先定义采样率，然后基于此创建时间数组、信号和滤波器的原因。）

您可以通过创建一个过滤器来验证这是否正确，该过滤器会扩展到您对[-5, 10) 感兴趣的时间范围，但采样率与原始信号相同。所以使用：

t = np.linspace(-10, 10, 1000)
t_response = t[t > -5.0]

在不同的时间范围内以相同的采样率生成信号和过滤器，因此卷积应该是正确的。这也意味着您需要修改每个数组的绘制方式。代码应该是：

ax4.scatter(t_response[::2], m_same[125:-125:2], label='Same') # Conv extends beyond by ((N - M) / 2) == 250 / 2 == 125 on each side
ax4.scatter(t_full[::2], m_full[::2], label='Full')
ax4.scatter(t_response, m, label='Analytic solution')

这会生成下面的图，其中解析、完整和相同的卷积很好地匹配。