Как сделать БПФ очень длинного сигнала, используя Numpy, чтобы найти PSD и корреляционную функцию - PullRequest
0 голосов
/ 17 февраля 2020

Я хочу найти PSD и корреляционную функцию стационарного процесса с очень долгим временем, потому что я ищу чрезвычайно точные результаты (dt = 10 ^ -3 и N = 10 ^ 9). Я хочу использовать метод fft и ifft, чтобы найти его, поскольку он самый быстрый в численном отношении. (теорема Винера-Хинчина)

Проблема в том, что fft, использующий numpy или scipy, не работает для таких больших значений N, поэтому мне интересно, возможно ли, например, нарезать сигнал, находя PSD срезов перед добавлением их снова и предварительно формируя ifft, чтобы найти корреляционную функцию. Или есть ли альтернативы библиотекам scipy или numpy fft, которые могут работать с такими длинными массивами, или есть какие-то хитрости для go в отношении этой проблемы.

Процесс, который я ищу, является стационарным процессом с экспоненциальным стационарным распределением, и это мой код процесса и функции корреляции:

def process(tau, var, N, dt): ## var and tau are parameters of the process
    v = np.random.laplace(0, np.sqrt(var / 2))
    gamma = (5 / ( 2 * np.sqrt(2))) * np.sqrt(var) / tau
    D = (5/4) * var / tau
    diff = np.sqrt(2 * D * dt)
    drift = gamma * dt
    v_data = np.zeros(N)
    v_data[0] = v
    for j in range(1, N):
        F = np.random.normal(0,1)
        if v >= 0:
            v = v - drift + diff * F
        elif v < 0:
            v = v + drift + diff*  F
        v_data[j] = v
    return v_data
def Corr_sim(tau, var):
    N = 10 ** 9
    dt =  10 ** - 3
    t = np.arange(N) * dt
    v_data = process(tau, var, N, dt)
## fft for such a large N always gives memory error
    psd = np.fft.fft(v_data) * np.conj(np.fft.fft(v_data)) / N 
    C = np.fft.ifft(psd).real - np.mean(v_data)**2
    return C
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...