고속 푸리에 변환

FFT는 여기로 연결됩니다. 다른 뜻에 대해서는 FFT (동음이의) 문서를 참고하십시오.

고속 푸리에 변환(高速 푸리에 變換, 영어: Fast Fourier Transform, FFT)은 이산 푸리에 변환(영어: Discrete Fourier Transform, DFT)과 그 역변환을 빠르게 수행하는 효율적인 알고리즘이다. FFT는 디지털 신호 처리에서 편미분 방정식의 근을 구하는 알고리즘에 이르기까지 많은 분야에서 사용한다.

신호 $x_{0},...,x_{N-1}$이 복소수라고 가정할 때, DFT는 다음과 같이 정의한다.

$f_{k}=\sum _{n=0}^{N-1}x_{n}e^{-{2\pi i \over N}kn}\qquad \qquad k=0,\dots ,N-1\qquad \qquad \qquad (1)$

이 식을 정의에 따라 계산하면 O(N ²)의 연산이 필요하지만, FFT를 이용하면 O(N log N)의 연산만으로 가능하다. N이 합성수일 경우 그 소인수 분해를 이용하여 연산 횟수를 줄일 수는 있지만, FFT를 사용하면 N이 소수일 경우에도 O(N log N)번의 연산 횟수를 보장한다.

$W_{N}=e^{-2\pi i/N}$ 을 회전자(Twiddle factor)라고 부른다. 1의 거듭제곱근(Root of unity), 단위근도 회전자의 다른 호칭이다. 이 회전자를 사용해 식(1)을 다시 정리하면,

$f_{k}=\sum _{n=0}^{N-1}x_{n}W_{N}^{kn}\qquad \qquad k=0,\dots ,N-1\qquad \qquad \qquad (2)$

회전자는 N에 관한 주기성으로 인해 다음과 같은 성질이 있음을 쉽게 알수있다.

$W_{N}^{k+N}=W_{N}^{k}\qquad \because (W_{N}^{k+N}=e^{-{2\pi i \over N}(k+N)}=e^{-{2\pi i \over N}k}e^{-2\pi i}=e^{-{2\pi i \over N}k}(1)=W_{N}^{k})\qquad (3a)$

$W_{N}^{k+N/2}=-W_{N}^{k}\quad \because (W_{N}^{k+N/2}=e^{-{2\pi i \over N}(k+{N \over 2})}=e^{-{2\pi i \over N}k}e^{-\pi i}=e^{-{2\pi i \over N}k}(-1)=-W_{N}^{k})\qquad (3b)$

$W_{a}^{b}=W_{na}^{nb}\qquad$n:자연수$\qquad \qquad (3c)$

(예 1)

$W_{4}^{6}=W_{4}^{2},\qquad W_{2}^{1}=W_{4}^{2}=W_{8}^{4},\qquad W_{4}^{(k+4/2)}=W_{4}^{k}W_{4}^{2}=-W_{4}^{k}$

쿨리-튜키 알고리즘

분할(Decimation)

어떤 길이 N인 수열을 다음과 같이 index가 짝수인 것과 홀수인 것들만 각각 모아서 두 개의 수열로 나누는 것이다.

$(x_{0},x_{1},x_{2},...,x_{N-2},x_{N-1})$→$(x_{0},x_{2},x_{4},...,x_{N-4},x_{N-2})+(x_{1},x_{3},x_{5},...,x_{N-3},x_{N-1})$

분할 대신 일정 부분을 뽑아낸다는 '솎음'으로 번역한 경우도 있는데 단어 번역으로는 더 정확하지만, FFT에서는 오히려 분할이라는 표현이 이해하는데 더 낫다.

다니엘슨-란초스 보조정리(Danielson-Lanczos Lemma)

1942년 위 두사람은 FFT의 기초가 되는 공식을 발견한다. 즉, 길이 N(단 짝수인 경우)의 DFT는 각각 길이가 N/2인 두 개의 DFT의 합으로 다시 쓸 수 있다는 것. 하나는 짝수 index 포인트들로 형성되고, 다른 것은 홀수 index 포인트들로 만들어진다. 이 분할(decimation)을 수식으로 표현해 보자.

식(2)의 DFT를 짝수번째 항, 홀수 번째 항끼리 묶으면 즉, 분할하면

$f_{k}={\Bigl [}x_{0}W_{N}^{0}+x_{2}W_{N}^{2k}+x_{4}W_{N}^{4k}+...+x_{(N-4)}W_{N}^{k(N-4)}+x_{(N-2)}W_{N}^{k(N-2)}{\Bigr ]}$

$+{\Bigl [}x_{1}W_{N}^{k}+x_{3}W_{N}^{3k}+x_{5}W_{N}^{5k}+...+x_{(N-3)}W_{N}^{k(N-3)}+x_{(N-1)}W_{N}^{k(N-1)}{\Bigr ]}$

길이가 짝수이므로 각 항의 갯수는 같다. 두번째 괄호 항에서 $W_{N}^{k}$을 밖으로 빼면

$f_{k}={\Bigl [}x_{0}W_{N}^{0}+x_{2}W_{N}^{2k}+x_{4}W_{N}^{4k}+...+x_{(N-4)}W_{N}^{k(N-4)}+x_{(N-2)}W_{N}^{k(N-2)}{\Bigr ]}$

$+W_{N}^{k}{\Bigl [}x_{1}+x_{3}W_{N}^{2k}+x_{5}W_{N}^{4k}+...+x_{(N-3)}W_{N}^{k(N-4)}+x_{(N-1)}W_{N}^{k(N-2)}{\Bigr ]}$

이를 다시 쓰면

$f_{k}=\sum _{n=0}^{N/2-1}x_{2n}W_{N}^{2kn}+W_{N}^{k}\sum _{n=0}^{N/2-1}x_{2n+1}W_{N}^{2kn}$

그런데, $W_{N}^{2}=e^{-4\pi i/N}=e^{-2\pi i/(N/2)}=W_{N/2}$ 가 성립하므로,

