고속 푸리에 변환

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

푸리에 해석(Fourier Analysis)은 신호를 원래의 영역(주로 시간 또는 공간)에서 주파수 영역으로 변환하거나 그 반대로 변환하는 과정이다. 이산 데이터에서 사용하는 DFT는 값의 시퀀스를 서로 다른 주파수 성분으로 분해하여 얻는다. 이 작업은 다양한 분야에서 유용하지만, 정의를 기반으로 직접 계산하는 방식은 너무 느려 실용적이지 않다. FFT는 DFT 행렬을 희소 행렬로 인수분해하여 이러한 변환을 빠르게 계산한다. 이를 통해 DFT 계산 복잡도를 정의를 직접 적용했을 때 ${\textstyle O(n^{2})}$ 에서 ${\textstyle O(n\log n)}$ 으로 줄어든다 (여기서 $n$ 은 데이터 크기). 이러한 속도 차이는 데이터 세트가 수천 또는 수백만 개에 이를 경우 특히 크다. 반올림 오류가 있는 경우에도, 많은 FFT 알고리즘은 DFT 정의를 직접 또는 간접적으로 사용하는 것보다 훨씬 더 정확하다. FFT 알고리즘은 복소수 산술부터 군론 및 정수론에 이르기까지 다양한 이론을 기반으로 한 여러 형태가 존재한다.

FFT는 공학, 음악, 과학, 수학 등 다양한 분야에서 널리 활용된다. 기본 개념은 1965년에 대중화되었으나, 일부 알고리즘은 카를 프리드리히 가우스에 의해 이미 1805년에 도출되었다.^[1]^[2] 그 뒤에도 제한된 형태의 FFT가 종종 발견되었음이 밝혀졌다. 1994년, 길버트 스트랭(Gilbert Strang)은 FFT를 "우리 시대에서 가장 중요한 수치 알고리즘"이라고 묘사했으며, IEEE 잡지 Computing in Science & Engineering에서 선정한 "20세기 최고의 알고리즘 10선"에 포함되었다.

가장 잘 알려진 FFT 알고리즘은 n의 소인수분해에 의존하지만, 모든 $n$ 에 대해 $O(n\log n)$ 복잡도를 갖는 FFT도 존재한다. 많은 FFT 알고리즘은 ${\textstyle e^{-2\pi i/n}}$ 이 $n$ 차 원시 단위근이라는 사실만을 활용하며, 이는 수론적 변환과 같은 유한체에서의 유사 변환에도 적용할 수 있다. 역 DFT는 지수에서 부호가 반대이고 $1/ n$ 계수가 추가된 형태로 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)

이 식을 정의에 따라 계산하면 ${\textstyle O(N^{2})}$ 의 연산이 필요하지만, FFT를 이용하면 ${\textstyle O(N\log N)}$ 의 연산만으로 가능하다. $N$ 이 합성수일 경우 그 소인수 분해를 이용하여 연산 횟수를 줄일 수는 있지만, FFT를 사용하면 $N$ 이 소수일 경우에도 ${\textstyle 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의 기초가 되는 공식을 발견한다.^[3] 즉, 길이 N(단 짝수인 경우)의 DFT는 각각 길이가 N/2인 두 개의 DFT의 합으로 다시 쓸 수 있다는 것. 하나는 짝수 index 포인트들로 형성되고, 다른 것은 홀수 index 포인트들로 만들어진다. 이 분할(decimation)을 수식으로 표현해 보자.^[4]

식(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)라고 한다.

역비트순(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과 같은 순서가 된다.

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

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

가장 일반적으로 사용되는 FFT 알고리즘은 쿨리-튜키 알고리즘이다. 이 방법은 1965년에 J. W. 쿨리와 J. W. 튜키가 발표하여 널리 알려졌다.^[5]

식(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)을 의미한다).

$={\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}}$

또한 식(3a)로 부터 $W^{4}=W^{0},\quad W^{6}=W^{2+4}=W^{2+0}=W^{2}$ 의 관계를 얻으므로

={\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)

={\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 나비를 볼 수 있다.

식(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 로 된 상수와 다름없어서 연산 횟수는 제로에 가깝다. 그래서, 다음과 같이 연산 횟수가 줄어든다.^[6]

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

$64^{2}=4096$ → $2(32^{2})+32=2080$

$32^{2}=1024$ → $2(16^{2})+16=528$

$16^{2}=256$ → $2(8^{2})+8=136$

… …

결국 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)으로부터 다음과 같이 구성된다.

다음 식(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)에서도 확인할 수 있다.

공통인수 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$

4) 기수에 따른 계산 효율^[7]

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

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

표의 결과만을 볼 때 radix-4 가 고속이 될 것 같지만, 실제로는 실장 방법에 의존한다. CPU로 실행하는 경우는 연산기가 1개 이므로, 연산 횟수가 적은 radix-4 쪽이 좀 더 고속이 될 것 같다. 그러나, FPGA로 구현하는 경우의 상충관계(trade-off)는 복잡하다. 나비연산을 복수의 연산기를 사용해 파이프라인으로 하면 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)일 때도 대응할 수 있는 보다 일반적인 알고리즘이다.^[8] 특징은 다음과 같다.

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)의 기본 아이디어를 재귀적으로 적용한다.

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

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

2. 회전자를 곱한다.

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

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

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

예) 위 (16)식에서 본 길이 N=24인 신호의 이산 푸리에 변환을 생각해 보자.

