07-6. Java로 FFT 알고리즘 충실히 구현하기

이번 페이지에서는 FFT 알고리즘을 Java 언어의 특성을 충실히 살리면서, 읽을 수 있는 수준의 알고리즘 구현을 해보고자 한다.

 

Java 언어의 특성을 살린다는 의미는, 복소수의 표현을 Complex라는 객체를 만들어서 사용함으로써 표현의 간단함과 확장성을 늘리고, 에러 처리를 Exception을 사용해서 하고, 메서드 오버로딩을 통해서 fft의 다양한 호출이 가능하게 한다는 의미.

 

읽을 수 있는 수준의 알고리즘이라는 것은, FFT의 구현을 각 의미있는 부분에 대해 별도의 메서드를 만들어서 구분되게 하고, 가능한 알고리즘에서 표현된 방법으로 코딩을 했다는 의미이다.

 

그러나, 위와 같이만 하면 속도에 손해를 보기에, 수행 속도에 지대한 영향을 끼치는 부분은 코드 가독성이 좀 떨어지더라도, 효율적인 코드로 대체했다.

 

---

1. 소스 코드 (테스트 코드 포함)

전체 소스 코드는 여기 참조.

 

2. 사용법

입력 데이터를 double형 배열로 만들고, FFT.fft를 호출하고, 그 결괏값을 Complex형 배열로 받으면 됨

 

1) 입력데이터 및 결괏값용 배열 생성

FFT의 입력은 double형 배열이다. 복소수가 아닌 double형이라는 것에 유의 (FFT는 일반 실수형 신호값에 대해 주파수 공간의 복소수형 값으로 변환한다.)

 

FFT의 반환값은 Complex형 배열이고, fft 메서드 내에서 해당 배열이 만들어져 리턴되기에, Complex[] X = null 로 한 후, 이 X에 반환값을 받으면 된다.

double[] x = {0,0,0,0,1,1,1,1};
Complex[] X = null;

2) FFT.fft 메서드 호출

fft는 static형 메서드이기에 FFT 클래스에 대한 인스턴스 생성 없이 그대로 호출하면 된다. 

 

X = FFT.fft(x);

• 입력값인 x는 fft 메서드 내에서 변경되지 않는다. 즉, x에 들어 있는 입력값이 그대로 유지된다. 
• 출력값인 X는 fft 메서드내에서 생성되어 리턴된다. 따라서 호출 전에 인스턴스를 만들 필요 없이 null로 지정한 후 fft 호출한 결과를 그대로 받으면 된다.

3) 결괏값 확인

FFT.fft의 결과는 Complex형 배열이고, Complex는 여기서 자체 정의한 클래스이다. Complex 클래스는 실수부 값인 re와 허수부 값인 im을 double형으로 가지고 있는 단순 클래스.

 

결괏값인 복소수를 "a+bi"형의 String으로 바꾸려면 Complex.toString 메서드를 이용한다. 

for(int i=0;i<X.length;++i) 
	System.out.print(X[i].toString()+", ");

• 만약 "a+bj"형태로 출력하려면 X[i].toString("j")로 하면 된다.

---

3. 테스트 코드 수행 결과

테스트용 샘플 코드는 Test_FFT.java에 있고, Test_FFT의 main 함수의 수행 결과는 다음과 같다.

 

** test_simple_fft **
4.0, -1.0+2.414213562373095i, 0, -1.0+0.4142135623730949i, 0, -0.9999999999999999-0.4142135623730949i, 0, -0.9999999999999997-2.414213562373095i, 
4.0, -1.0+2.41j, 0, -1.0+0.41j, 0, -1.0-0.41j, 0, -1.0-2.41j, 

** test_not_power_of_2 (auto zero padding) **
2.0, -1.707+0.707i, 1.0-1.0i, -0.293+0.707i, 0, -0.293-0.707i, 1.0+1.0i, -1.707-0.707i, 

** test_measure_time **
Execution Time(32kB): 21ms
[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): 22ms
[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, 

** test_use_file **
fft calculation time: 2ms
fft result is written to the file

 

• 위 테스트를 수행할 때, 테스트 파일을 읽어서 입력값으로 하는 테스트가 있기에 그 입력 파일이 존재해야 수행된다. 입력 파일은 input.txt 파일이고 프로그램 수행되는 폴더 내 test 폴더 안에 있어야 한다. 

---

4. 소스코드 설명

1) 사용된 클래스

3 개의 클래스를 만들어 사용한다.

 

• FFT: fft와 ifft 메서드를 가지고 있는 핵심 클래스
• Complex: 복소수형 데이터를 처리하기 위한 클래스. 실수부 re와 허수부 im을 가지는 단순 형태
• CMath: 복소수끼리의 연산이 구현된 클래스.

