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

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

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

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

当我观察到改变我的响应函数的填充域改变了卷积函数的结果时,我产生了混乱.这意味着当我设置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)的滤波器时,那就不再正确了.应用傅立叶变换并乘以相应的元素不再正确,因为相应元素处的频率不相同.

让我们用你的例子来具体化.信号的变换的第一个元素(即,np.fft.fftfreq(1000,np.diff(t).mean())[1])的频率约为0.05Hz.但是对于滤波器(r),第一元件的频率约为0.066Hz.因此,将这两个元素相乘会使两个不同频率的功率相乘. (这种微妙之处就是为什么你经常看到信号处理示例首先定义采样率,然后根据它创建时间数组,信号和滤波器.)

您可以通过创建一个在[-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')

这将生成下图,其中分析,完整和相同的卷积匹配良好.

标签：python,scipy,convolution