$f_{k}=\sum _{n=0}^{N/2-1}x_{2n}W_{N/2}^{kn}+W_{N}^{k}\sum _{n=0}^{N/2-1}x_{2n+1}W_{N/2}^{kn}$

$=E_{k}+W_{N}^{k}O_{k}\qquad \qquad k=0,\dots ,N/2-1\qquad \qquad \qquad$

여기서, $E_{k},O_{k}$ 는 위 수식에서 보였듯이 각각 짝수(Even)항과 홀수(Odd)항으로 구성된 N/2-DFT이다.

다음으로 f_k+N/2를 구해보자. 윗식으로부터

$f_{k+N/2}=\sum _{n=0}^{N/2-1}x_{2n}W_{N/2}^{(k+N/2)n}+W_{N}^{k+N/2}\sum _{n=0}^{N/2-1}x_{2n+1}W_{N/2}^{(k+N/2)n}$

$=\sum _{n=0}^{N/2-1}x_{2n}W_{N/2}^{kn}-W_{N}^{k}\sum _{n=0}^{N/2-1}x_{2n+1}W_{N/2}^{kn}\qquad \qquad$ from (3a),(3b)

$=E_{k}-W_{N}^{k}O_{k}\qquad \qquad k=0,\dots ,N/2-1\qquad \qquad \qquad$

f_k+N/2 의 실제 범위는 f_k 와 달리 k에 N/2가 추가되었으므로 $k=N/2,\dots ,N-1$ 임에 주의하자. 정리하면,

$f_{k}=E_{k}+W_{N}^{k}O_{k}\qquad \qquad k=0,\dots ,N/2-1\qquad \qquad \qquad$

$f_{k+N/2}=E_{k}-W_{N}^{k}O_{k}\qquad \qquad k=0,\dots ,N/2-1\qquad \qquad (4)$

식(4)의 중요성은 다음 쿨리-튜키 알고리즘에서 설명하겠다. 이 식을 그림(1)로 나타낼 수 있으며, 마치 나비 모양이라고 해서 이런 분할 과정을 나비연산(butterfly operation)이라 부르고, 그 그림을 나비(butterfly)라고 한다.

그림 1. 기본적인 기수 2 시분할 나비

역비트순(Bit-reversed Order)

만약 N이 2의 거듭제곱이면 즉 N=2ⁿ 이면 분할(decimation)에의해 길이가 N/2 인 두신호는 길이가 N/4 인 두신호로 쪼개어진다. 이런 방식으로 계속 분할을 진행하면 결국 길이가 2인(N=2) 신호로 나누어진다. N=2 인 DFT는 매우 쉽게 계산된다. 그림(2)와 같이 길이가 8인(N=8)인 신호를 생각하자. 처음 길이가 8인 신호는 다니엘슨-란초스 보조정리에의해 분할을 계속하면 stage 3과 같은 순서가 된다.

그림 2. N=8 일 때 분할에 기인한 역비트순

그런데 위의 표에서 처음 원래의 신호와 최종적으로 나누어진 신호의 순서를 이진수로 나타내면 그림에서와 같이 대응하는 신호의 모든 비트들을 거꾸로 배열한 것과 같다. 이것을 Bit reversal 이라 부른다. 이런 현상은 신호가 N=2ⁿ 라면 항상 나타나는 현상이다. 계산구조가 in-place(입력 신호 메모리 장소에 출력 값을 덮어쓰는 것. 다른 표현으로 원래의 샘플 메모리 안에서 FFT를 수행하고 추가적인 버퍼 메모리를 사용하지않는 것)이어야 한다는 요구조건을 포기한다면 입력자료와 출력 결과 둘 다 순서대로 나열되도록 하는 것이 가능하며, 이를 out-of-place라 한다.

쿨리-튜키 알고리즘(Cooley-Tukey algorithm)

가장 일반적으로 사용되는 FFT 알고리즘은 쿨리-튜키 알고리즘이다. 이 방법은 1965년에 J. W. 쿨리와 J. W. 튜키가 발표하여 널리 알려졌지만, 나중에 카를 프리드리히 가우스가 1805년에 같은 방법을 이미 개발하였으며, 그 뒤에도 제한된 형태의 FFT가 종종 발견되었음이 밝혀졌다.

식(2)의 경우 k =0 에서 N-1까지의 각 항의 계산에 N 회의 곱셈이 필요하기 때문에, 전체적으로 N ² 회의 곱셈이 필요하게 된다. 그러나, 식(4)의 길이 N/2 의 DFT는 솔직하게 계산하면 (N/2) ²=N ²/4 회의 곱셈으로 실행할 수 있으므로 이 분할(Decimation)에 의해 계산량은 약 절반으로 줄어든다. 또한, 이 분할을 2회, 3회,... 반복하면 즉, 재사용(reuse)하면 계산량은 약 1/4, 1/8,... 로 격감한다. 이것이 쿨리-튜키 FFT (정확하게는 기수 2 쿨리-튜키 FFT)의 기본적인 생각이다.

쿨리-튜키 알고리즘은 보통 크기 N을 재귀적으로 2등분하여 길이가 2인 신호가 얻어질 때까지 분할 정복 알고리즘(Divide and conquer algorithm)을 적용하기 때문에 N = 2ⁿ인 경우에 많이 적용된다.

보통의 프로그램 코드는 log N 개의 단계(그림(2)나 그림(3)의 stage에 해당)에 대한 재귀 호출을 수행하고, 각 단계에서 N / 2 회의 복소수 곱셈과 N 회의 복소수 가산을 수행한다. 따라서 복소수 곱셈 횟수는 (N / 2) * log N 로 감소하고 부동 소수점 연산의 양은 N * log N의 order가 될 수 있다. 이것은 쿨리-튜키 FFT의 일반적인 연산량이다. 일반적으로 N = n₁ n₂가 성립하는 n₁과 n₂는 같을 필요가 없으며, 따라서 N이 임의의 합성수일 때에도 적용 가능하다(아래 혼합 기수 FFT 를 참고). 쿨리-튜키 알고리즘은 DFT를 더 작은 크기의 DFT로 나누기 때문에, 뒤에 설명된 다른 FFT 알고리즘과 함께 사용될 수 있다.

