了解 Scipy 卷积
Understanding Scipy Convolution
我想了解 Scipy 提供的离散卷积与获得的分析结果之间的差异。我的问题是输入信号的时间轴和响应函数如何与离散卷积输出的时间轴相关联？
为了尝试回答这个问题，我考虑了一个具有分析结果的示例。我的输入信号是高斯信号，我的响应函数是阶跃函数的指数衰减。这两个信号卷积的解析结果是修正的高斯分布（https://en.wikipedia.org/wiki/Exponentially_modified_Gaussian_distribution).
Scipy给出三种卷积模式，"same"、"full"、"valid"。我应用了 "same" 和 "full" 卷积，并根据下面的分析解决方案绘制了结果。
你可以看到他们非常匹配。
需要注意的一个重要特征是，对于"full"离散卷积，生成的时域比输入信号时域大得多（参见。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），但对于"same" 离散卷积时域相同，输入域和响应域相同（我在这个例子中选择相同）。
当我观察到改变填充我的响应函数的域会改变卷积函数的结果时，我感到很困惑。这意味着当我设置 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()
问题的症结在于，在第一种情况下，您的信号具有相同的采样率，而在第二种情况下，它们却不同。
我觉得在频域考虑这个更容易，卷积就是乘法。当您创建具有相同时间轴 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')
这会生成下面的图，其中解析的、完整的和相同的卷积匹配得很好。
我想了解 Scipy 提供的离散卷积与获得的分析结果之间的差异。我的问题是输入信号的时间轴和响应函数如何与离散卷积输出的时间轴相关联？
为了尝试回答这个问题，我考虑了一个具有分析结果的示例。我的输入信号是高斯信号，我的响应函数是阶跃函数的指数衰减。这两个信号卷积的解析结果是修正的高斯分布（https://en.wikipedia.org/wiki/Exponentially_modified_Gaussian_distribution).
Scipy给出三种卷积模式，"same"、"full"、"valid"。我应用了 "same" 和 "full" 卷积，并根据下面的分析解决方案绘制了结果。
你可以看到他们非常匹配。
需要注意的一个重要特征是，对于"full"离散卷积，生成的时域比输入信号时域大得多（参见。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），但对于"same" 离散卷积时域相同，输入域和响应域相同（我在这个例子中选择相同）。
当我观察到改变填充我的响应函数的域会改变卷积函数的结果时，我感到很困惑。这意味着当我设置 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()
问题的症结在于，在第一种情况下，您的信号具有相同的采样率，而在第二种情况下，它们却不同。
我觉得在频域考虑这个更容易，卷积就是乘法。当您创建具有相同时间轴 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')
这会生成下面的图，其中解析的、完整的和相同的卷积匹配得很好。