소인수 알고리즘

소인수 FFT 알고리즘(Prime-factor FFT algorithm)이라고도 칭하며, 발견자의 이름을 사용해 Good–Thomas algorithm(1958/1963)으로도 불리는 소인수 알고리즘(Prime-factor Algorithm, PFA)은 중국인의 나머지 정리(Chinese Reminder Theorem, CRT)를 활용한 고속 푸리에 변환(FFT) 알고리즘이다. 이 알고리즘은 크기 N = N₁N₂인 신호를 2차원 이산 푸리에 변환(DFT) N = N₁ X N₂ 로 재표현(re-expression)한다. 다만, 여기서 N₁과 N₂는 서로소(relatively prime), 즉 gcd(N₁, N₂)=1이어야 한다.

이 알고리즘의 장점은

1. 회전자(twiddle factor)에 의한 곱셈 연산(multiplication)이 불필요하다. 따라서, 처리 속도가 빠르다.
2. 변환의 in-place, in-order version을 구현하는 것이 가능하다. 즉, 변환된 값들은 메모리에서 입력값들을 바꾼다. 그리고 입력값들이 시간순서일때 출력은 주파수 순서이다.

단점은

1. 앞서 언급한 것처럼 N₁과 N₂는 '서로 소'이어야한다.
2. 중국인의 나머지 정리(CRT)를 사용하기때문에 더 복잡한 색인 재지정(re-indexing)을 수행해야한다.

소인수 알고리즘은 쿨리-튜키 알고리즘의 일반적 형태인 혼합기수(mixed-radix)알고리즘과 결합하여 사용할 수 있어서, 위에서 언급한 특별한 경우만이 아닌 임의의 크기 N 에도 적용할 수 있다. 데이터 크기 N=1000인 경우에 처리 예를 다음 그림에서 보여준다.

Good-Thomas 알고리즘과 혼합기수 Cooley-Tukey 알고리즘을 이용한 데이터 수 N=1000의 FFT 처리 방식.

알고리즘

신호의 길이가 N 일 때 이산 푸리에 변환(DFT)은

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

신호의 길이가 N = r s인 경우를 생각해 봅니다. r 과 s는 서로소인 두 정수이다(위에서 언급한 N₁, N₂에 해당). 입력과 출력을 재색인화(re-indexing)하고 이를 DFT 공식에 대입하여 크기가 r X s인 2차원 DFT로 변환해 보자. 이와같이 일차원 DFT를 다차원으로 변환하는 것을 다중차원 사상(multi-dimensional mapping) 이라고 한다.

색인 재지정(re-indexing)

정수론의 이론을 이용하여 (식 1) DFT의 입력 쪽 색인 n과 출력 쪽 색인 k를 다음과 같이 나타낼 수 있다.

