我试图找到以下随机过程的相关函数: 其中β和D是常数,席(t)是高斯噪声项。你知道吗
在用欧拉方法模拟这个过程之后,我想找出这个过程的自相关函数。首先,我找到了相关函数的解析解,并且已经用相关函数的定义进行了模拟,两个结果非常接近(请看照片,相应的代码在本文末尾)。你知道吗
(图1)
现在我想用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()
那么,请有人指出我的代码中的错误在哪里,我如何才能更好地模拟它以获得正确的相关函数?你知道吗
我可以看到代码有两个问题。你知道吗
作为francis said in a comment,您需要从信号中减去平均值,以使自相关系数达到零。
用错误的x轴值绘制自相关函数。你知道吗
v_data
定义为:意思是
t
从0到1e4。但是:意思是用
t
从0到1e3绘制。您可能应该使用从情节上看,我认为把蓝线延长10倍会使它和红线很好地吻合。
相关问题 更多 >
编程相关推荐