Recrear datos de series temporales usando resultados de FFT sin usar ifft

Question

Analicé los datos sunspots.dat (a continuación) usando fft, que es un ejemplo clásico en esta área. Obtuve resultados de fft en partes reales e imaginarias. Luego traté de usar estos coeficientes (los primeros 20) para recrear los datos siguiendo la fórmula de la transformada de Fourier. Pensando partes reales corresponden a a_n y imaginario a B_N, tengo

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

Por alguna razón, sin embargo, mi argumento no refleja la generada por IFFT (utilizo el mismo número de coeficientes en ambos lados). Qué podría estar mal ?

datos:

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

Answer 1

Cuando llamó fft(wolfer), que le dijo a la transformación para asumir un periodo fundamental igual a la longitud de los datos. Para reconstruir los datos, debe usar funciones básicas del mismo período fundamental = 2*pi/N. Por la misma razón, su índice de tiempo xs tiene que extenderse sobre las muestras de tiempo de la señal original.

Otro error fue olvidarse de hacer la multiplicación compleja completa. Es más fácil pensar en esto como Y[omega]*exp(1j*n*omega/N).

Aquí está el código fijo. Tenga en cuenta que renombré i a ctr para evitar confusiones con sqrt(-1), y n a N para seguir la convención de procesamiento de señal habitual de utilizar la minúscula para una muestra y la mayúscula para la longitud total de la muestra. También importé __future__ division para evitar confusiones sobre la división de enteros.

olvidó añadir antes: Tenga en cuenta que SciPy de fft no divide por N después de haber acumulado. No dividí esto antes de usar Y[n]; deberías hacerlo si quieres recuperar los mismos números, en lugar de solo ver la misma forma.

Y, por último, tenga en cuenta que estoy sumando la gama completa de coeficientes de frecuencia. Cuando tracé np.abs(Y), parecía que había valores significativos en las frecuencias superiores, al menos hasta la muestra 70 más o menos. Pensé que sería más fácil entender el resultado sumando todo el rango, viendo el resultado correcto, luego recortando los coeficientes y viendo qué pasa.

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*pi/N) + 1j*np.sin(x*ctr*2*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()