БПФ, чтобы найти функцию автокорреляции - PullRequest
Новая
       49

БПФ, чтобы найти функцию автокорреляции

1 голос
/ 06 мая 2019

Я пытаюсь найти корреляционную функцию следующего стохастического процесса: enter image description here, где бета и D - постоянные, а xi (t) - гауссовский шумовой термин.

Послемоделируя этот процесс методом Эйлера, я хочу найти функцию автокорреляции этого процесса.Прежде всего, я нашел аналитическое решение для корреляционной функции и уже использовал определение корреляционной функции, чтобы смоделировать ее, и два результата были довольно близки (пожалуйста, смотрите фото, соответствующий код находится в конце этого поста).

figure 1, correlation function analytically and numerically

(рисунок 1)

Теперь я хочу использовать теорему Винера-Хинчина (используя fft), чтобы найтикорреляционная функция, взяв fft из реализаций, умножьте ее на сопряженную и затем найдите возмещение, чтобы получить корреляционную функцию.Но очевидно, что я получаю результаты, далекие от ожидаемой корреляционной функции, поэтому я почти уверен, что в коде есть что-то, что я неправильно понял, чтобы получить неправильные результаты ..

Вот мой код для решениястохастический процесс (который, я уверен, это правильно, хотя мой код может быть небрежным) и моя попытка найти автокорреляцию с помощью 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()

И это то, что я получу, используя этот метод для корреляциифункция, enter image description here

И это мой код для функции корреляции с использованием определения: 

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()

, поэтому кто-то, пожалуйста, укажет, где ошибка вмой код, и как я могу имитировать его лучше, чтобы получить правильную корреляционную функцию?

1 Ответ

1 голос
/ 08 мая 2019

Есть две проблемы с кодом, которые я вижу.

  1. Как сказал Фрэнсис в комментарии , вам нужно вычесть среднее значение из вашего сигнала вполучить автокорреляцию для достижения нуля.

  2. Вы строите свою функцию автокорреляции с неправильными значениями оси X.

    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 раз сделает ее достаточно хорошо выровненной по красной линии.

...