要结合算法导论理解,参考:http://blog.csdn.NET/fjssharpsword/article/details/53281889
FFT的迭代实现,可以实现并行电路,和比较网络中的比较器思想一样。具体代码如下:
package sk.mlib;
public class FFTItrator {
int n, m;
// Lookup tables. Only need to recompute when size of FFT changes.
double[] cos;
double[] sin;
double[] window;
public FFTItrator(int n) {
this.n = n;
this.m = (int)(Math.log(n) / Math.log(2));
// Make sure n is a power of 2
if(n != (1<<m))
throw new RuntimeException("FFT length must be power of 2");
// precompute tables
cos = new double[n/2];
sin = new double[n/2];
for(int i=0; i<n/2; i++) {
cos[i] = Math.cos(-2*Math.PI*i/n);
sin[i] = Math.sin(-2*Math.PI*i/n);
}
makeWindow();
}
protected void makeWindow() {
// Make a blackman window:
// w(n)=0.42-0.5cos{(2*PI*n)/(N-1)}+0.08cos{(4*PI*n)/(N-1)};
window = new double[n];
for(int i = 0; i < window.length; i++)
window[i] = 0.42 - 0.5 * Math.cos(2*Math.PI*i/(n-1)) + 0.08 * Math.cos(4*Math.PI*i/(n-1));
}
public double[] getWindow() {
return window;
}
/***************************************************************
* fft.c
* Douglas L. Jones
* University of Illinois at Urbana-Champaign
* January 19, 1992
* http://cnx.rice.edu/content/m12016/latest/
*
* fft: in-place radix-2 DIT DFT of a complex input
*
* input:
* n: length of FFT: must be a power of two
* m: n = 2**m
* input/output
* x: double array of length n with real part of data
* y: double array of length n with imag part of data
*
* Permission to copy and use this program is granted
* as long as this header is included.
****************************************************************/
public void fft(double[] x, double[] y)
{
int i,j,k,n1,n2,a;
double c,s,e,t1,t2;
// Bit-reverse 二进制位翻转置换
j = 0;
n2 = n/2;
for (i=1; i < n - 1; i++) {
n1 = n2;
while ( j >= n1 ) {
j = j - n1;
n1 = n1/2;
}
j = j + n1;
if (i < j) {
t1 = x[i];
x[i] = x[j];
x[j] = t1;
t1 = y[i];
y[i] = y[j];
y[j] = t1;
}
}
// FFT
n1 = 0;
n2 = 1;
for (i=0; i < m; i++) {
n1 = n2;
n2 = n2 + n2;
a = 0;
for (j=0; j < n1; j++) {
c = cos[a];
s = sin[a];
a += 1 << (m-i-1);
for (k=j; k < n; k=k+n2) {
t1 = c*x[k+n1] - s*y[k+n1];
t2 = s*x[k+n1] + c*y[k+n1];
x[k+n1] = x[k] - t1;
y[k+n1] = y[k] - t2;
x[k] = x[k] + t1;
y[k] = y[k] + t2;
}
}
}
}
// Test the FFT to make sure it's working
public static void main(String[] args) {
int N = 8;
FFTItrator fft = new FFTItrator(N);
double[] window = fft.getWindow();
double[] re = new double[N];
double[] im = new double[N];
// Impulse
re[0] = 1; im[0] = 0;
for(int i=1; i<N; i++)
re[i] = im[i] = 0;
beforeAfter(fft, re, im);
// Nyquist
for(int i=0; i<N; i++) {
re[i] = Math.pow(-1, i);
im[i] = 0;
}
beforeAfter(fft, re, im);
// Single sin
for(int i=0; i<N; i++) {
re[i] = Math.cos(2*Math.PI*i / N);
im[i] = 0;
}
beforeAfter(fft, re, im);
// Ramp
for(int i=0; i<N; i++) {
re[i] = i;
im[i] = 0;
}
beforeAfter(fft, re, im);
long time = System.currentTimeMillis();
double iter = 30000;
for(int i=0; i<iter; i++)
fft.fft(re,im);
time = System.currentTimeMillis() - time;
System.out.println("Averaged " + (time/iter) + "ms per iteration");
}
protected static void beforeAfter(FFTItrator fft, double[] re, double[] im) {
System.out.println("Before: ");
printReIm(re, im);
fft.fft(re, im);
System.out.println("After: ");
printReIm(re, im);
}
protected static void printReIm(double[] re, double[] im) {
System.out.print("Re: [");
for(int i=0; i<re.length; i++)
System.out.print(((int)(re[i]*1000)/1000.0) + " ");
System.out.print("]\nIm: [");
for(int i=0; i<im.length; i++)
System.out.print(((int)(im[i]*1000)/1000.0) + " ");
System.out.println("]");
}
}
/*
Result:
Before:
Re: [1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ]
Im: [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ]
After:
Re: [1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 ]
Im: [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ]
Before:
Re: [1.0 -1.0 1.0 -1.0 1.0 -1.0 1.0 -1.0 ]
Im: [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ]
After:
Re: [0.0 0.0 0.0 0.0 8.0 0.0 0.0 0.0 ]
Im: [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ]
Before:
Re: [1.0 0.707 0.0 -0.707 -1.0 -0.707 0.0 0.707 ]
Im: [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ]
After:
Re: [0.0 4.0 0.0 0.0 0.0 0.0 0.0 4.0 ]
Im: [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ]
Before:
Re: [0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 ]
Im: [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ]
After:
Re: [28.0 -4.0 -4.0 -4.0 -4.0 -3.999 -3.999 -3.999 ]
Im: [0.0 9.656 4.0 1.656 0.0 -1.656 -4.0 -9.656 ]
Averaged 5.333333333333334E-4ms per iteration
*/