반응형
This page introduce very fast FFT code for Java.
It will execute FFT calculation for 64 kBytes input data in 10ms.
Execution result of Test_FFFT class
** 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,
Source Code
// FFFT.java
// 2022.8.5 Created by HJ
package fourier;
public class FFFT {
public static void fft(double[] XR, double[] XI) throws Exception{
//1. check 2^n count -> error. no zero padding
int N = XR.length;
if ( (N !=0) && ((N & (N-1)) != 0) )
throw new RuntimeException("The input size should be power of 2");
if ( N != XI.length )
throw new RuntimeException("The size of Real and Imaginary array should be same");
//2. calculate W
int N2 = N >> 1; //N2 = N/2
double[] WR = new double[N2]; //real part
double[] WI = new double[N2]; //imaginary part
double T = 2 * Math.PI / N;
double theta=0;
for(int i=0;i<N2;++i) {
theta = T*i;
WR[i] = Math.cos(theta);
WI[i] = -Math.sin(theta);
}
//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;
}
}
//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;
}
}
}
}
// make "a+bi" style complex string
public static String getCplxStr(double a, double b) {
String s = "";
String im_ch = "i";
if(a==0 && b==0) {
s = "0";
}else if(a == 0) {
s = b + im_ch;
}else if(b == 0) {
s = "" + a;
}else if(b > 0) {
s = a + "+" + b + im_ch;
}else {
s = "" + a + b + im_ch;
}
return s;
}
}
//Test_FFFT.java
//2022.8.5 Created by HJ
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])+", ");
}
}
- The End -
Before: 07-5. 가장 빠른 Java용 FFT 구현해보기 |
Next: 07-6. Java로 FFT 알고리즘 충실히 구현하기 |
Next: 07-7. C로 짠 FFT Code |
반응형
'푸리에 변환, 신호 > 푸리에 변환의 모든 것' 카테고리의 다른 글
[07-6 Code] FFT program for Java (0) | 2022.08.06 |
---|---|
07-6. Java로 FFT 알고리즘 충실히 구현하기 (0) | 2022.08.05 |
07-5. 가장 빠른 Java용 FFT 구현해보기 (2) | 2022.08.05 |
[07-4 Code] Fastest FFT code for Excel VBA (0) | 2022.08.05 |
07-4. 엑셀에서 직접 FFT 프로그램 만들기 (7) | 2022.07.28 |