이 글은 다음 글을 조금 더 고쳐서 다시 쓴 것이다.

http://snowall.tistory.com/501

파이썬 구현 예제는 생각중이다.

왜 4번이 먼저 올라왔느냐는 중요한 문제가 아니다. 나도 모른다.


# Elementary Numerical analysis 4
# based on Python
# (C) 2013. Keehwan Nam, Dept. of physics, KAIST.
# snowall@gmail.com / snowall@kaist.ac.kr

# Fast Fourier Transform

# 푸리에 변환은 적분변환이다. 다시말해서, 적분을 잘 수행하면 푸리에 변환은 계산된다.

# 0부터 1까지 구간만 생각해 보자.

# F1(k) = Integrate f(x)cos(kx)dx from 0 to 1
# F2(k) = Integrate f(x)sin(kx)dx from 0 to 1

# 위의 두 수를 k를 정해놓고 적분하면 된다. 하지만 현실은 그렇게 만만하지 않은데, 1개의 k에 대해서만 푸리에 계수를 구하는 것이 아니라, k에 대한 함수를 구해야 하는 사태가 벌어지기 때문이다.

# 따라서, 다음과 같이 하자. 조금 바꾼다.

# F = F1+i*F2
# 여기서 i는 제곱하면 -1이 나오는 허수 단위이다.
# 이렇게 되면 F는 오일러 공식을 사용해서 다음과 같이 쓸 수 있다.
# F(k) = Integrate f(x)exp(ikx)dx from 0 to 1
# 이제, 이 함수를 잘 적분하면 되는데, 우리는 지금 수치해석을 공부하고 있으므로 수치적으로 적분하면 된다.

# 0부터 1까지를 N등분하면 dx = 1/N이고 k는 k/N로 생각할 수 있다.
# F(k) = Summation f(x)exp(ikx/N) over integer x from x = 0 to x = N
# 위의 공식을 이산 푸리에 변환Discrete Fourier Transform이라고 한다.

# N등분된 구간은 k가 N개 있으므로, N개의 k값에 대한 각각의 함수를 계산하기 위해서 N개의 값을 더해야 한다. 즉, N*N번의 계산이 필요하다.

# 이제 exp(ikx/N)을 아주 잘 관찰해야 한다.

# f(x)exp(ikx/N)를 계산할 때, 잘 관찰해보면 같은 값들이 반복되는 것을 알 수 있다. k에 대한 값을 계산한 후 k+1에 대한 값을 계산한다면, f(x)exp(i(k+1)x/N)=f(x)exp(ikx/N)exp(ix/N) 이다.
# http://snowall.tistory.com/501
# 이제, 가령 점이 4개만 있다 치고, exp(ikx/N)을 대충 퉁쳐서 w(k*x)라 치고 그냥 대입해보자.

# F(0)=f(0)w(0)+f(1)w(1*0)+f(2)w(1*0)+f(3)w(1*0)
# F(1)=f(0)w(0*1)+f(1)w(1*1)+f(2)w(2*1)+f(3)w(3*1)
# F(2)=f(0)w(0*2)+f(1)w(1*2)+f(2)w(2*2)+f(3)w(3*2)
# F(3)=f(0)w(0*3)+f(1)w(1*3)+f(2)w(2*3)+f(3)w(3*3)

# 뭔가 있어 보인다면, 계산을 완성하자.

# F(0)=f(0)w(0)+f(1)w(0)+f(2)w(0)+f(3)w(0)
# F(1)=f(0)w(0)+f(1)w(1)+f(2)w(2)+f(3)w(3)
# F(2)=f(0)w(0)+f(1)w(2)+f(2)w(0)+f(3)w(2)
# F(3)=f(0)w(0)+f(1)w(3)+f(2)w(2)+f(3)w(1)

# 잘 보면 반복되는 것들이 보인다는 사실을 알 수 있다.
# 따라서 반복되는 부분을 미리 계산하고 재활용하면 된다.

# E1=f(0)w(0)+f(2)w(0)
# E2=f(0)w(0)+f(2)w(2)
# E3=f(1)w(0)+f(3)w(0)
# E4=f(1)w(1)+f(3)w(3)

# F(0)=E1+E3
# F(1)=E2+E4
# F(2)=E1+E3*w(2)
# F(3)=E2+E4*w(2)

# 원래는 16번의 곱셈과 12번의 덧셈이 있었는데, 잘 고쳤더니 10번의 곱셈과 8번의 덧셈으로 줄어들었다.
# 만약 점의 수가 수십만개로 늘어난다면, 아마 그 효율은 매우 좋아질 것이다. 즉 N*N번의 연산이 N*logN번의 연산으로 확 줄어들게 된다.
# 이 알고리즘은 Cooley-Tukey 알고리즘의 특수한 경우이며, 가우스가 최초에 발견하고 쿨리와 튜키가 나중에 다시 발견했다.
# http://ko.wikipedia.org/wiki/고속_푸리에_변환

# 파이썬 구현은 생략한다.

# http://docs.scipy.org/doc/numpy/reference/routines.fft.html 파이썬에서는 numpy 패키지에서 제공한다.


by snowall 2013. 9. 3. 18:33