理解Scipy卷积

Question

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

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

Scipy给出了三种卷积模式,"相同","完整","有效".我应用了"相同"和"完整"卷积,并根据下面的分析解决方案绘制了结果.

你可以看到它们非常匹配.

需要注意的一个重要特性是,对于"完整"离散卷积,得到的时域远大于输入信号时域(参见.http://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)，这不再是真的。应用傅里叶变换并乘以相应元素不再正确，因为相应元素处的频率不相同。

让我们用你的例子把它具体化。信号( )的变换（即np.fft.fftfreq(1000, np.diff(t).mean())[1]）的第一个元素的频率约为。但对于滤波器( )，第一个元素的频率约为。因此，将这两个元素相乘就是将两个不同频率的功率相乘。（这种微妙之处就是为什么您经常看到信号处理示例首先定义采样率，然后基于此创建时间数组、信号和滤波器。）s0.05 Hzr0.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')

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