ifft를 사용하지 않고 FFT 결과를 사용하여 시계열 데이터 다시 만들기

Question

나는이 영역의 전형적인 예인 fft를 사용하여 sunspots.dat 데이터 (아래)를 분석했습니다. 나는 실제 및 imaginery 부품에서 fft의 결과를 얻었다. 그런 다음이 계수를 사용하여 푸리에 변환 공식을 따라 데이터를 다시 만들려고했습니다. 실제 부품 b_n을 a_n과 imaginery에 해당 생각, 내가 그러나 어떤 이유로

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 (나는 양쪽에 계수의 같은 번호를 사용)에 의해 생성 된 하나를 반영하지 않습니다 내 줄거리. 무엇이 잘못 될 수 있습니까?

은 데이터 : 당신이 fft(wolfer) 전화했을 때

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

Answer 1

, 당신은 데이터의 길이에 해당하는 기본주기를 가정하는 변환했다. 데이터를 재구성하려면 동일한 기본주기 = 2*pi/N의 기본 함수를 사용해야합니다. 같은 토큰으로, 시간 인덱스 xs은 원래 신호의 시간 샘플 범위를 벗어나지 않아야합니다.

또 다른 실수는 완전한 복소 곱셈을하는 것을 잊는 것입니다. 이것을 Y[omega]*exp(1j*n*omega/N)이라고 생각하는 것이 더 쉽습니다.

여기 고정 코드가 있습니다. 참고 sqrt(-1)과 n ~ N의 혼동을 피하기 위해 i을 으로 바꿨습니다. 샘플의 경우 소문자를 사용하고 전체 샘플 길이의 경우 대문자를 사용하는 일반적인 신호 처리 규칙을 따르십시오. 또한 정수 나누기에 대한 혼동을 피하기 위해 __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*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()