Воссоздание данных временного ряда с использованием результатов FFT без использования ifft - PullRequest
Новая
       12

Воссоздание данных временного ряда с использованием результатов FFT без использования ifft

5 голосов
/ 15 декабря 2010

Я проанализировал данные sunspots.dat (ниже), используя fft, который является классическим примером в этой области.Я получил результаты от FFT в реальных и мнимых частях.Затем я попытался использовать эти коэффициенты (первые 20), чтобы воссоздать данные, следуя формуле для преобразования Фурье.Думая, что реальные части соответствуют a_n, а воображаемые - b_n, у меня есть 

import numpy as np
from scipy import *
from matplotlib import pyplot as gplt
from scipy import fftpack

def f(Y,x):
    total = 0
    for i in range(20):
        total += Y.real[i]*np.cos(i*x) + Y.imag[i]*np.sin(i*x)
    return total

tempdata = np.loadtxt("sunspots.dat")

year=tempdata[:,0]
wolfer=tempdata[:,1]

Y=fft(wolfer)
n=len(Y)
print n

xs = linspace(0, 2*pi,1000)
gplt.plot(xs, [f(Y, x) for x in xs], '.')
gplt.show()

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

Данные:

http://linuxgazette.net/115/misc/andreasen/sunspots.dat

1 Ответ

12 голосов
/ 15 декабря 2010

Когда вы вызывали fft(wolfer), вы указали преобразованию принять фундаментальный период, равный длине данных.Чтобы восстановить данные, вы должны использовать базисные функции того же фундаментального периода = 2*pi/N.Точно так же, ваш временной индекс xs должен варьироваться по временным выборкам исходного сигнала.

Другая ошибка заключалась в том, что вы забыли выполнить полное комплексное умножение.Проще думать об этом как Y[omega]*exp(1j*n*omega/N).

Вот фиксированный код.Обратите внимание, что я переименовал i в ctr, чтобы избежать путаницы с sqrt(-1), и с n в N, чтобы следовать обычному соглашению об обработке сигналов, состоящем в использовании нижнего регистра для образца и верхнего регистра для общей длины выборки,Я также импортировал __future__ division, чтобы избежать путаницы с целочисленным делением.

забыл добавить ранее: Обратите внимание, что SciPy fft не делится на N после накопления.Я не делил это перед использованием Y[n];Вы должны это сделать, если хотите вернуть те же числа, а не просто увидеть одну и ту же форму.

И, наконец, обратите внимание, что я суммирую по всему диапазону частотных коэффициентов.Когда я построил график np.abs(Y), казалось, что на верхних частотах были значительные значения, по крайней мере, до выборки 70 или около того.Я подумал, что было бы легче понять результат, суммируя по всему диапазону, увидев правильный результат, затем проанализировав коэффициенты и посмотрев, что произойдет. 

from __future__ import division
import numpy as np
from scipy import *
from matplotlib import pyplot as gplt
from scipy import fftpack

def f(Y,x, N):
    total = 0
    for ctr in range(len(Y)):
        total += Y[ctr] * (np.cos(x*ctr*2*np.pi/N) + 1j*np.sin(x*ctr*2*np.pi/N))
    return real(total)

tempdata = np.loadtxt("sunspots.dat")

year=tempdata[:,0]
wolfer=tempdata[:,1]

Y=fft(wolfer)
N=len(Y)
print(N)

xs = range(N)
gplt.plot(xs, [f(Y, x, N) for x in xs])
gplt.show()
...