이 페이지는 Java로 구현하는 FFT, 그것도 수행 속도를 빠르게 하는 것에 목적을 둔 프로그래밍 코드를 소개한다.
전체 소스코드는 여기 참조
앞 페이지에서 엑셀로 만드는 FFT의 경우는, 속도는 느리지만 엑셀 자체에서 데이터를 다루거나 분석을 하고, 그래프를 그리고 등 부가적인 기능으로 해서, 필요성이 있어 FFT 프로그램을 만들어서 소개했다.
Java로는 만들 생각을 안했었다.
FFT가 주로 쓰이는 것이, 연구 목적으로 신호를 분석하거나(그래서 엑셀에서의 FFT가 유용), 실시간으로 신호를 FFT 변환을 해야 하는 것이기에, C나 C++을 이용해서 코드 수행을 빠르게 하는 게 중요하지, Java 같은 Virtual Machine 기반의 언어로는 유용성이 떨어질 것이기 때문이다.
그런데, 안드로이드에서라던가, 서버단에서, Java로 된 FFT 코드를 필요로 하는 경우도 있겠다는 생각이 들었다. 그리고선, Java로 된 FFT 코드들을 검색해봤는데, 역시 C나 C++에 비해 완성도가 좀 떨어지는 느낌이다.
C로 된 코드를 그대로 Java로 포팅한 느낌이 나는 코드이거나, Java의 특징을 살리지 못한 코드라던가, 딱 마음에 드는 코드가 없다.
해서, 두 종류의 Java로 된 FFT 코드를 작성했다.
하나는, 속도를 중시한, Java 코드로 낼 수 있는 최대한의 퍼포먼스를 내도록 짠 코드이고(이 코드는 이 페이지에서 설명), 다른 하나는 속도는 조금 떨어지지만 좀 Java 스럽게 객체를 사용하고 확장성이나 다른 코드에서 사용하기 편리하게 작성했다. (이 코드는 다다음 페이지에서 소개)
수행 속도
작성된 FFT 코드로 32 kBytes와 64 kBytes 크기의 데이터에 대해 FFT를 구해봤을 때, 5ms와 7ms가 소요되었다.
| Input Data Size | FFFT (이 페이지에서 설명할 FFT 코드) |
| 32 kBytes | 5ms |
| 64 kBytes | 7ms |
* 사용된 컴퓨터 사양: 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz
코드 작성 기본 방침
- 수행 속도를 최우선으로 한다.
- 그러나, 코드가 너무 지저분해지도록 옵티마이제이션을 하진 않는다.
- 클래스의 명칭은 Fastest FFT의 의미로, FFFT로 한다. (Fastest Fast Fourier Transform)
속도 향상을 위한 알고리즘 기법
1. $O(N\log{N})$인 쿨리-튜키 알고리즘을 Bottom-up 방식 사용
당연한 얘기겠지만, 재귀호출 사용 안 하고 Bottom-up 방식으로 구현
2. 콘볼루션 링인 $W^k$값을 한 번만 계산한 후 재활용
이것도 대부분의 구현에서 이렇게 하고 있는 것. 그러나, 알고리즘만 보고 처음 구현할 때는, 매 라운드마다 달라지는 N에 대해 $W_N^k$를 계산하는 방식을 취할 수도 있다.
여기서는 당연하게 $W^k$를 한 번 구한 후 사용하고, $W^0$ ~ $W^{\frac{N}{2}-1}$까지만 구해서 사용.
3. Bottom-up에서 가장 신경 쓰이는, 입력 데이터에 대한 재정렬 부분
N개의 입력 데이터에 대해서 순서 조정을 위해서 $N\log{N}$ 정도의 연산이 필요한데, 이 부분을 가능한 최적화를 한다. 최적화된 코드는 아래와 같다.
//3. calculate index for bottom-up style fft and arrange input data
int nBits = (int)(Math.log(N) / Math.log(2));
int k=0;
double tmp;
for(int num=1;num<(N-1);++num) { // 0 and last num is not changed. so [1,N-2]
int d = 0;
for(int i=0;i<nBits;++i) {
k = (int) (num >> (nBits-1-i));
d = d + ((k & 1) << i);
}
if(num < d) { //change each other
tmp = XR[num];
XR[num] = XR[d];
XR[d] = tmp;
}
}
- 기본 원리는 십진수 $\rightarrow$ 이진수 변환 -> 이진수에 대해 Reverse -> 십진수로 변환.
- 위와 같은 변환 과정에서, 수의 절반을 기준으로 해서 서로 교환되는 형태이기에, 절반 정도만(if num<d) 실제 교환이 이뤄지도록 했다.
- 그리고, 입력 데이터로 사용되는 XR 배열 하나를 가지고 교환이 이뤄지게 해서, 새로운 배열을 만들지 않아도 되게 했다.
속도 향상을 위한 프로그래밍 기법
1. static 메서드 하나로 작성
- 여러 메서드를 호출하면, 메서드 전환에 오버헤드 발생한다. 해서 하나의 메서드로.
- static 메서드로 하면, 클래스 인스턴스의 생성 없이 바로 호출 가능해서 편하고, 속도도 빠르다.
2. 객체 생성을 최소화하기 위해, 배열만 사용
FFT의 계산 결과는 복소수 형태이기에, 코드 가독성을 위해서는 Complex용 클래스를 만들어서 사용하는 것이 좋다. 그러나, 이렇게 하면 여러 개의 복소수의 표현을 위해서 계속해서 인스턴스를 만들어야 하기에 속도에는 도움이 안 된다.
해서, 여기서는 복소수의 표현을 실수부용 배열과 허수부용 배열을 만들어서 사용한다.
복소수 하나의 표현에 2개의 수가 필요하기에 코드가 복잡해지긴 하지만, 속도 향상에는 도움이 된다.
3. 배열의 생성을 최소화하기 위해, fft 호출할 때 파라미터로 사용한 배열을 출력값으로 사용한다.
일반적인 자바 프로그래밍 습관이라면 메서드의 리턴값으로 결과를 받게 코딩을 할 것이다.
double[] fft_result = FFFT.fft(double[] x_real, double[] x_imagine);위와 같이 파라미터로 입력 배열을 넘기고, 계산된 결과를 배열로 받는, 이런 형태가 일반적이다.
그런데, 위 코드는 아래와 같이 변경 가능하고, 아래 경우가 오버헤드가 적다.
double[] XR = ...
double[] XI = ...
FFFT.fft(XR, XI);이 경우는, fft를 호출하기 전에 XR에는 입력 데이터를, XI에는 0으로 초기화된 값을 넣어서 준비해뒀다가 fft 메서드의 파라미터로 배열을 넘기게 되고, fft내에서는 이 배열을 재활용하여 결괏값을 넣어서 전달한다.
Java에서 메서드를 호출할 때 파라미터는 모두 Call by Value이나, 배열의 경우에는 Call by Reference처럼 동작하기에, 위와 같은 트릭이 가능하다.
이렇게 함으로써 메모리의 사용을 최소화하게 되고, 이에 따라 속도 향상에 도움이 된다.
4. 나눗셈의 사용을 최대한 줄인다.
2의 배수승으로 나누는 경우, 나눗셈을 사용하지 않고 Shift 연산을 사용한다.
regionSize를 2로 나누는 경우,
half = regionSize >> 1;
5. FFT의 계산이 이뤄지는 메인 루틴 최적화
코드 가독성을 위해서는 블록 사이즈별 계산되는 부분을 별도의 메서드로 만들면 좋으나, 그냥 하나의 메서드 내에서 중첩 for 구문을 사용해서 계산이 이뤄지게 했다.
//5. calculate fft
double TR, TI; //T = X[j+half] * W[k]
int regionSize, kJump, half;
int blockEnd;
for(int loop=0; loop < nBits; ++loop) { //nBits == totalLoop
regionSize = 1 << (loop+1); //if N=8: 2 -> 4 -> 8
kJump = 1 << (nBits - loop -1); //if N=8: 4 -> 2 -> 1
half = regionSize >> 1;
for(int i=0; i<N; i += regionSize) {
blockEnd = i + half - 1;
k = 0;
for(int j=i; j<=blockEnd; ++j) {
//T = X[j+half] * W[k]
TR = WR[k] * XR[j+half] - WI[k] * XI[j+half];
TI = WI[k] * XR[j+half] + WR[k] * XI[j+half];
XR[j+half] = XR[j] - TR;
XI[j+half] = XI[j] - TI;
XR[j] = XR[j] + TR;
XI[j] = XI[j] + TI;
k += kJump;
}
}
}여기서도 새로운 배열을 만들어서 사용하지 않고, 메서드 호출할 때 파라미터로 넘어온 XR과 XI를 그대로 활용하도록 했다. 그 이유로, XR[j+half]를 XR[j]보다 먼저 계산되도록 했다.
수행 방법
입력으로 넣을 데이터를 배열로 만들고, FFFT.fft 메서드를 호출하면 끝.
1) 입력 및 출력으로 사용될 double형 배열 2개를 만든다. 실수부용과 허수부용
N = (int)Math.pow(2, 15);
double[] XR = new double[N];
double[] XI = new double[N];- 배열의 크기는 2의 지수승이라야 한다.
- 만약 데이터가 2의 지수승개가 아니라면, 그 데이터의 개수보다 큰 바로 위의 2의 지수승개의 배열을 만들고, 모자라는 데이터 부분은 0으로 채워 넣으면 된다. 예을 들어, 데이터의 개수가 14개라면, N=16으로해서 배열을 만들고, 모자라는 2개 부분은 0으로 넣으면 됨
2) 실수부용 배열에 입력 데이터를 넣는다.
허수부용 배열은 그대로 놔두면 되고, 그러면 초기값인 0이 전달된다. (FFT에서 입력 신호는 실수부 값만 있다. 그 결괏값이 실수부와 허수부의 복소수 형태)
for(int i=0;i<N;++i) XR[i] = Math.cos(2*Math.PI*0.004*i);- 위 코드는 정현파 신호값을 생성해서 넣는 코드.
- 실제 신호를 파일에서 읽어서 넣어도 되고, 다른 배열에 있는 값을 카피해 넣어도 된다.
- 신호값이 2의 지수승개가 아니면, 모자라는 부분을 0으로 넣어서 입력하면된다.
3) FFFT.fft 메서드를 호출
FFFT.fft(XR, XI);- XR에는 입력 신호값이 들어가 있어야 한다. 허수부용 배열인 XI의 값은 아무값이나 있어도 된다. 즉, 허수부용 배열은 FFT 계산할 때 참조하지 않는다.
- 수행이 완료되면, 결괏값이 XR과 XI에 넣어지게 된다. 즉, 입력으로 사용했단 XR의 데이터가 전부 날아가고 새로운 FFT 결괏값으로 채워지는 것에 유의
- XR에는 FFT결과의 실수부 값이, XI에는 허수부 값이 들어가게 됨
4) 결괏값의 출력
XR과 XI 배열을 그대로 사용하면 되는데, "2+3i" "-3i" "3"과 같이 보기 좋게 String으로 출력하려면, FFFT.getCplxStr 메서드를 사용하면 된다.
for(int i=0;i<N;++i)
System.out.println("["+i+"]"+FFFT.getCplxStr(XR[i],XI[i])+", ");테스트 코드
Test_FFFT.java에 테스트용 샘플 코드가 들어 있다.
//Test_FFFT.java
package test;
import fourier.FFFT;
public class Test_FFFT {
public static void main(String[] args) {
try {
test_simple_fft();
test_measure_time();
} catch (Exception e) {
e.printStackTrace();
}
}
// 1. test for simple a few data
private static void test_simple_fft() throws Exception{
int N=16;
double[] XR = {0,0,0,0,1,1,1,1,0,0,0,0,1,1,1,1}; //input
double[] XI = new double[N]; // initial value is zero
System.out.println("** test_simple_fft **");
//1. calculate FFT for a few data. (count = 8)
FFFT.fft(XR, XI);
for(int i=0;i<XR.length;++i) System.out.print(FFFT.getCplxStr(XR[i], XI[i])+", ");
}
// 2. measure execution time to 32k, 64k size data
private static void test_measure_time() throws Exception{
long beforeT, diffT;
int N = 0;
System.out.println("\n** test_measure_time **");
//1) 32kB data
N = (int)Math.pow(2, 15);
double[] XR32 = new double[N];
double[] XI32 = new double[N];
for(int i=0;i<N;++i) XR32[i] = Math.cos(2*Math.PI*0.004*i);
beforeT = System.currentTimeMillis();
FFFT.fft(XR32, XI32);
diffT = System.currentTimeMillis() - beforeT;
System.out.println("Execution Time(32kB): " + diffT + "ms");
for(int i=0;i<4;++i) System.out.println("["+i+"]"+FFFT.getCplxStr(XR32[i],XI32[i])+", ");
System.out.println("......");
for(int i=N-4;i<N;++i) System.out.println("["+i+"]"+FFFT.getCplxStr(XR32[i],XI32[i])+", ");
System.out.println("");
//2) 64kB data
N = (int)Math.pow(2, 16);
double[] XR64 = new double[N];
double[] XI64 = new double[N];
for(int i=0;i<N;++i) XR64[i] = Math.cos(2*Math.PI*0.004*i);
beforeT = System.currentTimeMillis();
FFFT.fft(XR64, XI64);
diffT = System.currentTimeMillis() - beforeT;
System.out.println("Execution Time(64kB): " + diffT + "ms");
for(int i=0;i<4;++i) System.out.println("["+i+"]"+FFFT.getCplxStr(XR64[i],XI64[i])+", ");
System.out.println("......");
for(int i=N-4;i<N;++i) System.out.println("["+i+"]"+FFFT.getCplxStr(XR64[i],XI64[i])+", ");
}
}
테스트 코드 수행 결과
위의 Test_FFFT.java의 main 함수의 수행 결과
** test_simple_fft **
8.0, 0, -2.0+4.82842712474619i, 0, 0, 0, -2.0+0.8284271247461898i, 0, 0, 0, -1.9999999999999998-0.8284271247461898i, 0, 0, 0, -1.9999999999999993-4.82842712474619i, 0,
** test_measure_time **
Execution Time(32kB): 5ms
[0]17.441665612932592,
[1]17.442678036016932+0.03054032554147773i,
[2]17.445716012601977+0.06109131908951394i,
[3]17.450781665925863+0.09166366108017976i,
......
[32764]17.457878538842998-0.1222680568280476i,
[32765]17.450781665924485-0.09166366107969492i,
[32766]17.445716012601984-0.06109131908950698i,
[32767]17.442678036016957-0.03054032554145092i,
Execution Time(64kB): 7ms
[0]31.474846045104538,
[1]31.475301314548666+0.058005922412330935i,
[2]31.4766672023788+0.11601690943639723i,
[3]31.478943947152842+0.17403802715017125i,
......
[65532]31.48213194657478-0.2320743446023137i,
[65533]31.47894394715285-0.1740380271501539i,
[65534]31.476667202378863-0.11601690943633183i,
[65535]31.47530131454865-0.05800592241228639i,- CPU나 메모리 상황에 따라, 그리고 수행할 때마다의 캐시 상황 때문에, Execution Time의 차이가 있을 수 있다.
- 위 수행에 사용한 컴퓨터 사양은, 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz
-끝-
| 이전글: [07-4] Fastest FFT code for Excel VBA |
| 다음글: [07-5] Fastest FFT code for Java |
| 다음다음글: 07-6. Java로 FFT 알고리즘 충실히 구현하기 |
'푸리에 변환, 신호 > 푸리에 변환의 모든 것' 카테고리의 다른 글
| 07-6. Java로 FFT 알고리즘 충실히 구현하기 (0) | 2022.08.05 |
|---|---|
| [07-5 Code] Fastest FFT code for Java (0) | 2022.08.05 |
| [07-4 Code] Fastest FFT code for Excel VBA (0) | 2022.08.05 |
| 07-4. 엑셀에서 직접 FFT 프로그램 만들기 (7) | 2022.07.28 |
| 07-3. 엑셀에서 FFT 계산하기 (엑셀 자체 FFT 기능) (0) | 2022.07.28 |
