고속 푸리에 변환(高速 푸리에 變換, 영어: Fast Fourier Transform, FFT)은 이산 푸리에 변환(Discrete Fourier Transform, DFT)과 그 역변환을 빠르게 수행하는 효율적인 알고리즘이다. FFT는 디지털 신호 처리에서 편미분 방정식의 근을 구하는 알고리즘에 이르기까지 많은 분야에서 사용한다.
푸리에 해석(Fourier Analysis)은 신호를 원래의 영역(주로 시간 또는 공간)에서 주파수 영역으로 변환하거나 그 반대로 변환하는 과정이다. 이산 데이터에서 사용하는 DFT는 값의 시퀀스를 서로 다른 주파수 성분으로 분해하여 얻는다. 이 작업은 다양한 분야에서 유용하지만, 정의를 기반으로 직접 계산하는 방식은 너무 느려 실용적이지 않다. FFT는 DFT 행렬을 희소 행렬로 인수분해하여 이러한 변환을 빠르게 계산한다. 이를 통해 DFT 계산 복잡도를 정의를 직접 적용했을 때 에서 으로 줄어든다 (여기서 n은 데이터 크기). 이러한 속도 차이는 데이터 세트가 수천 또는 수백만 개에 이를 경우 특히 크다. 반올림 오류가 있는 경우에도, 많은 FFT 알고리즘은 DFT 정의를 직접 또는 간접적으로 사용하는 것보다 훨씬 더 정확하다. FFT 알고리즘은 복소수 산술부터 군론 및 정수론에 이르기까지 다양한 이론을 기반으로 한 여러 형태가 존재한다.
FFT는 공학, 음악, 과학, 수학 등 다양한 분야에서 널리 활용된다. 기본 개념은 1965년에 대중화되었으나, 일부 알고리즘은 카를 프리드리히 가우스에 의해 이미 1805년에 도출되었다.[1][2] 그 뒤에도 제한된 형태의 FFT가 종종 발견되었음이 밝혀졌다. 1994년, 길버트 스트랭(Gilbert Strang)은 FFT를 "우리 시대에서 가장 중요한 수치 알고리즘"이라고 묘사했으며, IEEE 잡지 Computing in Science & Engineering에서 선정한 "20세기 최고의 알고리즘 10선"에 포함되었다.
가장 잘 알려진 FFT 알고리즘은 n의 소인수분해에 의존하지만, 모든 n에 대해 복잡도를 갖는 FFT도 존재한다. 많은 FFT 알고리즘은 이 n차 원시 단위근이라는 사실만을 활용하며, 이는 수론적 변환과 같은 유한체에서의 유사 변환에도 적용할 수 있다. 역 DFT는 지수에서 부호가 반대이고 1/n 계수가 추가된 형태로 DFT와 동일하기 때문에, 모든 FFT 알고리즘은 역변환에도 쉽게 적용할 수 있다.
이 식을 정의에 따라 계산하면 의 연산이 필요하지만, FFT를 이용하면 의 연산만으로 가능하다. N이 합성수일 경우 그 소인수 분해를 이용하여 연산 횟수를 줄일 수는 있지만, FFT를 사용하면 N이 소수일 경우에도 번의 연산 횟수를 보장한다.
을 회전자(Twiddle factor)라고 부른다. 1의 거듭제곱근(Root of unity), 단위근도 회전자의 다른 호칭이다. 이 회전자를 사용해 식(1)을 다시 정리하면,
회전자는 N에 관한 주기성으로 인해 다음과 같은 성질이 있음을 쉽게 알수있다.
n:자연수
(예 1)
쿨리-튜키 알고리즘
분할(Decimation)
어떤 길이 N인 수열을 다음과 같이 index가 짝수인 것과 홀수인 것들만 각각 모아서 두 개의 수열로 나누는 것이다.
→
분할 대신 일정 부분을 뽑아낸다는 '솎음'으로 번역한 경우도 있는데 단어 번역으로는 더 정확하지만, FFT에서는 오히려 분할이라는 표현이 이해하는데 더 낫다.
다니엘슨-란초스 보조정리(Danielson-Lanczos Lemma)
1942년 위 두사람은 FFT의 기초가 되는 공식을 발견한다.[3] 즉, 길이 N(단 짝수인 경우)의 DFT는 각각 길이가 N/2인 두 개의 DFT의 합으로 다시 쓸 수 있다는 것. 하나는 짝수 index 포인트들로 형성되고, 다른 것은 홀수 index 포인트들로 만들어진다.
이 분할(decimation)을 수식으로 표현해 보자.[4]
식(2)의 DFT를 짝수번째 항, 홀수 번째 항끼리 묶으면 즉, 분할하면
길이가 짝수이므로 각 항의 갯수는 같다. 두번째 괄호 항에서 을 밖으로 빼면
이를 다시 쓰면
그런데,
가 성립하므로,
여기서, 는 위 수식에서 보였듯이 각각 짝수(Even)항과 홀수(Odd)항으로 구성된 N/2-DFT이다.
다음으로 fk+N/2를 구해보자. 윗식으로부터
from (3a),(3b)
fk+N/2 의 실제 범위는 fk 와 달리 k에 N/2가 추가되었으므로 임에 주의하자. 정리하면,
식(4)의 중요성은 다음 쿨리-튜키 알고리즘에서 설명하겠다. 이 식을 그림(1)로 나타낼 수 있으며, 마치 나비 모양이라고 해서 이런 분할 과정을 나비연산(butterfly operation)이라 부르고, 그 그림을 나비(butterfly)라고 한다.
역비트순(Bit-reversed Order)
만약 N이 2의 거듭제곱, 즉 N=2n 이면 분할(decimation)에 의해 길이가 N/2 인 두신호는 길이가 N/4 인 두신호로 쪼개어진다. 이런 방식으로 계속 분할을 진행하면 결국 길이가 2인(N=2) 신호로 나누어진다. N=2 인 DFT는 매우 쉽게 계산된다.
그림(2)와 같이 길이가 8인(N=8)인 신호를 생각하자. 처음 길이가 8인 신호는 다니엘슨-란초스 보조정리에 의해 분할을 계속하면 stage 3과 같은 순서가 된다.
그런데 위의 표에서 처음 원래의 신호와 최종적으로 나누어진 신호의 순서를 이진수로 나타내면 그림에서와 같이 대응하는 신호의 모든 비트들을 거꾸로 배열한 것과 같다. 이것을 Bit reversal 이라 부른다. 이런 현상은 신호가 N=2n 라면 항상 나타나는 현상이다. 계산구조가 in-place(입력 신호 메모리 장소에 출력 값을 덮어쓰는 것. 다른 표현으로 원래의 샘플 메모리 안에서 FFT를 수행하고 추가적인 버퍼 메모리를 사용하지않는 것)이어야 한다는 요구조건을 포기한다면 입력자료와 출력 결과 둘 다 순서대로 나열되도록 하는 것이 가능하며, 이를 out-of-place라 한다.
쿨리-튜키 알고리즘(Cooley-Tukey algorithm)
가장 일반적으로 사용되는 FFT 알고리즘은 쿨리-튜키 알고리즘이다. 이 방법은 1965년에 J. W. 쿨리와 J. W. 튜키가 발표하여 널리 알려졌다.[5]
식(2)의 경우 k =0 에서 N-1까지의 각 항의 계산에 N 회의 곱셈이 필요하기 때문에, 전체적으로 N 2 회의 곱셈이 필요하게 된다. 그러나, 식(4)의 길이 N/2 의 DFT는 솔직하게 계산하면 (N/2) 2=N 2/4 회의 곱셈으로 실행할 수 있으므로 이 분할(Decimation)에 의해 계산량은 약 절반으로 줄어든다. 또한, 이 분할을 2회, 3회,... 반복하면 즉, 재사용(reuse)하면 계산량은 약 1/4, 1/8,... 로 격감한다. 이것이 쿨리-튜키 FFT (정확하게는 기수 2 쿨리-튜키 FFT)의 기본적인 생각이다.
쿨리-튜키 알고리즘은 보통 크기 N을 재귀적으로 2등분하여 길이가 2인 신호가 얻어질 때까지 분할 정복 알고리즘(Divide and conquer algorithm)을 적용하기 때문에 N = 2n인 경우에 많이 적용된다.
보통의 프로그램 코드는 log N 개의 단계(그림(2)나 그림(3)의 stage에 해당)에 대한 재귀 호출을 수행하고, 각 단계에서 N / 2 회의 복소수 곱셈과 N 회의 복소수 가산을 수행한다. 따라서 복소수 곱셈 횟수는 (N / 2) * log N 로 감소하고 부동 소수점 연산의 양은 N * log N의 order가 될 수 있다. 이것은 쿨리-튜키 FFT의 일반적인 연산량이다.
일반적으로 N = n1n2가 성립하는 n1과 n2는 같을 필요가 없으며, 따라서 N이 임의의 합성수일 때에도 적용 가능하다(아래 혼합 기수 FFT 를 참고). 쿨리-튜키 알고리즘은 DFT를 더 작은 크기의 DFT로 나누기 때문에, 뒤에 설명된 다른 FFT 알고리즘과 함께 사용될 수 있다.
푸리에 행렬(Fourier Matrix)
회전자로 이루어진 차수(order) N의 복소 방데르몽드 행렬((Vandermonde matrix)
을 푸리에 행렬(Fourier Matrix)이라고 한다.
식(2)를 아래와 같이 푸리에 행렬로 표현할 수 있다. 편의를 위해 지금부터는 회전자의 아래첨자 N을 생략하겠다. 예를들면, 이다.
식(3a)가 보여주는 회전자의 주기성으로 인해 부터 까지 모두 계산할 필요가 없이
까지만 계산하면 된다. 다음 기수 2 시분할 FFT의 식(6)에서 확인해 볼 수 있다.
기수 2 알고리즘(Radix-2 FFT Algorithm)
앞에 열거한 내용들을 기초하여 기수 2 FFT 알고리즘을 다루어본다. 이는 다른 알고리즘들을 다루는데 있어 가장 기초가 되는 알고리즘이다. 우선 기수의 뜻을 살펴보자.
기수(radix)의 의미
- 기수(基數)는 밑 또는 기저(radix, base)로도 부르고, 숫자 표현(진법 체계)의 기준이 되는 수다.
radix 는 라틴어로 뿌리(root)라는 뜻이다.
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 행렬식은 다음과 같다.
지수의 짝 홀을 기준으로 위의 식을 다음과 같이 열의 순서를 바꾸는 변형할 수 있다(이는 곧, 앞서 설명한 분할(decimation)을 의미한다).
또한 식(3a)로 부터 의 관계를 얻으므로
식(7) 오른쪽 행렬의 는 N=2 DFT 행렬 안에 있으므로 식(3c)를 이용하여 로 변경해서 고쳐쓰면
식(3b)를 이용하면 , , 이므로 다음과 같은 최종식을 얻는다.
따라서, N=4 인 DFT는 위 식(8)와 같이 N/2 즉,N=2 인 DFT 두 개를 사용해서 계산할 수 있고, 이는 다니엘슨-란초스 보조정리가 적용됨을 행렬을 이용해 보인 것이다. 그 결과로 입력과 출력이 서로 역비트순이 됨도 알 수 있다. 식(8)의 결과는 아래 그림(3)의 신호흐름도(signal flow graph)와 동일하며,나비도표(butterfly diagram)로도 불린다. 두 개의 N=2 DFT 나비를 볼 수 있다.
• 대체적으로 기수 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로 구현하는 경우의 상충관계(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=N1N2 ...Nm 일 때 사용한다.
N=2 n(즉, 기수 2)인 경우에 N=2, 4, 8,...이고, N=4 n 인 경우 N=4, 16, 64... 등의 갯수를 가질 수 있으나,
중간 갯수인 N=12 , N=15, N=24 등은 단일한 기수로는 담당할 수 없게 된다. 따라서 서로 다른 기수가 혼합된 여러 기수를 사용해 해결하는 것이다.
N=2 n×q s (q=2 m>2, 1<s≤m ) 이라면 기수 2 알고리즘 n 단계를 변환의 시작 또는 끝에 적용하고 나머지는 기수 q 알고리즘 s 단계를 수행한다.
예) N=2 (2m+1)=2 1×4 m
N=32=2 1×4 2 즉 기수2 한 단계와 기수4 두 단계를 수행한다.
2) 나비도표에 서로 다른 기수의 나비들이 존재한다.
이론
식(5)에서 DFT는 푸리에 행렬 T와 벡터 x, f 로서 표현됨을 알았다. 즉
T의 행렬요소는 이다.
신호수가 합성수 N이며, 다음과 같은 m개의 인수로 분해(decomposition)될 때
이고, N X N 행렬 T 는 다음과 같이 행렬곱으로 표현된다.
여기서, P는 치환행렬(permutation matrix), Fi는 n의 계수 ni에 해당하는 변환 단계(transform step) 행렬이다.
행렬 Fi는 각 행과 열에 ni개의 0이 아닌 요소만 있고 차원 ni의 n/ni 개의 square 부분행렬(submatrix)들로 분할될 수 있다. FFT 알고리즘의 기초가 되는 것이 바로 이 분할과 그에 따른 곱셈의 감소이다. 행렬 Fi는 다음과 같이 더 분해될 수 있다.
여기서 Ri는 회전자의 대각행렬이고 Ti는 ni 개의 동일한 square 부분행렬로 분할될 수 있으며 각 행렬은 차원 ni의 복소 푸리에 변환이다.
신호의 길이가 합성수 N=N1N2 이라고 하자. 두 수의 곱으로 표현되므로 m=2 인 경우에 해당한다. 이 신호의 DFT는 식(2)로부터
신호 x 를 N1 x N2 이차원 배열로 재배치하는 것을 생각해 볼 수 있다.
N1 =3, N2 =8 인 신호를 고려하자.
으로 간주한다. 를 각각 위 (17)식 h의 행 번호, 열 번호, 를 각각 아래 (18)식 H의 행 번호, 열 번호이라 할 때,
이라 하자. 그러면, 위의 이차원 배열을
와 같이 쓸 수 있다. 배열 h 에 대응하는 DFT는
을 생각할 수 있다. 위의 DFT를 이차원 배열로 나타내면
두 이차원 배열을 구성한 순서가 (17)식은 행우선순서(row-major order), 출력에 해당하는 (18)식은 열우선순서(column-major order)로 다름을 알 수 있고 이를 전치(transpose)되었다고 한다.
이러한 모든 전치의 최종 결과는 입력(DIF) 또는 출력(DIT) 인덱스의 비트 반전에 해당한다. 그 이유는 잠시 뒤 명확해 질 것이다.
그 전에 위에 언급한 의 범위로부터 (16)식과 같이 의 범위가 이 되는 지 확인해 보자.
<증명>
이라할 때
의 최소값은 일 때 이며,
의 최대값은 일 때 이므로
이다.
따라서 이다. <증명 끝>
이제 원래의 DFT (16)식을 다음과 같이 쓸 수 있다.
윗 식 세 개의 회전자 곱에서 처음 두 회전자들 때문에 이차원 DFT를 연상하게 한다. 마지막 회전자는 이들 전체 이차원 DFT에 대한 일종의 회전자이다. 이 사실은 조금 전에 두 이차원 배열을 구성한 순서를 다르게 한 즉, 전치(transpose)를 만든 이유이다.
윗 식을 다시 쓰면
괄호 안은 배열의 신호수 개의 각 열(column)에 대한 DFT이며, 이를 라 하자. 그러면,
이는 먼저 회전자를 곱한 후에 계산한 신호수 개의 행들에 대한 DFT 식이다. 이것과 이차원 DFT와의 유일한 차이점은 회전자에 의한 추가적인 곱셈이 있다는 것이다.
만약 합성수였던 원래 신호의 길이 N=N1N2에서 N1 이나 N2가 다시 합성수라면, 분해(decomposition)의 기본 아이디어를 재귀적으로 적용한다.
↑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.
↑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. JSTOR2003354.