푸리에 변환이란?
푸리에 변환(Fourier transform, FT)은 시간이나 공간에 대한 함수를 시간 또는 공간 주파수 성분으로 분해하는 변환을 말한다.
푸리에 변환은 일종의 적분 변환으로, 리만 이상적분이기 때문에 복잡한 적분 이론을 요구하는 응용분야에서는 적합하지 않을 수 있으며, 푸리에 변환은 파동과 양자역학뿐 아닌 공간이나 운동량 또는 둘 모두를 함수로 표현할 때 파동 공식 표현이 중요한 분야에서 자연스럽게 사용하고 있다.
일반적으로 푸리에 공식이 적용가능한 함수는 복소수이며, 벡터값을 가질 수 있다.
복소수에 푸리에 공식을 적용하는 이유
푸리에 급수/변환의 목적은 어떤 함수 f(t)에 대해 주어진 시간 축에서 다수의 사인파 또는 코사인파로 분해하는 것이다.
현실에서 보통 다루는 함수 f(t)는 보통 실수값 함수이다.
하지만 푸리에 급수/변환 에서는 복소지수함수 e^iωt를 기본 벡터를 사용한다. 이유는 복소지수함수를 이용해 사인과 코사인을 동시에 표현할 수 있기 때문이다.
오일러 공식
푸리에 변환 / 이산 푸리에 변환(DFT)
푸리에 변환
- 입력: 연속시간 함수
- 출력: 연속 주파수 함수
- 아날로그 연속신호
이산 푸리에 변환
- 입력: 유한개의 이산샘플
- 출력: 유한개
- 컴퓨터에서 다루는 디지털데이터(디지털신호)
여기서 \(- \frac{2\pi i}{N} kn\) 는 회전 각도를 나타냄
DFT(이산 푸리에 변환) 코드
import numpy as np
def dft(x):
N = len(x)
X = np.zeros(N, dtype=complex)
for k in range(N):
for n in range(N):
angle = -2j * np.pi * k * n / N
X[k] += x[n] * np.exp(angle)
return X
- 보통 입력값 x는 시간영역의 신호 이고 N의 길이를 가진 복소수 또는 실수 배열이다
- DFT의 output담을 입력값과 같은 길이의 복소수 배열(dtype = complex)을 생성한다
- 이중 for문으로 푸리에 변환 공식을 적용하여 X에 넣은후 return한다.
- 이중 for문으로 구성되어있기 때문에 시간복잡도는 O(N^2)이다.
고속 푸리에 변환 (FFT)
- FFT는 이산 푸리에 변환과 그 역변환을 빠르게 수행하는 알고리즘이다
- DFT는 이중for문으로 구성되었기 때문에 시간복잡도가 O(N^2) 나오는 반면,
- FFT는 분할정복(Divide and Conquer방식을 사용해서 O(NlogN)연산만으로 가능하다.
고속 푸리에 변환 과정 모든 $n \in [0, N - 1]$ 을 짝수/홀수로 나눈다
- 짝수: n = 2m
- 홀수: n = 2m + 1
- 범위: m = 0 ~ N/2 - 1
이와 같이 쪼개진다
\(X[k] = \sum_{m=0}^{N/2 - 1} x[2m] \cdot e^{-2\pi i k (2m)/N} + \sum_{m=0}^{N/2 - 1} x[2m + 1] \cdot e^{-2\pi i k (2m+1)/N}\)
여기서 짝수항은 문제가 없지만 홀수항은 회전 인자 (\(W_N^k = e^{-2\pi i k / N}\)) 를 곱해야 한다
\(O[k] = \sum_{m=0}^{N/2 - 1} x[2m + 1] \cdot e^{-2\pi i k m / (N/2)} \cdot e^{-2\pi i k / N}\)
이처럼 뒤에 회전인사를 곱해줘야 한다
- \(\text{짝수 인덱스: } x[2m] \Rightarrow \text{회전각: } \frac{2\pi i k (2m)}{N} = \frac{2\pi i k m}{N/2} \Rightarrow \text N/2 \text{ 길이의 DFT이다}\)
- \[\text{홀수 인덱스: } x[2m+1] \Rightarrow \text{회전각: } \frac{2\pi i k (2m+1)}{N} = \frac{2\pi i k m}{N/2} + \frac{2\pi i k}{N} \Rightarrow \text{뒤쪽 항 때문에 } W_N^k \text{가 곱해진다}\]
결론적으로 홀수항은 전체기준에서 한칸 밀려있기 때문에 회전인자를 넣어 보정해준다
FFT(고속 푸리에 변환) 코드
import numpy as np
def FFT(x):
x = np.asarray(x, dtype=np.complex128)
N = len(x)
power_2 = 2 ** int(np.ceil(np.log2(N)))
if N != power_2:
x = np.pad(x, (0, power_2 - N))
return FFT_recursive(x)
def FFT_recursive(x):
N = len(x)
if N <= 1:
return x
even = FFT_recursive(x[::2])
odd = FFT_recursive(x[1::2])
W = np.exp(-2j * np.pi * np.arange(N) / N)
return np.concatenate([
even + W[:N // 2] * odd,
even - W[:N // 2] * odd
])
전반적인 코드 설명
- FFT의 핵심은 FFT_recursive에서 even과 odd로 나누는 것이다(분할정복)
- 재귀적으로 반씩 쪼개야하기 때문에 배열의 길이가 2^k이여야 하며, 아닐경우 zere-pending을 수행한다
- zere_pending이란 배열을 조건에 맞춘후 나머지 부분을 0으로 매꾸는 방식이다.
FFT 알고리즘 세부분석
Data, Modules Loading
import numpy as np
FFT 분석
def FFT(x)
def FFT(x):
x = np.asarray(x, dtype=np.complex128)
N = len(x)
power_2 = 2 ** int(np.ceil(np.log2(N)))
if N != power_2:
x = np.pad(x, (0, power_2 - N))
return FFT_recursive(x)
- 입력값 x의 배열을 복소수로 형변환(dtype=np.complex128)한다(FFT는 복소수 연산이기 때문)
- 2^k이고 N >= len(x)인 최소 N을 찾는다 (x의 배열의 길이를 2^k로 맞춰야 하기 때문)
- 배열 x의 길이가 2^k가 아닐경우 빈곳에 0으로 다 채운다 (zero-pading)
- 이후 핵심인 재귀 FFT를 수행
def FFT_recursive(x)
def FFT_recursive(x):
N = len(x)
if N <= 1:
return x
even = FFT_recursive(x[::2])
odd = FFT_recursive(x[1::2])
W = np.exp(-2j * np.pi * np.arange(N) / N)
return np.concatenate([
even + W[:N // 2] * odd,
even - W[:N // 2] * odd
])
세부 분석
N = len(x)
if N <= 1:
return x
- 1일때 까지 계속 반으로 쪼갠다
even = FFT_recursive(x[::2])
odd = FFT_recursive(x[1::2])
- 짝수 부분과 홀수 부분으로 나눈다
W = np.exp(-2j * np.pi * np.arange(N) / N)
- 회전인자 부분이다
return np.concatenate([
even + W[:N // 2] * odd,
even - W[:N // 2] * odd
])
- 마지막으로 짝수와 홀수를 결합해준다
- 이때 홀수항에는 회전인자를 곱해준다
- 복소지수함수의 대칭구조로 인해 +/-를 붙여 반영해햐한다
복소지수함수의 대칭 구조
- \(X[k] = E[k] + W_N^k \cdot O[k]\)
DFT대로면 이런식으로 나와야 한다
- \(X[k + N/2] = E[k] - W_N^k \cdot O[k]\)
하지만 FFT구조에서 상하대칭관계로 인해 이런식으로 나온다.
- \(e^{-2\pi i (k + N/2)/N} = -e^{-2\pi i k/N}\)
- 따라서 위상 보정 후, 앞/뒤로 합칠때 대칭구조를 반영해야 하기때문에 +/-가 붙는다