고속 푸리에 변환

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

신호 이 복소수라고 가정할 때, DFT는 다음과 같이 정의한다.

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

 을 회전자(Twiddle factor)라고 부른다. 1의 거듭제곱근(Root of unity), 단위근도 회전자의 다른 호칭이다. 이 회전자를 사용해 식(1)을 다시 정리하면,
 
회전자는 N에 관한 주기성으로 인해 다음과 같은 성질이 있음을 쉽게 알수있다.

n:자연수
(예 1)

쿨리-튜키 알고리즘

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

→

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

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

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

길이가 짝수이므로 각 항의 갯수는 같다. 두번째 괄호 항에서 을 밖으로 빼면

이를 다시 쓰면

그런데,

가 성립하므로,

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

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

 from (3a),(3b)

fk+N/2 의 실제 범위는 f<sub>k<sub> 와 달리 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. 튜키가 발표하여 널리 알려졌지만, 나중에 카를 프리드리히 가우스가 1805년에 같은 방법을 이미 개발하였으며, 그 뒤에도 제한된 형태의 FFT가 종종 발견되었음이 밝혀졌다.

식(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 = n1 n2가 성립하는 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)라는 뜻이다.

(예) 10진수의 기수(base)는 10이다.
10 1, 10 2 , 10 3 의 밑(base)은 10이다.

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

(예) N=2 n 는 기수 2(radix-2)이다. N=4 n 는 기수 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 행렬식은 다음과 같다.
 
지수의 짝 홀을 기준으로 위의 식을 다음과 같이 열의 순서를 바꾸는 변형할 수 있다(이는 곧, 앞서 설명한 분할(decimation)을 의미한다).
또한 식(3a)로부터  다음의 관계를 얻으므로

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

윗 식 오른쪽 행렬의  는 N=2 DFT 행렬 안에 있으므로 식(3c)를 이용하여  로 변경해서 고쳐쓰면

식(3b)를 이용하면  ,  ,  이므로 다음과 같은 최종식을 얻는다.

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

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

여기서, I2: 2x2 단위행렬, T2: 회전자로 구성된 2x2 대각행렬, O2: 2x2 Null 행렬  C4: 4x4 column 치환행렬(permutation matrix)

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

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

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

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

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

 → 

 → 

 → 

 → 

…             …

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

 → 
 → 
…         …
 → 
――――――――――――――――――――――――
 → 

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

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

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

→

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

윗 식 둘째 항은 DFT가 아니다.  그런데,  이므로,

이제  에 대하여

이라고 하자. 그러면 짝수  에 대하여

이고, 홀수  에 대하여

여기서  를 이용하였다.  이제 위의 두 식은 각각 길이가 N/2인 과 의 이산 푸리에 변환(DFT)이 되었다.

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

이를 길이가 2인 신호들이 얻어질 때까지 반복한다(분할정복 알고리즘). 기본적인 주파수 분할 나비는 식(13)으로부터 다음과 같이 구성된다. 
다음 식(15)은 N=4일때 기수 2 주파수분할(DIF) FFT이다.  식(6)의 시분할(DIT) 방식과 달리  출력 f 의 원소들이 분할된 것을 볼 수 있다.

다음은 (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  인 경우이다. 그러나 N 은 다양한 기수를 가질 수 있다. 즉, 일반적으로 기수 r 일 때, N=r n이다. 그림(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는

(예)

N=16 이고 기수 r=2 이면 stage 수 이다.

N=16 이고 기수 r=4 이면 stage 수 이다.

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

4)  기수에 따른 계산 효율

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

표의 결과만을 볼 때 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=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)의 기본 아이디어를 재귀적으로 적용한다.

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

1. N2 크기의 N1 DFT를 수행한다.

2. 회전자를 곱한다.

3. N1 크기의 N2 DFT를 수행한다.

일반적으로 N1 또는 N2는 기수이다. N1이 기수이면 DIT(Decimation in Time) 알고리즘이라고 하고, N2가 기수이면 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

푸리에 해석학
디지털 신호 처리
이산 변환