푸리에 행렬(Fourier Matrix)

회전자로 이루어진 차수(order) N의 복소 방데르몽드 행렬((Vandermonde matrix)

$F_{jk}=e^{-2\pi i(j-1)(k-1)/N}$

을 푸리에 행렬(Fourier Matrix)이라고 한다. 식(2)를 아래와 같이 푸리에 행렬로 표현할 수 있다. 편의를 위해 지금부터는 회전자의 아래첨자 N을 생략하겠다. 예를들면, $W^{2}=W_{N}^{2}$이다.

${\begin{bmatrix}f_{0}\\f_{1}\\f_{2}\\f_{3}\\.\\.\\f_{N-1}\end{bmatrix}}={\begin{bmatrix}W^{0}&W^{0}&W^{0}&.&.&.&W^{0}\\W^{0}&W^{1}&W^{2}&.&.&.&W^{N-1}\\W^{0}&W^{2}&W^{4}&.&.&.&W^{2(N-1)}\\W^{0}&W^{3}&W^{6}&.&.&.&W^{3(N-1)}\\.&.&.&.&&&.&\\.&.&.&&.&&.&\\W^{0}&W^{N-1}&W^{2(N-1)}&.&.&.&W^{(N-1)^{2}}\\\end{bmatrix}}{\begin{bmatrix}x_{0}\\x_{1}\\x_{2}\\x_{3}\\.\\.\\x_{N-1}\end{bmatrix}}\qquad \qquad \qquad (5)$

식(3a)가 보여주는 회전자의 주기성으로 인해 $W^{0}$부터 $W^{(N-1)^{2}}$까지 모두 계산할 필요가 없이 $W^{(N-1)}$까지만 계산하면 된다. 다음 기수 2 시분할 FFT의 식(6)에서 확인해 볼 수 있다.

기수 2 알고리즘(Radix-2 FFT Algorithm)

앞에 열거한 내용들을 기초하여 기수 2 FFT 알고리즘을 다루어본다. 이는 다른 알고리즘들을 다루는데 있어 가장 기초가 되는 알고리즘이다. 우선 기수의 뜻을 살펴보자.

기수(radix)의 의미

- 기수(基數)는 밑 또는 기저(radix, base)로도 부르고, 숫자 표현(진법 체계)의 기준이 되는 수다. radix 는 라틴어로 뿌리(root)라는 뜻이다.

(예) 10진수의 기수(base)는 10이다. 10 ¹, 10 ² , 10 ³ 의 밑(base)은 10이다.

- 거듭제곱 r ⁿ의 기수(radix)는 r이다. 여기서, n는 지수(exponent)이다.

(예) N=2 ⁿ는 기수 2(radix-2)이다. N=4 ⁿ는 기수 4(radix-4)이다.

시분할(DIT) FFT와 주파수분할(DIF) FFT

FFT를 수행하는 방식에는 두가지 알고리즘이 있다.

1) 입력(일반적으로 시간 domain)을 나누는 방식 - 시분할(DIT, Decimation in Time) FFT

2) 출력(일반적으로 주파수 domain)을 나누는 방식 - 주파수분할(DIF, Decimation in Frequency) FFT

다음은 이들에 대해 설명한다.

1. 기수 2 시분할 FFT(Radix-2 Decimation in Time FFT)

N = 4 일 때, 기수 2 시분할 FFT를 유도해 보자. 우선 식(5)로부터 N = 4 DFT 행렬식은 다음과 같다.

${\begin{bmatrix}f_{0}\\f_{1}\\f_{2}\\f_{3}\end{bmatrix}}={\begin{bmatrix}W^{0}&W^{0}&W^{0}&W^{0}\\W^{0}&W^{1}&W^{2}&W^{3}\\W^{0}&W^{2}&W^{4}&W^{6}\\W^{0}&W^{3}&W^{6}&W^{9}\end{bmatrix}}{\begin{bmatrix}x_{0}\\x_{1}\\x_{2}\\x_{3}\end{bmatrix}}\qquad \qquad \qquad (5a)$

지수의 짝 홀을 기준으로 위의 식을 다음과 같이 열의 순서를 바꾸는 변형할 수 있다(이는 곧, 앞서 설명한 분할(decimation)을 의미한다). 또한 식(3a)로부터 다음의 관계를 얻으므로

$W^{4}=W^{0},\quad W^{6}=W^{2+4}=W^{2+0}=W^{2}$${\begin{bmatrix}f_{0}\\f_{1}\\f_{2}\\f_{3}\end{bmatrix}}={\begin{bmatrix}W^{0}&W^{0}&W^{0}&W^{0}\\W^{0}&W^{2}&W^{1}&W^{3}\\W^{0}&W^{4}&W^{2}&W^{6}\\W^{0}&W^{6}&W^{3}&W^{9}\end{bmatrix}}{\begin{bmatrix}x_{0}\\x_{2}\\x_{1}\\x_{3}\end{bmatrix}}={\begin{bmatrix}W^{0}&W^{0}&W^{0}W^{0}&W^{0}W^{0}\\W^{0}&W^{2}&W^{1}W^{0}&W^{1}W^{2}\\W^{0}&W^{0}&W^{2}W^{0}&W^{2}W^{0}\\W^{0}&W^{2}&W^{3}W^{0}&W^{3}W^{2}\end{bmatrix}}{\begin{bmatrix}x_{0}\\x_{2}\\x_{1}\\x_{3}\end{bmatrix}}\qquad (6)$

식(6)은 다음과 같이 다시 쓸 수 있다.

${\begin{bmatrix}f_{0}\\f_{1}\\f_{2}\\f_{3}\end{bmatrix}}={\begin{bmatrix}1&0&W^{0}&0\\0&1&0&W^{1}\\1&0&W^{2}&0\\0&1&0&W^{3}\end{bmatrix}}\,{\begin{bmatrix}W^{0}&W^{0}&0&0\\W^{0}&W^{2}&0&0\\0&0&W^{0}&W^{0}\\0&0&W^{0}&W^{2}\end{bmatrix}}{\begin{bmatrix}x_{0}\\x_{2}\\x_{1}\\x_{3}\end{bmatrix}}\qquad$