$f_{p}=\sum _{k=0}^{23}f(k)e^{{-i2\pi p \over 24}k}=\sum _{k=0}^{N-1}f(k)W_{24}^{pk}\qquad (0\leq p\leq 23)$

위의 식을 다음과 같이 두 개의 12 point data blocks 으로 쓸 수 있다.

$f_{p}=\sum _{m=0}^{1}\sum _{n=0}^{11}f(2n+m)W_{24}^{(2n+m)p}$

또한 위의 12 point block을 두 개의 6 point data blocks 으로 쓸 수 있다.

$f_{p}=\sum _{s=0}^{1}\sum _{m=0}^{1}\sum _{n=0}^{5}f(4n+2m+s)W_{24}^{(4n+2m+s)p}$

마지막으로 6 point block은 두 개의 3 point data blocks 으로 쓸 수 있다.

$f_{p}=\sum _{t=0}^{1}\sum _{s=0}^{1}\sum _{m=0}^{1}\sum _{n=0}^{2}f(8n+4m+2s+t)W_{24}^{(8n+4m+2s+t)p}$

위의 DFT가 바로 길이 24인 신호의 혼합 기수 FFT 알고리즘을 위한 식이다. 한 단계(stage)의 기수3 DFT와 세 단계의 기수2 DFT를 볼 수 있다.

그 외의 알고리즘

Prime Factor Algorithm (Good-Thomas algorithm)
Bruun's FFT algorithm
Winograd FT algorithm
Rader's FFT algorithm
Bluestein's FFT algorithm (chirp-z algorithm)
hexagonal FFT algorithm
Bailey's FFT algorithm

응용

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

같이 보기

참고 문헌

↑ Heideman, Michael T.; Johnson, Don H.; Burrus, Charles Sidney (1984). “Gauss and the history of the fast Fourier transform” (PDF). 《IEEE ASSP Magazine》 1 (4): 14–21. CiteSeerX 10.1.1.309.181. doi:10.1109/MASSP.1984.1162257. S2CID 10032502. 2013년 3월 19일에 원본 문서 (PDF)에서 보존된 문서.
↑ Gauss, Carl Friedrich (1866). “Nachlass: Theoria interpolationis methodo nova tractata” (PDF). 《Königliche Gesellschaft der Wissenschaften, Göttingen》. Werke band 3: 265–327. 2006년 3월 15일에 원본 문서에서 보존된 문서.
↑ Danielson, G. C.; Lanczos, C (1945). “Some Improvements in Practical Fourier Analysis and Their Application to X-ray Scattering form Liquids”. 《J. Frank. Inst.》 233 (4 & 5): 365–380 & 435–452.
↑ 최, 행진 (2007). 《물리학자 푸리에와 고속 푸리에 변환》. 교우사.
↑ Cooley, James W.; Tukey, John W. (1965). “An algorithm for the machine calculation of complex Fourier series”. 《Math. Comput.》 19 (90): 297–301. doi:10.2307/2003354. JSTOR 2003354.
↑ 장, 준원. “Fourier Matrix & Fast Fourier Transform”.
↑ zakii. “高速フーリエ変換 Fast Fourie Transform (FFT)”.
↑ Singleton, Richard C. (1969). “An algorithm for computing the mixed radix fast Fourier transform”. 《IEEE Trans. Audio Electroac.》 17 (2): 93–103.
↑ wikipedia.org. “Cooley–Tukey FFT algorithm”.

[Heideman_Johnson_Burrus_1984-1] Heideman, Michael T.; Johnson, Don H.; Burrus, Charles Sidney (1984). “Gauss and the history of the fast Fourier transform” (PDF). 《IEEE ASSP Magazine》 1 (4): 14–21. CiteSeerX 10.1.1.309.181. doi:10.1109/MASSP.1984.1162257. S2CID 10032502. 2013년 3월 19일에 원본 문서 (PDF)에서 보존된 문서.

[Gauss_1866-2] Gauss, Carl Friedrich (1866). “Nachlass: Theoria interpolationis methodo nova tractata” (PDF). 《Königliche Gesellschaft der Wissenschaften, Göttingen》. Werke band 3: 265–327. 2006년 3월 15일에 원본 문서에서 보존된 문서.

[Danielson_Lanczos_1942-3] Danielson, G. C.; Lanczos, C (1945). “Some Improvements in Practical Fourier Analysis and Their Application to X-ray Scattering form Liquids”. 《J. Frank. Inst.》 233 (4 & 5): 365–380 & 435–452.

[최_2007-4] 최, 행진 (2007). 《물리학자 푸리에와 고속 푸리에 변환》. 교우사.

[Cooley_Tukey_1965-5] Cooley, James W.; Tukey, John W. (1965). “An algorithm for the machine calculation of complex Fourier series”. 《Math. Comput.》 19 (90): 297–301. doi:10.2307/2003354. JSTOR 2003354.

[장_2020-6] 장, 준원. “Fourier Matrix & Fast Fourier Transform”.

[zakii_31-7] zakii. “高速フーリエ変換 Fast Fourie Transform (FFT)”.

[Singleton_1969-8] Singleton, Richard C. (1969). “An algorithm for computing the mixed radix fast Fourier transform”. 《IEEE Trans. Audio Electroac.》 17 (2): 93–103.

[Cooley-Tukey_FFT-9] wikipedia.org. “Cooley–Tukey FFT algorithm”.

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]