• 복소수 연산을 Complex 클래스의 메서드로 처리할 수도 있으나, 별도의 CMath를 두고 그 내부의 static 메서드로 구현했다. Complex 클래스를 단순화하기 위함이고, static 메서드로 처리함으로써 조금이라도 속도 향상을 꾀하고자 함이다.

2) fft와 ifft

FFT는 fft 메서드로, IFFT는 ifft 메서드로 수행 가능하다.

 

fft 메서드의 입력은 double형 배열이고, ifft의 입력은 Complex형 배열

 

fft 및 ifft를 수행할 때 round_digit를 지정하면, 지정된 자릿수에서 반올림한 결괏값을 얻을 수 있다. 이를 지정하지 않으면 java에서의 double형에서 지원하는 최대 자릿수까지 표시된다.

 

public static Complex[] fft(double[] src) throws Exception{
    return FFT.fft(src, -1);
}

public static Complex[] fft(double[] src, int roundDigit) throws Exception{
    Complex[] complexSrc = new Complex[src.length];

    for(int i=0;i<src.length;++i) complexSrc[i] = new Complex(src[i],0);

    Complex[] X = fft_forward(complexSrc,  false);

    // round up at roundDigit
    if(roundDigit >= 0) {
        double d = Math.pow(10.0, roundDigit);
        for(int i=0;i<X.length;++i) {
            X[i] = new Complex(Math.round(X[i].re * d)/d, Math.round(X[i].im * d)/d);
        }
    }		
    return X;
}

public static double[] ifft(Complex[] src) throws Exception{
    return FFT.ifft(src,-1);
}

public static double[] ifft(Complex[] src, int roundDigit) throws Exception{
    Complex[] complexArr = fft_forward(src,  true);

    double[] x = new double[complexArr.length];

    if(roundDigit >= 0) {
        double d = Math.pow(10.0, roundDigit);  // round up at roundDigit
        for(int i=0;i<x.length;++i) {
            x[i] = Math.round(complexArr[i].re * d) / d ;
        } 
    }else {
        for(int i=0;i<x.length;++i) {
            x[i] = complexArr[i].re;
        } 
    }		

    return x;
}

 

3) FFT 계산의 핵심인 fft_forward 메서드

FFT 알고리즘의 핵심이 구현된 메서드로, fft와 ifft에서 공용으로 사용될 수 있도록 작성했다. 

ifft의 경우는 fft_forward의 isInverse=False로해서 호출한다.

 

private static Complex[] fft_forward(Complex[] src, boolean isInverse) throws Exception {
    int N = src.length;

    //1. check src data size. If it's not 2^n size then make it to be 2^n size with zero padding
    int fillCount = 0;
    fillCount = checkPowerOfTwo(N);
    if(fillCount > 0) {
        src = makePowerOfTwoSize(src, fillCount);
    }else if(fillCount < 0) {
        throw new Exception("Something wrong: while calculateing filling count from the input data");
    }

    //2. calculate loop count: 2-> 1, 4->2, 8->3, 16->4		
    N = src.length; //input array can be changed. so, read one more time.
    int totalLoop = (int)FFT.log2(N); //it is guaranteed that N is 2^n number, therefore log2(N) will be int value. 

    //3. arrange src index to calculate FFT as bottom-up style
    Complex[] arrangedSrc = arrange4fft(src, totalLoop);

    //4. calculate W: 0 to (N/2-1)
    Complex[] W = calculateW(N, isInverse); //for fft. for ifft:calculateW(N, true)

    //5. use 2-dimensional array to save memory space. X[0, ] <-> X[1, ]
    Complex[][] X = new Complex[2][N];
    System.arraycopy(arrangedSrc, 0, X[0], 0, N);

    //6. calculate FFT by bottom-up style
    int srcIdx=0, tgtIdx=0;
    int kJump=0, regionSize =0;
    for(int curLoop=0; curLoop<totalLoop;++curLoop) {
        tgtIdx = (tgtIdx+1) % 2;
        srcIdx = (tgtIdx+1) % 2;
        regionSize = 1 << (curLoop+1);			//if N=8: 2 -> 4 -> 8
        kJump = 1 << (totalLoop - curLoop - 1);	//if N=8: 4 -> 2 -> 1
        int i=0;
        while(i<N) {
            regionFFT(X, srcIdx, tgtIdx, i, regionSize, kJump, W);
            i += regionSize;
        }
    }

    //7. find final output index in the X array
    int resultIdx = ((totalLoop & 1) == 0) ? 0 : 1;

    //8. in case of IFFT, divide N
    if(isInverse) {
        for(int i=0;i<N;++i) X[resultIdx][i] = CMath.divideR(X[resultIdx][i], N);			
    }

    //9. return the result as 1d array
    return X[resultIdx];
}

 

-끝-

 

 

이전글: [07-5] Fastest FFT code for Java

다음글: [07-6 Code] FFT program for Java

다음다음글: 07-7. C로 짠 FFT Code