$={\begin{bmatrix}1&0&W^{0}&0\\0&1&0&W^{1}\\1&0&W^{2}&0\\0&1&0&W^{3}\end{bmatrix}}\,{\begin{bmatrix}{\begin{bmatrix}W^{0}&W^{0}\\W^{0}&W^{2}\end{bmatrix}}{\begin{bmatrix}x_{0}\\x_{2}\end{bmatrix}}\\{\begin{bmatrix}W^{0}&W^{0}\\W^{0}&W^{2}\end{bmatrix}}{\begin{bmatrix}x_{1}\\x_{3}\end{bmatrix}}\end{bmatrix}}\qquad \qquad \qquad (7)$

식(7) 오른쪽 행렬의 $W^{2}$ 는 N=2 DFT 행렬 안에 있으므로 식(3c)를 이용하여 $W_{4}^{2}\rightarrow W_{2}^{1}$ 로 변경해서 고쳐쓰면

$={\begin{bmatrix}1&0&W_{4}^{0}&0\\0&1&0&W_{4}^{1}\\1&0&W_{4}^{2}&0\\0&1&0&W_{4}^{3}\end{bmatrix}}\,{\begin{bmatrix}{\begin{bmatrix}W_{2}^{0}&W_{2}^{0}\\W_{2}^{0}&W_{2}^{1}\end{bmatrix}}{\begin{bmatrix}x_{0}\\x_{2}\end{bmatrix}}\\{\begin{bmatrix}W_{2}^{0}&W_{2}^{0}\\W_{2}^{0}&W_{2}^{1}\end{bmatrix}}{\begin{bmatrix}x_{1}\\x_{3}\end{bmatrix}}\end{bmatrix}}$

식(3b)를 이용하면 $W_{4}^{2}\rightarrow -W_{4}^{0}$, $W_{4}^{3}\rightarrow -W_{4}^{1}$, $W_{2}^{1}\rightarrow -W_{2}^{0}$ 이므로 다음과 같은 최종식을 얻는다.

${\begin{bmatrix}f_{0}\\f_{1}\\f_{2}\\f_{3}\end{bmatrix}}={\begin{bmatrix}1&0&W_{4}^{0}&0\\0&1&0&W_{4}^{1}\\1&0&-W_{4}^{0}&0\\0&1&0&-W_{4}^{1}\end{bmatrix}}\,{\begin{bmatrix}{\begin{bmatrix}W_{2}^{0}&W_{2}^{0}\\W_{2}^{0}&-W_{2}^{0}\end{bmatrix}}{\begin{bmatrix}x_{0}\\x_{2}\end{bmatrix}}\\{\begin{bmatrix}W_{2}^{0}&W_{2}^{0}\\W_{2}^{0}&-W_{2}^{0}\end{bmatrix}}{\begin{bmatrix}x_{1}\\x_{3}\end{bmatrix}}\end{bmatrix}}\qquad \qquad \qquad (8)$

따라서, N=4 인 DFT는 위 식(8)와 같이 N/2 즉,N=2 인 DFT 두 개를 사용해서 계산할 수 있고, 이는 다니엘슨-란초스 보조정리가 적용됨을 행렬을 이용해 보인 것이다. 그 결과로 입력과 출력이 서로 역비트순이 됨도 알 수 있다. 식(8)의 결과는 아래 그림(3)의 신호흐름도(signal flow graph)와 동일하며,나비도표(butterfly diagram)로도 불린다. 두 개의 N=2 DFT 나비를 볼 수 있다.

(그림 3) N=4 인 신호에대한 기수 2 시분할(DIT) FFT

식(5a)의 푸리에 행렬 F₄는 (7)식으로부터 다음과 같이 표현된다.

$\ F_{4}={\begin{bmatrix}I_{2}&T_{2}\\I_{2}&-T_{2}\end{bmatrix}}\,{\begin{bmatrix}F_{2}&O_{2}\\O_{2}&F_{2}\end{bmatrix}}C_{4}\qquad (9)$

여기서, I₂: 2x2 단위행렬, T₂: 회전자로 구성된 2x2 대각행렬, O₂: 2x2 Null 행렬 C₄: 4x4 column 치환행렬(permutation matrix)

식(9)에서 분할(decimation)을 수행하는 C₄가 식(7)에서 볼 수 없는 것은 이미 식(6)에서 시행했기 때문이다. (9)식을 값으로 표현하면 다음과 같다.

$\ F_{4}={\begin{bmatrix}1&0&1&0\\0&1&0&-i\\1&0&-1&0\\0&1&0&i\end{bmatrix}}\,{\begin{bmatrix}1&1&0&0\\1&-1&0&0\\0&0&1&1\\0&0&1&-1\end{bmatrix}}\,{\begin{bmatrix}1&0&0&0\\0&0&1&0\\0&1&0&0\\0&0&0&1\end{bmatrix}}={\begin{bmatrix}1&1&1&1\\1&-i&-1&i\\1&-1&1&i\\1&i&-1&-i\end{bmatrix}}\qquad (10)$

(9)식을 일반화하면, N=2m 일 때

$\ F_{2m}={\begin{bmatrix}I_{m}&T_{m}\\I_{m}&-T_{m}\end{bmatrix}}\,{\begin{bmatrix}F_{m}&O_{m}\\O_{m}&F_{m}\end{bmatrix}}C_{2m}\qquad (11)$

윗 식의 N/2-DFT 행렬 F_m 을 F₂ 가 될 때까지 재귀적으로(recursively) 분할하는 과정을 반복할 수 있다. 행렬 F_m 을 한번 더 분할할 수 있다면 다음과 같이 될 것이다.

