|
/Users/lyon/j4p/src/math/transforms/FFTSimple1d.java
|
1 /***************************************************************************
2 // @(#)vsFFT1D.java
3 //
4 // Perform a non recursive in place FFT using doubles.
5 //
6 // @version 1.0, 4 JULY 1997
7 // @author D. Lyon, V. Silva
8 ****************************************************************************/
9 package math.transforms;
10
11 import math.Mat1;
12 import math.MathUtils;
13 import utils.Timer;
14
15 /**
16 * A 1D FFT for integral power of two FFTs
17 */
18 public class FFTSimple1d {
19 private final int FORWARD_FFT = -1;
20 private final int REVERSE_FFT = 1;
21 private int direction = FORWARD_FFT;
22 private static final float twoPI = (float) (2 * Math.PI);
23 // These arrays hold the complex 1D FFT data.
24 private float realData[] = null;
25 private float imaginaryData[] = null;
26
27 public static float[] getPsd(float audioWaveForm[]) {
28 float i[] = new float[audioWaveForm.length];
29 FFTSimple1d f = new FFTSimple1d();
30 f.computeForwardFFT(audioWaveForm, i);
31 float[] magnitudeSpectrum = f.getMagnitudeSpectrum();
32 return magnitudeSpectrum;
33 }
34
35 public void visuallyTest() {
36 System.out.println("Testing FFT engine in 1D class...");
37 int N = 8;
38 float[] in_r = new float[N];
39 float[] in_i = new float[N];
40 float[] holder_r = new float[N];
41 float[] holder_i = new float[N];
42
43 // Create a test ramp signal.
44 for (int i = 0; i < N; i++) {
45 in_r[i] = (float) i / N;
46 in_i[i] = (float) 0;
47 }
48 computeForwardFFT(in_r, in_i);
49
50 // Copy to new array because IFFT will destroy the FFT results.
51 for (int i = 0; i < N; i++) {
52 holder_r[i] = in_r[i];
53 holder_i[i] = in_i[i];
54 }
55 for (int i = 0; i < N; i++) {
56 in_r[i] = in_r[i];
57 in_i[i] = in_i[i];
58 }
59 computeBackwardFFT(in_r, in_i);
60 float[] mag = getMagnitudeSpectrum(in_r, in_i);
61 System.out.println("i x1[i] re[i] im[i] tv[i]");
62 for (int i = 0; i < N; i++) {
63 System.out.println(i +
64 "\t" + i +
65 "\t" + holder_r[i] +
66 "\t\t" +
67 holder_i[i] +
68 "\t\t" + mag[i]);
69 }
70 }
71
72 public void timeFFT(int n) {
73 // Time a 256K FFT
74 float realPart[] = new float[n];
75 float imaginaryPart[] = new float[n];
76
77 // Create a test ramp signal.
78
79 synthesizeRamp(realPart, imaginaryPart);
80 //MathUtils.print(realPart);
81
82 utils.Timer t1 = new utils.Timer();
83 runExperiment(t1, realPart, imaginaryPart);
84 }
85
86 public static void print(Object o[]) {
87 for (int i = 0; i < o.length; i++)
88 System.out.println(o[i]);
89 }
90
91 /**
92 * Destroy the input data with a linear ramp.
93 *
94 * @param realPart input datas real component
95 * @param imaginaryPart input datas' imaginary component.
96 */
97 public static final void synthesizeRamp(float[] realPart, float[] imaginaryPart) {
98 int n = realPart.length;
99 for (int i = 0; i < n; i++) {
100 realPart[i] = (float) i / n;
101 imaginaryPart[i] = (float) 0;
102 }
103 }
104
105 private void runExperiment(Timer t1,
106 float[] realPart,
107 float[] imaginaryPart) {
108 t1.start();
109 computeForwardFFT(realPart, imaginaryPart);
110 t1.stop();
111 t1.print("did an fft-radix 2 of length" + realPart.length +
112 "using float in:");
113 }
114
115 private static final void window(float r[], float w[]) {
116 for (int i = 0; i < r.length; i++)
117 r[i] = r[i] * w[i];
118 }
119
120 /**
121 * 1D FFT utility functions.
122 */
123 public void swap(int i, int numBits) {
124 float tempr;
125 int j = bitr(i, numBits);
126 tempr = realData[j];
127 realData[j] = realData[i];
128 realData[i] = tempr;
129 tempr = imaginaryData[j];
130 imaginaryData[j] = imaginaryData[i];
131 imaginaryData[i] = tempr;
132 }
133
134 // swap Zi with Zj
135 private void swapInt(int i, int j) {
136 float tempr;
137 int ti;
138 int tj;
139 ti = i - 1;
140 tj = j - 1;
141 // System.out.print("swapInt " + ti+","+tj+"\t");
142 // System.out.println(Integer.toBinaryString(ti)+","+Integer.toBinaryString(tj));
143 // System.out.println(Integer.toBinaryString(ti)+"bitr="+Integer.toBinaryString(bitr(ti)));
144 tempr = realData[tj];
145 realData[tj] = realData[ti];
146 realData[ti] = tempr;
147 tempr = imaginaryData[tj];
148 imaginaryData[tj] = imaginaryData[ti];
149 imaginaryData[ti] = tempr;
150 }
151
152 /**
153 * get the location of the maximum partial
154 *
155 * @param realPart
156 * @param imaginaryPart
157 * @return location of max(realPart[i],imaginaryPart[i])
158 */
159 public static final int getMaxPSDLocation(float realPart[],
160 float imaginaryPart[]) {
161 float max;
162 max = -0.99e30f;
163 int location = 0;
164 for (int i = 0; i < realPart.length; i++)
165 if (getMagnitude(realPart[i], imaginaryPart[i]) > max) {
166 max = realPart[i];
167 location = i;
168 }
169 return location;
170 }
171
172 /**
173 * Compute the sum of the squares of a complex number
174 *
175 * @param r real part
176 * @param imag imaginary part
177 * @return r*r + imag * imag.
178 */
179 public static final float getMagnitude(float r, float imag) {
180 return r * r + imag * imag;
181 }
182
183 float getMaxValue(float in[]) {
184 float max;
185 max = -0.99e30f;
186 for (int i = 0; i < in.length; i++)
187 if (in[i] > max)
188 max = in[i];
189 return max;
190 }
191
192 private void bitReverse2() {
193 /* bit reversal */
194 int n = realData.length;
195 int j = 1;
196 int k;
197 for (int i = 1; i < n; i++) {
198 if (i < j) {
199 swapInt(i, j);
200 } //if
201 k = n / 2;
202 while (k >= 1 && k < j) {
203 j = j - k;
204 k = k / 2;
205 }
206 j = j + k;
207 } // for
208 }
209
210 private int bitr(int j, int numBits) {
211 int ans = 0;
212 for (int i = 0; i < numBits; i++) {
213 ans = (ans << 1) + (j & 1);
214 j = j >> 1;
215 }
216 return ans;
217 }
218
219 public void computeForwardFFT(float[] in_r, float[] in_i) {
220 direction = FORWARD_FFT;
221 fft(in_r, in_i);
222 }
223
224 public void computeBackwardFFT(float[] in_r, float[] in_i) {
225 direction = REVERSE_FFT;
226 fft(in_r, in_i);
227 }
228
229 /**
230 * FFT engine.
231 */
232 public void fft(float _realPart[], float _imaginaryPart[]) {
233 int id;
234 // radix 2 number if sample, 1D of course.
235 int localN;
236 float wtemp, Wjk_r, Wjk_i, Wj_r, Wj_i;
237 float theta, tempr, tempi;
238 int numBits = MathUtils.getLogBase2(_realPart.length);
239 // Truncate input data to a power of two
240 int length = 1 << numBits; // length = 2**nu
241 int n = length;
242
243 // Copy passed references to variables to be used within
244 // transforms.fft routines & utilities
245 realData = _realPart;
246 imaginaryData = _imaginaryPart;
247 bitReverse2();
248 for (int m = 1; m <= numBits; m++) {
249 // localN = 2^m;
250 localN = 1 << m;
251 Wjk_r = 1;
252 Wjk_i = 0;
253 theta = twoPI / localN;
254 Wj_r = (float) Math.cos(theta);
255 Wj_i = (float) (direction * Math.sin(theta));
256 int nby2 = localN / 2;
257 for (int j = 0; j < nby2; j++) {
258 // This is the FFT innermost loop
259 // Any optimizations that can be made here will yield
260 // great rewards.
261 for (int k = j; k < n; k += localN) {
262 id = k + nby2;
263 tempr = Wjk_r * realData[id] - Wjk_i * imaginaryData[id];
264 tempi = Wjk_r * imaginaryData[id] + Wjk_i * realData[id];
265
266 // Zid = Zi -C
267 realData[id] = realData[k] - tempr;
268 imaginaryData[id] = imaginaryData[k] - tempi;
269 realData[k] += tempr;
270 imaginaryData[k] += tempi;
271 }
272
273 // (eq 6.23) and (eq 6.24)
274
275 wtemp = Wjk_r;
276 Wjk_r = Wj_r * Wjk_r - Wj_i * Wjk_i;
277 Wjk_i = Wj_r * Wjk_i + Wj_i * wtemp;
278 }
279 }
280 }
281
282 /**
283 * @return gets the real component from FFT.
284 */
285 public float[] getRealData() {
286 return realData;
287 }
288
289 /**
290 * @return gets the imaginary component from FFT.
291 */
292 public float[] getImaginaryData() {
293 return imaginaryData;
294 }
295
296 public float[] getMagnitudeSpectrum() {
297 float magnitudeSpectrum[] =
298 getMagnitudeSpectrum(realData, imaginaryData);
299 double s = 0.0;
300 for (int i = 0; i < magnitudeSpectrum.length; i++)
301 s = s + magnitudeSpectrum[i];
302 double v = s;
303 scale(magnitudeSpectrum, (float) (1.0 / v));
304 return magnitudeSpectrum;
305 }
306 public static void scale(final float[] a,
307 final float k) {
308 if (a == null) return;
309 for (int i = 0; i < a.length; i++) {
310 a[i] = a[i] * k;
311 }
312 }
313 /**
314 * Compute the power spectral density of the input
315 * arrays
316 *
317 * @param in_r real part of an fft
318 * @param in_i imaginary part of an fft
319 * @return the psd.
320 */
321 public float[] getMagnitudeSpectrum(float in_r[], float in_i[]) {
322 int N = in_r.length;
323 float[] mag = new float[N];
324 for (int i = 0; i < in_r.length; i++)
325 mag[i] = (float) Math.sqrt(in_r[i] * in_r[i] + in_i[i] * in_i[i]);
326 return (mag);
327 }
328
329 public static void testFFT() {
330 FFTSimple1d f = new FFTSimple1d();
331 for (int i = 0; i < 10; i++)
332 f.timeFFT(1024 * 1024);
333 }
334
335 public static void main(String args[]) {
336 testFFT();
337 }
338 }
339