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

Analicé los datos de sunspots.dat (abajo) usando fft, que es un ejemplo clásico en esta área. Obtuve resultados de FFT en partes reales e imaginarias. Luego intenté 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 e imaginarios 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() 

Sin embargo, por alguna razón, mi gráfica no refleja la generada por ifft (uso el mismo número de coeficientes en ambos lados). Qué podría estar mal ?

Datos:

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

Cuando llamó a fft(wolfer) , le dijo a la transformación que asumiera un período 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 De la misma manera, su índice de tiempo xs tiene que extenderse a lo largo de 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. Nota Cambié el nombre de i a ctr para evitar confusiones con sqrt(-1) , y n a N para seguir la convención usual de procesamiento de señal de usar la minúscula para una muestra y la mayúscula para la longitud total de la muestra. También __future__ division para evitar confusiones acerca de la división entera.

Olvidé añadir antes: tenga en cuenta que los pies de SciPy no se dividen por N después de acumularse. No dividí esto antes de usar Y[n] ; debería hacerlo si desea recuperar los mismos números, en lugar de simplemente ver la misma forma.

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

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