${\begin{bmatrix}F_{m}&O_{m}\\O_{m}&F_{m}\end{bmatrix}}\,$$\ ={\begin{bmatrix}I_{m/2}&T_{m/2}&0&0\\I_{m/2}&-T_{m/2}&0&0\\0&0&I_{m/2}&T_{m/2}\\0&0&I_{m/2}&-T_{m/2}\end{bmatrix}}\,{\begin{bmatrix}F_{m/2}&0&0&0\\0&F_{m/2}&0&0\\0&0&F_{m/2}&0\\0&0&0&F_{m/2}\end{bmatrix}}C_{2m}\,$

그럼, (11)식의 N=2m 의 푸리에 행렬 F_2m을 2개의 N=m 인 행렬 F_m로 나누었을 때 연산 횟수의 변화를 살펴보도록 하자.

F_m 행렬 2 개는 2 m² 의 연산 횟수을 가지며, 식(11)의 맨 왼쪽 행렬은 m 개의 대각성분을 가지므로 연산횟수 m 개이다. C_2m 은 0,1 로 된 상수와 다름없어서 연산 횟수는 제로에 가깝다. 그래서, 다음과 같이 연산 횟수가 줄어든다.

$(2m)^{2}$ → $2(m^{2})+m$

$64^{2}$ → $2(32^{2})+32$

$32^{2}$ → $2(16^{2})+16$

$16^{2}$ → $2(8^{2})+8$

… …

결국 N=4,8,...,64 일 때의 연산횟수의 변화를 살펴보면 다음과 같으며, 맨 아래에 앞서 언급했던 연산횟수 식인 (12)식이 얻어진다.

$4^{2}$ → $2+2$

$8^{2}$ → $2(2+2)+4$

… $\qquad$ …

$64^{2}$ → $2(2(2(2(2+2)+4)+8)+16)+32$

――――――――――――――――――――――――

$N^{2}$ → $(N/2)log_{2}N\qquad \qquad \qquad (12)$

이 (12)식은 아래 (표.1)에서와 같은 방식으로도 유도될 수도 있다.

2. 기수 2 주파수분할 FFT(Radix-2 Decimation in Frequency FFT)

Partition

어떤 길이 N인 수열을 다음과 같이 단순히 index 순서를 기준으로 처음 절반과 나중 절반으로 나눈 것을 Partition 이라 한다.

$(x_{0},x_{1},x_{2},...,x_{N-2},x_{N-1})$→$(x_{0},x_{1},x_{2},...,x_{{\frac {N}{2}}-2},x_{{\frac {N}{2}}-1})+(x_{\frac {N}{2}},x_{{\frac {N}{2}}+1},x_{{\frac {N}{2}}+2},...,x_{N-2},x_{N-1})$

이론

앞서 언급했듯이 주파수분할은 출력신호를 나누는 방식이다. 위의 partition 방식으로 (2)식의 DFT를 나누면 다음과 같다.

$f_{k}=\sum _{n=0}^{N-1}x_{n}W_{N}^{kn}=\sum _{n=0}^{N/2-1}x_{n}W_{N}^{kn}+\sum _{n=N/2}^{N-1}x_{n}W_{N}^{kn}$

$=\sum _{n=0}^{N/2-1}x_{n}W_{N}^{kn}+W_{N}^{kN/2}\sum _{n=0}^{N/2-1}x_{n+N/2}W_{N}^{kn}$

윗 식 둘째 항은 DFT가 아니다. 그런데, $W_{N}^{kN/2}=(-1)^{k}$ 이므로,

$f_{k}=\sum _{n=0}^{N/2-1}x_{n}W_{N}^{kn}+(-1)^{k}\sum _{n=0}^{N/2-1}x_{n+N/2}W_{N}^{kn}=\sum _{n=0}^{N/2-1}{\Bigl [}x_{n}+(-1)^{k}x_{n+N/2}{\Bigr ]}W_{N}^{kn}$

이제 $0\leq n\leq {\frac {N}{2}}-1$ 에 대하여

$g_{n}=x_{n}+x_{n+N/2},\quad h_{n}=(x_{n}-x_{n+N/2})W_{N}^{n}\quad \qquad (13)$

이라고 하자. 그러면 짝수 $k=2r$ 에 대하여

$f_{2r}=\sum _{n=0}^{N/2-1}g_{n}W_{N}^{2rn}=\sum _{n=0}^{N/2-1}g_{n}W_{N/2}^{rn}\qquad \qquad (14a)$

이고, 홀수 $k=2r+1$ 에 대하여

$f_{2r+1}=\sum _{n=0}^{N/2-1}h_{n}W_{N}^{2rn}=\sum _{n=0}^{N/2-1}h_{n}W_{N/2}^{rn}\qquad \qquad (14b)$

여기서 $W_{N}^{2}=W_{N/2}$ 를 이용하였다. 이제 위의 두 식은 각각 길이가 N/2인 $g_{n}$과 $h_{n}$의 이산 푸리에 변환(DFT)이 되었다.

위와 같은 방식으로 다음과 같이 나눌 수 있다.

$g_{1n}=g_{n}+g_{n+N/4},\quad g_{2n}=(g_{n}-g_{n+N/4})W_{N/2}^{n}$

$h_{1n}=h_{n}+h_{n+N/4},\quad h_{2n}=(h_{n}-h_{n+N/4})W_{N/2}^{n}$

이를 길이가 2인 신호들이 얻어질 때까지 반복한다(분할정복 알고리즘). 기본적인 주파수 분할 나비는 식(13)으로부터 다음과 같이 구성된다.

기수 2 주파수분할 나비

다음 식(15)은 N=4일때 기수 2 주파수분할(DIF) FFT이다. 식(6)의 시분할(DIT) 방식과 달리 출력 f 의 원소들이 분할된 것을 볼 수 있다.

${\begin{bmatrix}f_{0}\\f_{2}\\f_{1}\\f_{3}\end{bmatrix}}={\begin{bmatrix}W^{0}&W^{0}&0&0\\W^{0}&W^{2}&0&0\\0&0&W^{0}&W^{1}\\0&0&W^{0}&W^{3}\end{bmatrix}}{\begin{bmatrix}1&0&W^{0}&0\\0&1&0&W^{0}\\1&0&W^{2}&0\\0&1&0&W^{2}\end{bmatrix}}\,{\begin{bmatrix}x_{0}\\x_{1}\\x_{2}\\x_{3}\end{bmatrix}}$

