|
/Users/lyon/j4p/src/math/transforms/DFT.java
|
1 package math.transforms;
2
3
4
5 import sound.Oscillator;
6
7 import java.awt.*;
8
9 public class DFT extends Frame {
10
11 private double r_data[] = null;
12 private double i_data[] = null;
13
14
15 private String title = "FFT :";
16
17 // if forward = true then perform a forward fft
18 // else perform a backward fft;
19 private boolean forward = true;
20 private int nu; // number of bits needed to represent the index for r_data
21 private static final float twoPI = (float) (2 * Math.PI);
22
23 public DFT(int N) {
24 r_data = new double[N];
25 i_data = new double[N];
26 }
27
28 public DFT() {
29 }
30
31
32
33 public void graphs(String t) {
34 setTitle(t);
35 }
36
37 public void setTitle(String t) {
38 title = t;
39 }
40
41 // swap Zi with Zj
42 private void swapInt(int i, int j) {
43 double tempr;
44 int ti;
45 int tj;
46 ti = i - 1;
47 tj = j - 1;
48 tempr = r_data[tj];
49 r_data[tj] = r_data[ti];
50 r_data[ti] = tempr;
51 tempr = i_data[tj];
52 i_data[tj] = i_data[ti];
53 i_data[ti] = tempr;
54 }
55
56 public static double getMaxValue(double in[]) {
57 double max;
58 max = -0.99e30;
59 for (int i = 0; i < in.length; i++)
60 if (in[i] > max)
61 max = in[i];
62 return max;
63
64 }
65
66 private void normalizeAndTruncateInput(double in[]) {
67 double data_max = getMaxValue(in);
68 /* copy over normalized input */
69 for (int i = 0; i < r_data.length; i++) {
70 r_data[i] = in[i];
71 }
72 }
73
74 private void bitReverse2() {
75 System.out.println("fft: bit reversal");
76 /* bit reversal */
77 int n = r_data.length;
78 int j = 1;
79
80 int k;
81
82 for (int i = 1; i < n; i++) {
83
84 if (i < j) swapInt(i, j);
85 k = n / 2;
86 while (k >= 1 && k < j) {
87
88 j = j - k;
89 k = k / 2;
90 }
91 j = j + k;
92 } // for
93 }
94
95 private int bitr(int j) {
96 int ans = 0;
97 for (int i = 0; i < nu; i++) {
98 ans = (ans << 1) + (j & 1);
99 j = j >> 1;
100 }
101 return ans;
102 }
103
104 public void reverseFFT(double in_r[], double in_i[]) {
105 forward = false;
106 forwardFFT(in_r, in_i);
107 forward = true;
108 //centering(in_r);
109 }
110
111 private void centering(double r[]) {
112 int s = 1;
113 for (int i = 0; i < r.length; i++) {
114 s = -s;
115 r[i] *= s;
116 }
117 }
118
119 public void forwardFFT(double in_r[], double in_i[]) {
120 int id;
121
122 int localN;
123 double wtemp, Wjk_r, Wjk_i, Wj_r, Wj_i;
124 double theta, tempr, tempi;
125 int ti, tj;
126
127 int numBits = log2(in_r.length);
128 if (forward) {
129 //centering(in_r);
130 }
131
132
133 // Truncate input data to a power of two
134 int length = 1 << numBits; // length = 2**nu
135 int n = length;
136 int nby2;
137
138 // Copy passed references to variables to be used within
139 // fft routines & utilities
140 r_data = in_r;
141 i_data = in_i;
142
143 bitReverse2();
144 for (int m = 1; m <= numBits; m++) {
145 // localN = 2^m;
146 localN = 1 << m;
147
148 nby2 = localN / 2;
149 Wjk_r = 1;
150 Wjk_i = 0;
151
152 theta = Math.PI / nby2;
153
154 // for recursive comptutation of sine and cosine
155 Wj_r = Math.cos(theta);
156 Wj_i = -Math.sin(theta);
157 if (forward == false) {
158 Wj_i = -Wj_i;
159 }
160
161
162 for (int j = 0; j < nby2; j++) {
163 // This is the FFT innermost loop
164 // Any optimizations that can be made here will yield
165 // great rewards.
166 for (int k = j; k < n; k += localN) {
167 id = k + nby2;
168 tempr = Wjk_r * r_data[id] - Wjk_i * i_data[id];
169 tempi = Wjk_r * i_data[id] + Wjk_i * r_data[id];
170
171 // Zid = Zi -C
172 r_data[id] = r_data[k] - tempr;
173 i_data[id] = i_data[k] - tempi;
174 r_data[k] += tempr;
175 i_data[k] += tempi;
176 }
177
178 // (eq 6.23) and (eq 6.24)
179
180 wtemp = Wjk_r;
181
182 Wjk_r = Wj_r * Wjk_r - Wj_i * Wjk_i;
183 Wjk_i = Wj_r * Wjk_i + Wj_i * wtemp;
184 }
185 }
186 // normalize output of fft.
187 if (forward)
188 for (int i = 0; i < r_data.length; i++) {
189 r_data[i] = r_data[i] / (double) n;
190 i_data[i] = i_data[i] / (double) n;
191 }
192 in_r = r_data;
193 in_r = i_data;
194 }
195
196 public static int log2(int n) {
197 return (int) (Math.log(n) / Math.log(2.0));
198 }
199
200 public static double[] arrayCopy(double[] in) {
201 int N = in.length;
202 double out[] = new double[N];
203 for (int i = 0; i < N; i++) {
204 out[i] = in[i];
205 }
206 return out;
207 }
208
209 public double[] dft(double v[]) {
210 int N = v.length;
211
212 double t_img, t_real;
213
214 double twoPikOnN;
215 double twoPijkOnN;
216
217 // how many bits do we need?
218 N = log2(N);
219 //Truncate input data to a power of two
220 // length = 2**(number of bits).
221 N = 1 << N;
222
223 double twoPiOnN = 2 * Math.PI / N;
224 // We truncate to a power of two so that
225 // we can compare execution times with the FFT.
226 // DFT generally does not need to truncate its input.
227
228 r_data = new double[N];
229 i_data = new double[N];
230 double psd[] = new double[N];
231
232
233 System.out.println("Executing DFT on " + N + " points...");
234
235 for (int k = 0; k < N; k++) {
236
237 twoPikOnN = twoPiOnN * k;
238
239
240 for (int j = 0; j < N; j++) {
241 twoPijkOnN = twoPikOnN * j;
242 r_data[k] += v[j] * Math.cos(twoPijkOnN);
243 i_data[k] -= v[j] * Math.sin(twoPijkOnN);
244 }
245 }
246 normalizeData();
247 return (getPowerSpectralDensity());
248 }
249
250 private void normalizeData() {
251 int N = r_data.length;
252 for (int k = 0; k < N; k++) {
253 r_data[k] /= N;
254 i_data[k] /= N;
255 }
256
257 }
258
259 public double[] getReal() {
260 return r_data;
261 }
262
263 public double[] getImaginary() {
264 return i_data;
265 }
266
267
268 public double[] getPowerSpectralDensity() {
269 double[] psd = new double[r_data.length];
270 for (int k = 0; k < r_data.length; k++) {
271 psd[k] =
272 r_data[k] * r_data[k] +
273 i_data[k] * i_data[k];
274
275 }
276 return psd;
277 }
278
279 // assume that r_data and i_data are
280 // set. Also assume that the real
281 // value is to be returned
282 public double[] ifft() {
283 int i, m, j, id;
284 int N; // the radix 2 number of samples
285 double wtemp, wr, wpr, wpi, wi, theta, tempr, tempi;
286
287 // length is the number of input samples
288 int length = r_data.length;
289
290 // how many bits do we need?
291 nu = (int) (Math.log(length) / Math.log(2.0));
292 //Truncate input data to a power of two
293 length = 1 << nu; // length = 2**nu
294
295 int ti, tj;
296
297 int n = length;
298
299 for (m = 1; m <= nu; m++) {
300
301 // k = 2^m;
302 N = 1 << m;
303 theta = twoPI / N;
304 // theta = - 2Pi/(2^m)
305
306 wr = 1.0;
307 wi = 0.0;
308 wpr = Math.cos(theta);
309
310 // ifft uses - sin(theta);
311 wpi = -Math.sin(theta);
312
313
314 for (j = 1; j <= N / 2; j++) {
315 for (i = j; i <= n; i = i + N) {
316
317 id = i + N / 2;
318 tempr = wr * r_data[id - 1] - wi * i_data[id - 1];
319 tempi = wr * i_data[id - 1] + wi * r_data[id - 1];
320
321 // Zid-1 = Zi-1 - C(tempr,tempi)
322 r_data[id - 1] = r_data[i - 1] - tempr;
323 i_data[id - 1] = i_data[i - 1] - tempi;
324 r_data[i - 1] += tempr;
325 i_data[i - 1] += tempi;
326 }
327 wtemp = wr;
328 // W = W * WP
329 // W = (wr + i wi) (wpr + i wpi)
330 // W = wr * wpr - wi * wpi + i (wi * wpr + wr * wpi)
331 wr = wr * wpr - wi * wpi;
332 wi = wi * wpr + wtemp * wpi;
333 }
334 }
335 return (arrayCopy(r_data));
336 }
337
338
339 // assume that r_data and i_data are
340 // set. Also assume that the real
341 // value is to be returned
342 public double[] idft() {
343 int N = r_data.length;
344 double twoPiOnN = 2 * Math.PI / N;
345
346 double twoPikOnN;
347 double twoPijkOnN;
348
349 double v[] = new double[N];
350
351 System.out.println("Executing IDFT on " + N + " points...");
352
353
354 for (int k = 0; k < N; k++) {
355 twoPikOnN = twoPiOnN * k;
356 for (int j = 0; j < N; j++) {
357 twoPijkOnN = twoPikOnN * j;
358 v[k] += r_data[j] * Math.cos(twoPijkOnN)
359 - i_data[j] * Math.sin(twoPijkOnN);
360
361 }
362 }
363 return (v);
364 }
365
366 public static void printArray(double[] v, String title) {
367 System.out.println(title);
368 for (int i = 0; i < v.length; i++) {
369 System.out.println("v[" + i + "]=" + v[i]);
370 }
371 }
372
373
374 public void printArrays(String title) {
375 System.out.println(title);
376 for (int i = 0; i < r_data.length; i++) {
377 System.out.println("[" + i + "]=("
378 + r_data[i] + "," + i_data[i] + ")");
379 }
380
381 }
382
383 public void printReal(String title) {
384 System.out.println(title);
385 for (int i = 0; i < r_data.length; i++) {
386 System.out.println(
387 "r_data[" + i + "]="
388 + r_data[i]);
389 }
390
391 }
392
393 public static void main(String args[]) {
394 Oscillator.testAudioDft();
395
396 }
397
398 public static void testPSD() {
399
400
401 DFT f = new DFT();
402
403 int N = 8;
404 int numBits = log2(N);
405
406 double x1[] = new double[N];
407 for (int j = 0; j < N; j++)
408 x1[j] = j;
409
410
411 double[] in_r = new double[N];
412 double[] in_i = new double[N];
413
414 // copy test signal.
415 in_r = arrayCopy(x1);
416
417 f.forwardFFT(in_r, in_i);
418 f.printArrays("After the FFT");
419 double psd[] = f.getPowerSpectralDensity();
420 DFT.printArray(psd, "The psd");
421 }
422
423 public static void timeFFT() {
424 int N = 2048;
425 DFT f = new DFT(N);
426 int numBits = log2(N);
427
428 double x1[] = new double[N];
429 for (int j = 0; j < N; j++)
430 x1[j] = j;
431
432
433 double[] in_r = new double[N];
434 double[] in_i = new double[N];
435
436 double[] fftResult_r = new double[N];
437 double[] fftResult_i = new double[N];
438
439 // copy test signal.
440 in_r = arrayCopy(x1);
441
442
443 f.forwardFFT(in_r, in_i);
444
445
446 System.out.print("Time for " + N + "point fft");
447
448
449 // Copy to new array because IFFT will
450 // destroy the FFT results.
451 fftResult_r = arrayCopy(in_r);
452 fftResult_i = arrayCopy(in_i);
453
454
455 f.reverseFFT(in_r, in_i);
456
457 System.out.print("Time for " + N + "point ifft");
458
459
460
461 }
462
463 public static void testFFT() {
464 System.out.println("Starting 1D FFT test...");
465 DFT f = new DFT();
466
467 int N = 8;
468 int numBits = log2(N);
469
470 double x1[] = new double[N];
471 for (int j = 0; j < N; j++)
472 x1[j] = j;
473
474
475 double[] in_r = new double[N];
476 double[] in_i = new double[N];
477
478 double[] fftResult_r = new double[N];
479 double[] fftResult_i = new double[N];
480
481 // copy test signal.
482 in_r = arrayCopy(x1);
483
484 f.forwardFFT(in_r, in_i);
485
486 // Copy to new array because IFFT will
487 // destroy the FFT results.
488 fftResult_r = arrayCopy(in_r);
489 fftResult_i = arrayCopy(in_i);
490
491
492 f.reverseFFT(in_r, in_i);
493
494 System.out.println("j\tx1[j]\tre[j]\tim[j]\tv[j]");
495 for (int i = 0; i < N; i++) {
496 System.out.println(
497 i + "\t" +
498 x1[i] + "\t" +
499 fftResult_r[i] + "\t" +
500 fftResult_i[i] + "\t" +
501 in_r[i]);
502 }
503
504 }
505
506 public static void print(double o[]){
507 for (int i=0; i < o.length; i++)
508 System.out.println(o[i]);
509 }
510 public static void testDFT() {
511
512 int N = 8;
513 DFT f = new DFT(N);
514
515 double v[];
516 double x1[] = new double[N];
517 for (int j = 0; j < N; j++)
518 x1[j] = j;
519
520 // take dft
521 f.dft(x1);
522 v = f.idft();
523 System.out.println("j\tx1[j]\tre[j]\tim[j]\t v[j]");
524 for (int j = 0; j < N; j++)
525 System.out.println(
526 j + "\t" +
527 x1[j] + "\t" +
528 f.r_data[j] + "\t" +
529 f.i_data[j] + "\t" +
530 v[j]);
531
532 }
533
534 }
535