00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 package com.meapsoft;
00024
00025
00033 public class FFT {
00034
00035 int n, m;
00036
00037
00038 double[] cos;
00039 double[] sin;
00040
00041 double[] window;
00042
00043 public FFT(int n) {
00044 this.n = n;
00045 this.m = (int)(Math.log(n) / Math.log(2));
00046
00047
00048 if(n != (1<<m))
00049 throw new RuntimeException("FFT length must be power of 2");
00050
00051
00052 cos = new double[n/2];
00053 sin = new double[n/2];
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066 for(int i=0; i<n/2; i++) {
00067 cos[i] = Math.cos(-2*Math.PI*i/n);
00068 sin[i] = Math.sin(-2*Math.PI*i/n);
00069 }
00070
00071 makeWindow();
00072 }
00073
00074 protected void makeWindow() {
00075
00076
00077 window = new double[n];
00078 for(int i = 0; i < window.length; i++)
00079 window[i] = 0.42 - 0.5 * Math.cos(2*Math.PI*i/(n-1))
00080 + 0.08 * Math.cos(4*Math.PI*i/(n-1));
00081 }
00082
00083 public double[] getWindow() {
00084 return window;
00085 }
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107 public void fft(double[] x, double[] y)
00108 {
00109 int i,j,k,n1,n2,a;
00110 double c,s,e,t1,t2;
00111
00112
00113
00114 j = 0;
00115 n2 = n/2;
00116 for (i=1; i < n - 1; i++) {
00117 n1 = n2;
00118 while ( j >= n1 ) {
00119 j = j - n1;
00120 n1 = n1/2;
00121 }
00122 j = j + n1;
00123
00124 if (i < j) {
00125 t1 = x[i];
00126 x[i] = x[j];
00127 x[j] = t1;
00128 t1 = y[i];
00129 y[i] = y[j];
00130 y[j] = t1;
00131 }
00132 }
00133
00134
00135 n1 = 0;
00136 n2 = 1;
00137
00138 for (i=0; i < m; i++) {
00139 n1 = n2;
00140 n2 = n2 + n2;
00141 a = 0;
00142
00143 for (j=0; j < n1; j++) {
00144 c = cos[a];
00145 s = sin[a];
00146 a += 1 << (m-i-1);
00147
00148 for (k=j; k < n; k=k+n2) {
00149 t1 = c*x[k+n1] - s*y[k+n1];
00150 t2 = s*x[k+n1] + c*y[k+n1];
00151 x[k+n1] = x[k] - t1;
00152 y[k+n1] = y[k] - t2;
00153 x[k] = x[k] + t1;
00154 y[k] = y[k] + t2;
00155 }
00156 }
00157 }
00158 }
00159
00160
00161
00162
00163
00164 public static void main(String[] args) {
00165 int N = 8;
00166
00167 FFT fft = new FFT(N);
00168
00169 double[] window = fft.getWindow();
00170 double[] re = new double[N];
00171 double[] im = new double[N];
00172
00173
00174 re[0] = 1; im[0] = 0;
00175 for(int i=1; i<N; i++)
00176 re[i] = im[i] = 0;
00177 beforeAfter(fft, re, im);
00178
00179
00180 for(int i=0; i<N; i++) {
00181 re[i] = Math.pow(-1, i);
00182 im[i] = 0;
00183 }
00184 beforeAfter(fft, re, im);
00185
00186
00187 for(int i=0; i<N; i++) {
00188 re[i] = Math.cos(2*Math.PI*i / N);
00189 im[i] = 0;
00190 }
00191 beforeAfter(fft, re, im);
00192
00193
00194 for(int i=0; i<N; i++) {
00195 re[i] = i;
00196 im[i] = 0;
00197 }
00198 beforeAfter(fft, re, im);
00199
00200 long time = System.currentTimeMillis();
00201 double iter = 30000;
00202 for(int i=0; i<iter; i++)
00203 fft.fft(re,im);
00204 time = System.currentTimeMillis() - time;
00205 System.out.println("Averaged " + (time/iter) + "ms per iteration");
00206 }
00207
00208 protected static void beforeAfter(FFT fft, double[] re, double[] im) {
00209 System.out.println("Before: ");
00210 printReIm(re, im);
00211 fft.fft(re, im);
00212 System.out.println("After: ");
00213 printReIm(re, im);
00214 }
00215
00216 protected static void printReIm(double[] re, double[] im) {
00217 System.out.print("Re: [");
00218 for(int i=0; i<re.length; i++)
00219 System.out.print(((int)(re[i]*1000)/1000.0) + " ");
00220
00221 System.out.print("]\nIm: [");
00222 for(int i=0; i<im.length; i++)
00223 System.out.print(((int)(im[i]*1000)/1000.0) + " ");
00224
00225 System.out.println("]");
00226 }
00227 }