$\quad ={\begin{bmatrix}W^{0}&W^{0}&0&0\\W^{0}&-W^{0}&0&0\\0&0&W^{0}&W^{1}\\0&0&W^{0}&-W^{1}\end{bmatrix}}{\begin{bmatrix}1&0&W^{0}&0\\0&1&0&W^{0}\\1&0&-W^{0}&0\\0&1&0&-W^{0}\end{bmatrix}}\,{\begin{bmatrix}x_{0}\\x_{1}\\x_{2}\\x_{3}\end{bmatrix}}\qquad (15)$

다음은 (15)식을 표현한 신호흐름도이다. 시분할 FFT의 N=2 DFT 는 (그림 3)과 달리 마지막 stage에 있으며 이는 행렬(15)에서도 확인할 수 있다.

N=4 인 신호에대한 주파수분할 FFT의 나비도표

공통인수 FFT(Common-Factor FFT)

이와 대응하는 상반된 개념은 소인수 알고리즘 즉 PFA(Prime-Factor Algorithm)이라 할 수 있다.

1. 기수 r FFT (Radix-r FFT)

지금까지 논의에서는 기수 2(radix-2) FFT 만을 다루었다. 즉 신호 수 N=2 ⁿ 인 경우이다. 그러나 N 은 다양한 기수를 가질 수 있다. 즉, 일반적으로 기수 r 일 때, N=r ⁿ이다. 그림(5)를 보면서 다음의 특징들을 살펴보십시오

1)각 stage에서 볼 때, 사용하는 N/r-point DFT(즉 butterfly)는 동일하다.

2) 기수 r FFT에서 나비(butterfly)는 한 stage에서 N/r 개가 된다.

(예) N=16 이고 기수 r=2 이면 stage당 나비 수는 16/2=8 개.

3) FFT 나비도표(butterfly diagram)에서 stage의 수 s는

$s=log_{r}N$

(예)

N=16 이고 기수 r=2 이면 stage 수 $s=log_{2}16=log_{2}2^{4}=4$이다.

N=16 이고 기수 r=4 이면 stage 수 $s=log_{4}16=log_{4}4^{2}=2$이다.

결국, 기수 4의 stage수 s는 기수 2 일 때 보다 반으로 감소한다.

$\because \quad s=log_{4}N=log_{2}N/log_{2}4=log_{2}N/2$

그림 5. 기수 4와 기수 2인 경우에 stage수와 butterfly수/stage 비교

4) 기수에 따른 계산 효율

• 대체적으로 기수 4는 대규모 변환에서 기수 2보다 약 20% 더 효율적이다(아래 표 참고).

표 1. 기수 r 에 따른 복소수 곱셈 수

기수	stage 수 (a)	stage 내의 나비(butterfly)수 (b)	butterfly 연산 안에서의 복소수 곱셈의 수 (c)	전체 복소수 곱셈의 수 (a)x(b)x(c)
radix-2	log2N	N/2	1 (회전자 계수 중복으로 1/2만 연산)	(N/2)log2N
radix-4	log4N=(1/2)log2N	N/4	3	(3N/8)log2N

표의 결과만을 볼 때 radix-4 가 고속이 될 것 같지만, 실제로는 실장 방법에 의존한다. CPU로 실행하는 경우는 연산기가 1개 이므로, 연산 횟수가 적은 radix-4 쪽이 좀 더 고속이 될 것 같다. 그러나, FPGA로 구현하는 경우의 트레이드 오프는 복잡하다. 나비연산을 복수의 연산기를 사용해 파이프라인으로 하면 radix-2와 radix-4 모두 같은 1 클럭이 되어버리므로, radix-4가 4배 고속이 될 것 같아 보이지만 1개의 메모리에서는 1개씩 밖에 데이터를 읽을 수 없기 때문에 여기가 병목이 되어, radix-4는 4 클럭, radix-2는 2 클럭이 소요되어, radix-4는 radix-2의 2배 밖에 빨라지지 않는다. 그래서, 보통은 메모리를 분할하여 복수의 데이터를 동시에 읽도록 고안하지만, 메모리 읽기를 병렬화하는 것에 맞추어, 연산기도 radix-2의 나비 연산을 4병렬화시키면, 연산기 수도 radix-4와 거의 같아져 차이가 없어진다. 그렇게 생각하면, radix-2 쪽이 병렬화 설계가 더 간단한 만큼, 유리할지도 모른다.

• 때때로 기수 8이 사용되지만 더 긴 기수 나비는 일반적이지 않다. 효율성 개선이 크지않고, 추가된 복잡성이 상당하기 때문이다(특히 하드웨어 구현의 경우).

2. 혼합 기수 FFT (Mixed-radix FFT)

혼합된 기수 고속 푸리에 변환 알고리즘은 1969년 Richard C. Singleton에 의하여 처음 고안된 알고리즘으로 쿨리-튜키 알고리즘의 변종으로 신호의 길이가 임의의 합성수(composite number)일 때도 대응할 수 있는 보다 일반적인 알고리즘이다. 특징은 다음과 같다.

1) N≠r ⁿ 일 때도 포함해서 일반적으로 합성수 N=N₁N₂ ...N_m 일 때 사용한다.

N=2 ⁿ(즉, 기수 2)인 경우에 N=2, 4, 8,...이고, N=4 ⁿ 인 경우 N=4, 16, 64... 등의 갯수를 가질 수 있으나, 중간 갯수인 N=12 , N=15, N=24 등은 단일한 기수로는 담당할 수 없게 된다. 따라서 서로 다른 기수가 혼합된 여러 기수를 사용해 해결하는 것이다.

N=2 ⁿ×q ^s (q=2 ^m>2, 1<s≤m ) 이라면 기수 2 알고리즘 n 단계를 변환의 시작 또는 끝에 적용하고 나머지는 기수 q 알고리즘 s 단계를 수행한다.