1) 굿스 매핑(Good's mapping)

$n=(n_{1},n_{0})\equiv sn_{0}+rn_{1}({\bmod {N}})\qquad (0\leq n_{0}\leq r-1,\quad 0\leq n_{1}\leq s-1)\qquad \qquad \qquad (2)$

최대공약수 표현 정리에 의하면 각 n에 대해 위 식을 만족하는 정수 n₁과 n₂가 있고, 합동이론에 따라 위의 식을 만족하고 위의 범위에 속하는 정수 n₁과 n₂가 유일하다는 것을 쉽게 알 수 있다. 이를 루리타니안 매핑(Ruritanian mapping) 또는 굿스 매핑(Good's mapping)이라 한다.

2) CRT 매핑(CRT mapping)

$k\equiv k_{0}({\bmod {r}})\qquad (0\leq k_{0}\leq r-1)$

$k\equiv k_{1}({\bmod {s}})\qquad (0\leq k_{1}\leq s-1)\qquad \qquad \qquad (3)$

중국인의 나머지 정리(CRT)에 따르면 모든 쌍 (k₁, k₂)에 대해 (식 3)의 두 방정식을 만족하는 0 과 N-1 사이에 고유한 k가 존재한다. 이를 CRT 매핑이라고 하며, 이를 수식으로 표현하면 다음과 같다.

$k=(k_{1},k_{0})\equiv [ss^{-1}({\bmod {r}})k_{0}+rr^{-1}({\bmod {s}})k_{1}]{\bmod {N}}\qquad \qquad \qquad (4)$

n과 k를 위와 같이 (식 3)과 (식 4)로 나타낼 수 있는 이유는r 과 s 가 '서로 소'이기 때문이다. 위의 전개와는 반대로 입력 n을 CRT 매핑하고 출력 k를 Good's 매핑 할 수도 있다.

재표현(re-expression)

위의 두 식 (식 2)와 (식 4)를 (식 1)에 대입하면 식의 회전자(twiddle factor)는 다음 (식 5)와 같은 형태를 가지게 된다.

$W_{rs}^{[sn_{0}+rn_{1}]{\bmod {N}}[ss^{-1}({\bmod {r}})k_{0}+rr^{-1}({\bmod {s}})k_{1}]{\bmod {N}}}$

회전자 자체가 mod N 의 성질을 가지고 있으므로 지수(exponent)에 나타나는 mod N 은 생략할 수 있다. 즉,

$=W_{rs}^{[sn_{0}+rn_{1}][ss^{-1}({\bmod {r}})k_{0}+rr^{-1}({\bmod {s}})k_{1}]}$

$=W_{rs}^{sn_{0}ss^{-1}({\bmod {r}})k_{0}+sn_{0}rr^{-1}({\bmod {s}})k_{1}+rn_{1}ss^{-1}({\bmod {r}})k_{0}+rn_{1}rr^{-1}({\bmod {s}})k_{1}}$

어떤 임의의 수 m에 대해 $W_{rs}^{mrs}=1$ 이므로 윗 식 지수항의 두번째, 세번째 항은 생략 할 수 있다.

$=W_{rs}^{sn_{0}ss^{-1}({\bmod {r}})k_{0}+rn_{1}rr^{-1}({\bmod {s}})k_{1}}$

$=W_{r}^{n_{0}ss^{-1}({\bmod {r}})k_{0}}W_{s}^{n_{1}rr^{-1}({\bmod {s}})k_{1}}$

r과 s는 '서로 소'이므로 각각의 역수 $r^{-1},s^{-1}$ 이 존재하며 $ss^{-1}\equiv 1{\bmod {r}},\quad rr^{-1}\equiv 1{\bmod {s}}$ 이기 때문에, 윗식의 $W_{r}^{ss^{-1}{\bmod {r}}}=W_{r}^{1},\quad W_{s}^{rr^{-1}{\bmod {s}}}=W_{s}^{1}$이다. 따라서,

$=W_{r}^{n_{0}k_{0}}W_{s}^{n_{1}k_{1}}$ $\qquad \qquad (5)$

(식 2),(식 4), (식 5)로부터 (식 1)의 이산 푸리에 변환(DFT)은 최종적으로 다음과 같이 2차원으로 재표현된다.

$f(k_{1},k_{0})=\sum _{n_{1}=0}^{s-1}\sum _{n_{0}=0}^{r-1}x(n_{1},n_{0})W_{r}^{n_{0}k_{0}}W_{s}^{n_{1}k_{1}}\qquad \qquad (6)$

고속 푸리에 변환 위키 페이지에서 볼 수 있는 혼합 기수(mixed-radix) FFT의 다음(식 20)의 회전자 $W_{N_{1}N_{2}}^{k_{1}n_{2}}$가 위의 (식 6)에서는 생략된 것을 확인할 수 있다.

$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}}$

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

예제

r=3, s=4 인 신호의 소인수 FFT 를 적용해보자. 이 경우 (식 2)는

$n=(n_{1},n_{0})\equiv 4n_{0}+3n_{1}({\bmod {1}}2)\qquad (0\leq n_{0}\leq 2,\quad 0\leq n_{1}\leq 3)$

이며, (식 4)에서 $rr^{-1}({\bmod {s}})=33^{-1}({\bmod {4}})=9,\qquad ss^{-1}({\bmod {r}})=44^{-1}({\bmod {3}})=4$ 가 성립하므로 k 는

$k=(k_{1},k_{0})\equiv 4k_{0}+9k_{1}({\bmod {1}}2)\qquad (0\leq k_{0}\leq 2,\quad 0\leq k_{1}\leq 3)$

이 된다. (식 6)은

$f(k_{1},k_{0})=\sum _{n_{1}=0}^{3}\sum _{n_{0}=0}^{2}x(n_{1},n_{0})W_{3}^{n_{0}k_{0}}W_{4}^{n_{1}k_{1}}$

그럼 $(n_{1},n_{0})$이 변함에 따라 $x(n)=x(n_{1},n_{0})$이 어떻게 2차원 도표로 변하는 지, $(k_{1},k_{0})$에 따라 $f(k)=f(k_{1},k_{0})$ 는 어떻게 되는 지를 표로 만들어 보자.

1차원 신호 데이터 원본

n	0	1	2	3	4	5	6	7	8	9	10	11
xn	x0	x1	x2	x3	x4	x5	x6	x7	x8	x9	x10	x11

Good's mapping

x(n1,n0)		n0
		0	1	2
n1	0	x(0)	x(4)	x(8)
	1	x(3)	x(7)	x(11)
	2	x(6)	x(10)	x(2)
	3	x(9)	x(1)	x(5)

CRT mapping

f(k1,k0)		k0
		0	1	2
k1	0	f(0)	f(4)	f(8)
	1	f(9)	f(1)	f(5)
	2	f(6)	f(10)	f(2)
	3	f(3)	f(7)	f(11)