FFT 求自相关函数

Question

我试图找到以下随机过程的相关函数：其中 beta 和 D 是常数，xi(t) 是高斯噪声项。

用欧拉法模拟这个过程后，想求出这个过程的自相关函数。首先，我已经找到了相关函数的解析解，并且已经使用相关函数的定义来模拟它，两个结果非常接近（请看照片，相应的代码在这个末尾post).

（图一）

现在我想使用 Wiener-Khinchin 定理（使用 fft）通过取实现的 fft 求相关函数，将它乘以它的共轭然后求 ifft 得到相关函数.但显然我得到的结果与预期的相关函数相去甚远，所以我很确定我在代码中误解了一些东西来得到这个错误的结果。

这是我解决随机过程的代码（我确信它是正确的，尽管我的代码可能很草率）和我尝试用 fft 找到自相关：

N = 1000000
dt=0.01
gamma = 1
D=1
v_data = []
v_factor = math.sqrt(2*D*dt)
v=1
for t in range(N):
        F = random.gauss(0,1)
        v = v - gamma*dt + v_factor*F
        if v<0: ###boundary conditions.
            v=-v
        v_data.append(v)


def S(x,dt):  ### power spectrum 
    N=len(x)
    fft=np.fft.fft(x)
    s=fft*np.conjugate(fft)
 #   n=N*np.ones(N)-np.arange(0,N) #divide res(m) by (N-m)
    return s.real/(N)

c=np.fft.ifft(S(v_data,0.01))  ### correlation function 
t=np.linspace(0,1000,len(c))

plt.plot(t,c.real,label='fft method')
plt.xlim(0,20)
plt.legend()
plt.show()

这就是我使用这种方法获得的相关函数，

这是我使用定义的相关函数的代码：

def c_theo(t,b,d): ##this was obtained by integrating the solution of the SDE 
    I1=((-t*d)+((d**2)/(b**2))-((1/4)*(b**2)*(t**2)))*special.erfc(b*t/(2*np.sqrt(d*t)))
    I2=(((d/b)*(np.sqrt(d*t/np.pi)))+((1/2)*(b*t)*(np.sqrt(d*t/np.pi))))*np.exp(-((b**2)*t)/(4*d))
    return I1+I2 

## this is the correlation function that was plotted in the figure 1 using the definition of the autocorrelation. 
Ntau = 500
sum2=np.zeros(Ntau)
c=np.zeros(Ntau)
v_mean=0

for i in range (0,N):
    v_mean=v_mean+v_data[i]
v_mean=v_mean/N
for itau in range (0,Ntau):
    for i in range (0,N-10*itau):
            sum2[itau]=sum2[itau]+v_data[i]*v_data[itau*10+i]
    sum2[itau]=sum2[itau]/(N-itau*10)    
    c[itau]=sum2[itau]-v_mean**2

t=np.arange(Ntau)*dt*10
plt.plot(t,c,label='numericaly')
plt.plot(t,c_theo(t,1,1),label='analyticaly')
plt.legend()
plt.show()

请有人指出我的代码中的错误在哪里，我怎样才能更好地模拟它以获得正确的相关函数？

Answer 1

我看到的代码有两个问题。

与一样，您需要从信号中减去均值以使自相关达到零。
您使用错误的 x-axis 值绘制自相关函数。

v_data 定义为：
```
 N = 1000000   % 1e6
 dt = 0.01     % 1e-2
```
意味着 t 从 0 到 1e4。然而：
```
 t = np.linspace(0,1000,len(c))
```
意味着你用 t 从 0 到 1e3 绘图。您可能应该用
定义 t
```
 t = np.arange(N) * dt
```
看图，我想说将蓝线拉伸 10 倍会使它与红线很好地对齐。

FFT 求自相关函数

FFT to find autocorrelation function

python

signal-processing

numpy

fft

correlation