예) N=2 ^(2m+1)=2 ¹×4 ^m

N=32=2 ¹×4 ² 즉 기수2 한 단계와 기수4 두 단계를 수행한다.

2) 나비도표에 서로 다른 기수의 나비들이 존재한다.

이론

식(5)에서 DFT는 푸리에 행렬 T와 벡터 x, f 로서 표현됨을 알았다. 즉

$f=Tx$

T의 행렬요소는 $t_{jk}=e^{(i2\pi jk/n)}$이다.

신호수가 합성수 N이며, 다음과 같은 m개의 인수로 분해(decomposition)될 때

$N=\prod _{i=1}^{m}N_{i}$

이고, N X N 행렬 T 는 다음과 같이 행렬곱으로 표현된다.

$T=PF_{m}F_{m-1}F....F_{2}F_{1}$

여기서, P는 치환행렬(permutation matrix), F_i는 n의 계수 n_i에 해당하는 변환 단계(transform step) 행렬이다.

행렬 F_i는 각 행과 열에 n_i개의 0이 아닌 요소만 있고 차원 n_i의 n/n_i 개의 square 부분행렬(submatrix)들로 분할될 수 있다. FFT 알고리즘의 기초가 되는 것이 바로 이 분할과 그에 따른 곱셈의 감소이다. 행렬 F_i는 다음과 같이 더 분해될 수 있다.

$F_{i}=R_{i}T_{i}$

여기서 R_i는 회전자의 대각행렬이고 T_i는 n_i 개의 동일한 square 부분행렬로 분할될 수 있으며 각 행렬은 차원 n_i의 복소 푸리에 변환이다.

신호의 길이가 합성수 N=N₁N₂ 이라고 하자. 두 수의 곱으로 표현되므로 m=2 인 경우에 해당한다. 이 신호의 DFT는 식(2)로부터

$f_{k}=\sum _{n=0}^{N_{1}N_{2}-1}x_{n}W_{N_{1}N_{2}}^{kn}\qquad \qquad k=0,\dots ,N-1\qquad \qquad \qquad (16)$

신호 x 를 N₁ x N₂ 이차원 배열로 재배치하는 것을 생각해 볼 수 있다.

N₁ =3, N₂ =8 인 신호를 고려하자.

$h={\begin{bmatrix}x_{0}&x_{1}&x_{2}&x_{3}&x_{4}&x_{5}&x_{6}&x_{7}\\x_{8}&x_{9}&x_{10}&x_{11}&x_{12}&x_{13}&x_{14}&x_{15}\\x_{16}&x_{17}&x_{18}&x_{19}&x_{20}&x_{21}&x_{22}&x_{23}\\\end{bmatrix}}\qquad \qquad \qquad (17)$

으로 간주한다. $n_{1},n_{2}$ 를 각각 위 (17)식 h의 행 번호, 열 번호, $k_{1},k_{2}$ 를 각각 아래 (18)식 H의 행 번호, 열 번호이라 할 때, $n=n_{1}N_{2}+n_{2},\quad k=k_{2}N_{1}+k_{1}$ 이라 하자. 그러면, 위의 이차원 배열을

$h(n_{1},n_{2})=x(n_{1}N_{2}+n_{2})=x_{n}\qquad (0\leq n_{1}\leq N_{1}-1,\quad 0\leq n_{2}\leq N_{2}-1)$

와 같이 쓸 수 있다. 배열 h 에 대응하는 DFT는

$H(k_{1},k_{2})=f(k_{2}N_{1}+k_{1})=f_{k}\qquad (0\leq k_{1}\leq N_{1}-1,\quad 0\leq k_{2}\leq N_{2}-1)$

을 생각할 수 있다. 위의 DFT를 이차원 배열로 나타내면

$H={\begin{bmatrix}f_{0}&f_{3}&f_{6}&f_{9}&f_{12}&f_{15}&f_{18}&f_{21}\\f_{1}&f_{4}&f_{7}&f_{10}&f_{13}&f_{16}&f_{19}&f_{22}\\f_{2}&f_{5}&f_{8}&f_{11}&f_{14}&f_{17}&f_{20}&f_{23}\\\end{bmatrix}}\qquad \qquad \qquad (18)$

두 이차원 배열을 구성한 순서가 (17)식은 행우선순서(row-major order), 출력에 해당하는 (18)식은 열우선순서(column-major order)로 다름을 알 수 있고 이를 전치(transpose)되었다고 한다. 이러한 모든 전치의 최종 결과는 입력(DIF) 또는 출력(DIT) 인덱스의 비트 반전에 해당한다. 그 이유는 잠시 뒤 명확해 질 것이다.

그 전에 위에 언급한 $n_{1},n_{2}$의 범위로부터 (16)식과 같이 $n$의 범위가 $0\leq n\leq N_{1}N_{2}-1$ 이 되는 지 확인해 보자.

<증명>

$n=n_{1}N_{2}+n_{2}$ 이라할 때

$n$의 최소값은 $n_{1}=n_{2}=0$ 일 때 $n=0$이며,

$n$의 최대값은 $n_{1}=N_{1}-1\quad and\quad n_{2}=N_{2}-1$ 일 때 이므로

$n=n_{1}N_{2}+n_{2}=(N_{1}-1)N_{2}+(N_{2}-1)=N_{1}N_{2}-1\quad$이다.

따라서 $0\leq n\leq N_{1}N_{2}-1$이다. <증명 끝>

이제 원래의 DFT (16)식을 다음과 같이 쓸 수 있다.

$H(k_{1},k_{2})=\sum _{n_{1}=0}^{N_{1}-1}\sum _{n_{2}=0}^{N_{2}-1}h(n_{1},n_{2})W_{N_{1}N_{2}}^{(k_{2}N_{1}+k_{1})(n_{1}N_{2}+n_{2})}$

$=\sum _{n_{1}=0}^{N_{1}-1}\sum _{n_{2}=0}^{N_{2}-1}h(n_{1},n_{2})W_{N_{1}N_{2}}^{k_{2}n_{1}N_{1}N_{2}}W_{N_{1}}^{k_{1}n_{1}}W_{N_{2}}^{k_{2}n_{2}}W_{N_{1}N_{2}}^{k_{1}n_{2}}$

