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