$=\sum _{n_{1}=0}^{N_{1}-1}\sum _{n_{2}=0}^{N_{2}-1}h(n_{1},n_{2})W_{N_{1}}^{k_{1}n_{1}}W_{N_{2}}^{k_{2}n_{2}}W_{N_{1}N_{2}}^{k_{1}n_{2}}\quad (\because \ \ W_{N_{1}N_{2}}^{k_{2}n_{1}N_{1}N_{2}}=e^{-2\pi i(k_{2}n_{1})}=1)\qquad (19)$

윗 식 세 개의 회전자 곱에서 처음 두 회전자들 때문에 이차원 DFT를 연상하게 한다. 마지막 회전자는 이들 전체 이차원 DFT에 대한 일종의 회전자이다. 이 사실은 조금 전에 두 이차원 배열을 구성한 순서를 다르게 한 즉, 전치(transpose)를 만든 이유이다.

윗 식을 다시 쓰면

$H(k_{1},k_{2})=\sum _{n_{1}=0}^{N_{1}-1}W_{N_{1}N_{2}}^{k_{1}n_{2}}{\Biggl [}\sum _{n_{2}=0}^{N_{2}-1}h(n_{1},n_{2})W_{N_{2}}^{k_{2}n_{2}}{\Biggr ]}W_{N_{1}}^{k_{1}n_{1}}$

괄호 안은 배열의 신호수 $N_{2}$ 개의 각 열(column)에 대한 DFT이며, 이를 $h^{'}(n_{1},n_{2})$ 라 하자. 그러면,

$H(k_{1},k_{2})=\sum _{n_{1}=0}^{N_{1}-1}W_{N_{1}N_{2}}^{k_{1}n_{2}}h^{'}(n_{1},n_{2})W_{N_{1}}^{k_{1}n_{1}}\qquad \qquad (20)$

이는 먼저 회전자를 곱한 후에 계산한 신호수 $N_{1}$ 개의 행들에 대한 DFT 식이다. 이것과 이차원 DFT와의 유일한 차이점은 회전자에 의한 추가적인 곱셈이 있다는 것이다.

만약 합성수였던 원래 신호의 길이 N=N₁N₂에서 N₁ 이나 N₂가 다시 합성수라면, 분해(decomposition)의 기본 아이디어를 재귀적으로 적용한다.

일반적인 분해(decomposition)를 위한 쿨리-튜키 FFT의 기본 단계는 1d(1차원) DFT를 2d(2차원) DFT와 동일한 것으로 재해석하는 것으로 볼 수 있다. 길이가 N = N1N2인 1d 입력 배열은 열 우선 순서(column-major order)로 저장된 2d N1×N2 행렬로 재해석된다. 처음에는 N2 방향(비연속 방향)을 따라 더 작은 1d DFT를 수행한 다음 회전자를 곱하고 마지막으로 N1 방향을 따라 1d DFT를 수행한다. 전치(transpose) 단계는 여기에 표시된 것처럼 중간에 수행하거나 시작 또는 끝에서 수행할 수 있다. 이 과정을 더 작은 변환에 대해 재귀적으로 수행한다.

지금까지의 과정을 단계별로 정리하면 다음과 같다.

1. N₂ 크기의 N₁ DFT를 수행한다.

2. 회전자를 곱한다.

3. N₁ 크기의 N₂ DFT를 수행한다.

일반적으로 N₁ 또는 N₂는 기수이다. N₁이 기수이면 DIT(Decimation in Time) 알고리즘이라고 하고, N₂가 기수이면 DIF(Sande-Tukey 알고리즘이라고도 함)라고 한다.

과정의 역변환을 만드는 것은 그리 어렵지않다. 각 행의 IDFT(역이산 푸리에 변환)을 구한 다음, 결과 배열에 회전자의 공액(conjugate)을 곱한다. 마지막으로 각 열의 IDFT를 구한다.

그 외의 알고리즘

• Prime Factor Algorithm (PFA)
• Bruun's FFT algorithm
• Rader's FFT algorithm
• Bluestein's FFT algorithm

응용

• 스펙트럼 분석기
• OFDM 변복조기
• CT 스캐너, MRI 등
• MP3 압축방식

참고 문헌

• G. C. Danielson and C. Lanczos, "Some Improvements in Practical Fourier Analysis and Their Application to X-ray Scattering form Liquids", J. Frank. Inst., vol. 233, 4 & 5, 365-380 & 435-452(1942).
• J. W. Cooley and J. W. Tukey, Math. of Comput. 19, 90, 297 (1965).
• Richard C. Singleton, "An algorithm for computing the mixed radix fast Fourier transform", IEEE Trans. Audio Electroac. vol. 17, 2, 93-103(1969).
• Carl Friedrich Gauss, "Nachlass: Theoria interpolationis methodo nova tractata", Werke band 3, 265–327 (Königliche Gesellschaft der Wissenschaften, Göttingen, 1866). See also M. T. Heideman, D. H. Johnson, and C. S. Burrus, "Gauss and the history of the fast Fourier transform," IEEE ASSP Magazine vol. 1, 4, 14–21(1984).
• H. V. Sorensen, D. L. Jones, M. T. Heideman, and C. S. Burrus, "Real-valued fast Fourier transform algorithms," IEEE Trans. Acoust. Speech Sig. Processing ASSP-35, 849–863 (1987).
• 최행진, "물리학자 푸리에와 고속 푸리에 변환", 교우사 (2007).
• zakii, 高速フーリエ変換 Fast Fourie Transform (FFT), http://zakii.la.coocan.jp/fourie/31_fft_algorithm.htm
• 장준원, Fourier Matrix & Fast Fourier Transform, https://blog.naver.com/kaoara/221905762357
• Cooley-Tukey FFT algorithm, Wikipedia.org https://en.wikipedia.org/wiki/Cooley-Tukey_FFT_algorithm