|
/Users/lyon/j4p/src/math/transforms/FFT.java
|
1 /***************************************************************************
2 // @(#)vsFFT.java
3 //
4 // Perform a non recursive in place FFT.
5 //
6 // @version 1.0, 1 JULY 1997
7 // @author D. Lyon, V. Silva
8 ****************************************************************************/
9 package math.transforms;
10
11 import ip.vs.ColorUtils;
12 import ip.vs.ImageUtils;
13 import math.Mat2;
14 import utils.Timer;
15
16
17 public class FFT {
18 private final int FORWARD_FFT = -1;
19 private final int REVERSE_FFT = 1;
20 private float direction = (float) FORWARD_FFT;
21
22 static final float twoPI = (float) (2 * Math.PI);
23 private int N;
24 private int numBits;
25 private int width, height;
26 private float minPSD = 9999999;
27 private float maxPSD = -9999999;
28
29 ColorUtils CU = new ColorUtils();
30 ImageUtils iUtils = new ImageUtils();
31 public boolean DisplayLogPSD = false;
32
33 // These arrays hold the complex image data.
34 public float cR_r[][]; // Red real
35 public float cG_r[][]; // Green real
36 public float cB_r[][]; // Blue real
37 public float cR_i[][]; // Red imaginary
38 public float cG_i[][]; // Green imaginary
39 public float cB_i[][]; // Blue imaginary
40
41 // These arrays hold the complex 1D FFT data.
42 // These will hold either the full 1D FFT data, or the
43 // data from performing an FFT on one image row.
44 float r_data[] = null;
45 float i_data[] = null;
46
47 public void printStats() {
48 Mat2.printStats("cR_r", cR_r);
49 Mat2.printStats("cG_r", cG_r);
50 Mat2.printStats("cB_r", cB_r);
51 Mat2.printStats("cR_i", cR_i);
52 Mat2.printStats("cG_i", cG_i);
53 Mat2.printStats("cB_i", cB_i);
54
55 }
56
57 public void normalize() {
58 float sf = 1.0f / (cR_r.length * cR_r[0].length);
59 Mat2.scale(cR_r, sf);
60 Mat2.scale(cG_r, sf);
61 Mat2.scale(cB_r, sf);
62 Mat2.scale(cR_i, sf);
63 Mat2.scale(cG_i, sf);
64 Mat2.scale(cB_i, sf);
65 }
66
67 public void complexMult(FFT fft) {
68 for (int x = 0; x < cR_r.length; x++)
69 for (int y = 0; y < cR_r[0].length; y++) {
70 cR_r[x][y] = cR_r[x][y] * fft.cR_r[x][y];
71 cG_r[x][y] = cG_r[x][y] * fft.cG_r[x][y];
72 cB_r[x][y] = cB_r[x][y] * fft.cB_r[x][y];
73
74 cR_i[x][y] = cR_i[x][y] * fft.cR_i[x][y];
75 cG_i[x][y] = cG_i[x][y] * fft.cG_i[x][y];
76 cB_i[x][y] = cB_i[x][y] * fft.cB_i[x][y];
77 }
78
79 }
80
81 private void init() {
82 minPSD = 9999999;
83 maxPSD = -9999999;
84 }
85
86
87 public int[] getPsd() {
88 //-------------------------------------------------------------
89 // Take PSD on complex RGB values.
90 //-------------------------------------------------------------
91 // Magnitude of result.
92 float[] magnitudeR = magnitudeSpectrum(cR_r, cR_i);
93 float[] magnitudeG = magnitudeSpectrum(cG_r, cG_i);
94 float[] magnitudeB = magnitudeSpectrum(cB_r, cB_i);
95
96 double scaleFactor = 100; //255.0/Math.log(256);
97
98
99 System.out.println("Max psd = " + maxPSD);
100
101 scaleFactor = 255.0 / (Math.log(1 + maxPSD));
102 //System.out.println("Scalefactor = "+scaleFactor);
103
104
105
106 // Adjust 2-D FFT data so that the minimum PSD is
107 // based at a value close to a black pixel.
108
109 for (int i = 0; i < N; i++) {
110 magnitudeR[i] = (float)
111 (scaleFactor * Math.log(1 + magnitudeR[i]));
112 magnitudeG[i] = (float)
113 (scaleFactor * Math.log(1 + magnitudeG[i]));
114 magnitudeB[i] = (float)
115 (scaleFactor * Math.log(1 + magnitudeB[i]));
116 }
117
118 //System.out.println("Minimum PSD value: " + minPSD);
119
120 // Convert to single ARGB int array and return.
121 return (CU.imagetoInt(magnitudeR, magnitudeG, magnitudeB));
122 }
123
124 public int[] getPhaseImage() {
125 //-------------------------------------------------------------
126 // Take PSD on complex RGB values.
127 //-------------------------------------------------------------
128 // Magnitude of result.
129 float[] magnitudeR = getPhaseImage(cR_r, cR_i);
130 float[] magnitudeG = getPhaseImage(cG_r, cG_i);
131 float[] magnitudeB = getPhaseImage(cB_r, cB_i);
132
133 double scaleFactor = 100; //255.0/Math.log(256);
134
135
136 System.out.println("Max psd = " + maxPSD);
137
138 //scaleFactor = 255.0/(Math.log(1+maxPSD));
139 scaleFactor = 255;
140 System.out.println("Scalefactor = " + scaleFactor);
141
142
143
144 // Adjust 2-D FFT data so that the minimum PSD is
145 // based at a value close to a black pixel.
146
147 for (int i = 0; i < N; i++) {
148 magnitudeR[i] = (float)
149 (scaleFactor * Math.log(1 + magnitudeR[i]));
150 magnitudeG[i] = (float)
151 (scaleFactor * Math.log(1 + magnitudeG[i]));
152 magnitudeB[i] = (float)
153 (scaleFactor * Math.log(1 + magnitudeB[i]));
154 }
155
156 System.out.println("Minimum PSD value: " + minPSD);
157
158 // Convert to single ARGB int array and return.
159 return (CU.imagetoInt(magnitudeR, magnitudeG, magnitudeB));
160 }
161
162 public int[] forward2dFFT(float[] imageData_R, float[] imageData_G,
163 float[] imageData_B, int imageWidth, int imageHeight) {
164 init();
165
166 // save image size to locals.
167 width = imageWidth;
168 height = imageHeight;
169
170 // We know the size of the image now, allocate
171 // the required arrays to hold the complex image data.
172 cR_r = new float[height][width];
173 cR_i = new float[height][width];
174 cG_r = new float[height][width];
175 cG_i = new float[height][width];
176 cB_r = new float[height][width];
177 cB_i = new float[height][width];
178
179 // Total number of pixels.
180 N = width * height;
181
182 // Number of bits in one direction, i.e. 256x256 image m = 8
183 // assuming image is square, which has already been checked.
184 numBits = log2(width);
185
186 //-------------------------------------------------------------
187 // This section performs a the first forward FFT on RGB values.
188 //-------------------------------------------------------------
189 //vsTimer t1 = new vsTimer();
190 //t1.clear();
191 //t1.mark();
192
193 // FFT on RGB values 1st set.
194 // NOTE!! cX_r and cX_i are modified, they are used to hold
195 // the returned FFT complex data.
196 // RED
197 copyFloatToComplex(cR_r, cR_i, imageData_R);
198 copyFloatToComplex(cG_r, cG_i, imageData_G);
199 copyFloatToComplex(cB_r, cB_i, imageData_B);
200 // Do FFT by row.
201 for (int i = 0; i < height; i++) {
202 // Red
203 forwardFFT(cR_r[i], cR_i[i]);
204 // Green
205 forwardFFT(cG_r[i], cG_i[i]);
206 // Blue
207 forwardFFT(cB_r[i], cB_i[i]);
208 }
209
210 //System.out.println("Triple FFT on rows:");
211 // Stop the timer and report.
212 //t1.record();
213 //t1.report();
214
215 //-------------------------------------------------------------
216 // This section rotates complex RGB arrays CW 90 degrees.
217 //-------------------------------------------------------------
218 //t1.clear();
219 //t1.mark();
220 // Rotate array 90 degrees, returns a reference to a
221 // a new array, array that is passed in is lost.
222 cR_r = Rotate90(cR_r);
223 cR_i = Rotate90(cR_i);
224 cG_r = Rotate90(cG_r);
225 cG_i = Rotate90(cG_i);
226 cB_r = Rotate90(cB_r);
227 cB_i = Rotate90(cB_i);
228
229 //System.out.println("Rotate 90 degrees CW:");
230 // Stop the timer and report.
231 //t1.record();
232 //t1.report();
233
234 //-------------------------------------------------------------
235 // This section performs a the second forward FFT on RGB values.
236 //-------------------------------------------------------------
237 //t1.clear();
238 //t1.mark();
239
240 // FFT on complex data from 1st FFT.
241 // NOTE!! cX_r and cX_i that are passed in are destroyed.
242 // They are used to hold the returned FFT complex data.
243 // Do FFT by row.
244 for (int i = 0; i < height; i++) {
245 // Red
246 forwardFFT(cR_r[i], cR_i[i]);
247 // Green
248 forwardFFT(cG_r[i], cG_i[i]);
249 // Blue
250 forwardFFT(cB_r[i], cB_i[i]);
251 }
252
253 ///System.out.println("Triple FFT on columns:");
254 // Stop the timer and report.
255 //t1.record();
256 //t1.report();
257 return getPsd();
258
259 }
260
261 public int[] forward2dFFT(short[] imageData_R, short[] imageData_G,
262 short[] imageData_B, int imageWidth, int imageHeight) {
263 init();
264
265 // save image size to locals.
266 width = imageWidth;
267 height = imageHeight;
268
269 // We know the size of the image now, allocate
270 // the required arrays to hold the complex image data.
271 cR_r = new float[height][width];
272 cR_i = new float[height][width];
273 cG_r = new float[height][width];
274 cG_i = new float[height][width];
275 cB_r = new float[height][width];
276 cB_i = new float[height][width];
277
278 // Total number of pixels.
279 N = width * height;
280
281 // Number of bits in one direction, i.e. 256x256 image m = 8
282 // assuming image is square, which has already been checked.
283 numBits = log2(width);
284
285 //-------------------------------------------------------------
286 // This section performs a the first forward FFT on RGB values.
287 //-------------------------------------------------------------
288 //vsTimer t1 = new vsTimer();
289 //t1.clear();
290 //t1.mark();
291
292 // FFT on RGB values 1st set.
293 // NOTE!! cX_r and cX_i are modified, they are used to hold
294 // the returned FFT complex data.
295 // RED
296 copyShortToComplex(cR_r, cR_i, imageData_R);
297 copyShortToComplex(cG_r, cG_i, imageData_G);
298 copyShortToComplex(cB_r, cB_i, imageData_B);
299
300 // Do FFT by row.
301 for (int i = 0; i < height; i++) {
302 // Red
303 forwardFFT(cR_r[i], cR_i[i]);
304 // Green
305 forwardFFT(cG_r[i], cG_i[i]);
306 // Blue
307 forwardFFT(cB_r[i], cB_i[i]);
308 }
309
310 //System.out.println("Triple FFT on rows:");
311 // Stop the timer and report.
312 //t1.record();
313 //t1.report();
314
315 //-------------------------------------------------------------
316 // This section rotates complex RGB arrays CW 90 degrees.
317 //-------------------------------------------------------------
318 //t1.clear();
319 //t1.mark();
320 // Rotate array 90 degrees, returns a reference to a
321 // a new array, array that is passed in is lost.
322 cR_r = Rotate90(cR_r);
323 cR_i = Rotate90(cR_i);
324 cG_r = Rotate90(cG_r);
325 cG_i = Rotate90(cG_i);
326 cB_r = Rotate90(cB_r);
327 cB_i = Rotate90(cB_i);
328
329 //System.out.println("Rotate 90 degrees CW:");
330 // Stop the timer and report.
331 //t1.record();
332 //t1.report();
333
334 //-------------------------------------------------------------
335 // This section performs a the second forward FFT on RGB values.
336 //-------------------------------------------------------------
337 //t1.clear();
338 //t1.mark();
339
340 // FFT on complex data from 1st FFT.
341 // NOTE!! cX_r and cX_i that are passed in are destroyed.
342 // They are used to hold the returned FFT complex data.
343 // Do FFT by row.
344 for (int i = 0; i < height; i++) {
345 // Red
346 forwardFFT(cR_r[i], cR_i[i]);
347 // Green
348 forwardFFT(cG_r[i], cG_i[i]);
349 // Blue
350 forwardFFT(cB_r[i], cB_i[i]);
351 }
352
353 //System.out.println("Triple FFT on columns:");
354 // Stop the timer and report.
355 //t1.record();
356 //t1.report();
357 return getPsd();
358 }
359
360 int log2(double n) {
361 return (int) (Math.log(n) / Math.log(2));
362 }
363
364 public int[] reverse2dFFT() {
365 init();
366
367 // Total number of pixels.
368 N = width * height;
369
370 // Number of bits in one direction, i.e. 256x256 image m = 8
371 // assuming image is square, which has already been checked.
372 numBits = log2(width);
373
374 //-------------------------------------------------------------
375 // This section performs the first reverse FFT on
376 // the complex RGB values obtained from a forward FFT
377 // on an image.
378 //-------------------------------------------------------------
379 //vsTimer t1 = new vsTimer();
380 //t1.clear();
381 //t1.mark();
382
383 // Do IFFT by row.
384 for (int i = 0; i < height; i++) {
385 // Red
386 reverseFFT(cR_r[i], cR_i[i]);
387 // Green
388 reverseFFT(cG_r[i], cG_i[i]);
389 // Blue
390 reverseFFT(cB_r[i], cB_i[i]);
391 }
392
393 //System.out.println("Triple IFFT on rows:");
394 // Stop the timer and report.
395 //t1.record();
396 //t1.report();
397
398 //-------------------------------------------------------------
399 // This section performs a rotate CW 90 degrees to
400 // prepare the complex data for the second reverse FFT.
401 //-------------------------------------------------------------
402 //t1.clear();
403 //t1.mark();
404 // Rotate array 90 degrees, returns a reference to a
405 // a new array, array that is passed in is lost.
406 cR_r = Rotate90(cR_r);
407 cR_i = Rotate90(cR_i);
408 cG_r = Rotate90(cG_r);
409 cG_i = Rotate90(cG_i);
410 cB_r = Rotate90(cB_r);
411 cB_i = Rotate90(cB_i);
412
413 //System.out.println("Rotate 90 degrees CW:");
414 // Stop the timer and report.
415 //t1.record();
416 //t1.report();
417
418 //-------------------------------------------------------------
419 // This section performs the second reverse FFT on
420 // the complex RGB values obtained from a forward FFT
421 // on an image.
422 //-------------------------------------------------------------
423 //t1.clear();
424 //t1.mark();
425
426 // IFFT on complex data from 1st IFFT.
427 // NOTE!! cX_r and cX_i that are passed in are destroyed.
428 // They are used to hold the returned FFT complex data.
429 // Do FFT by row.
430 for (int i = 0; i < height; i++) {
431 // Red
432 reverseFFT(cR_r[i], cR_i[i]);
433 // Green
434 reverseFFT(cG_r[i], cG_i[i]);
435 // Blue
436 reverseFFT(cB_r[i], cB_i[i]);
437 }
438
439 //System.out.println("Triple IFFT on columns:");
440 // Stop the timer and report.
441 //t1.record();
442 //t1.report();
443
444
445 //-------------------------------------------------------------
446 // We now rotate the image 180 degrees otherwise it would be
447 // upside down. This is a due to the process we used to
448 // apply the FFT to columns (i.e. we rotated 90 degrees twice).
449 //-------------------------------------------------------------
450 rotateInPlace180(cR_r);
451 rotateInPlace180(cR_i);
452 rotateInPlace180(cG_r);
453 rotateInPlace180(cG_i);
454 rotateInPlace180(cB_r);
455 rotateInPlace180(cB_i);
456
457 //-------------------------------------------------------------
458 // Take PSD - power spectral density giving us original
459 // RGB values back.
460 //-------------------------------------------------------------
461 float[] realR = magnitudeSpectrum(cR_r, cR_i);
462 float[] realG = magnitudeSpectrum(cG_r, cG_i);
463 float[] realB = magnitudeSpectrum(cB_r, cB_i);
464
465 // Convert to single ARGB int array and return.
466 return (CU.imagetoInt(realR, realG, realB));
467 }
468
469
470 /**
471 * A way to visually test the 1D FFT on a small amount of data.
472 */
473 public void TestFFT() {
474 init();
475 System.out.println("Testing FFT engine in 2D class...");
476
477 numBits = 3;
478 N = 8;
479
480 float[] in_r = new float[N];
481 float[] in_i = new float[N];
482
483 float[] holder_r = new float[N];
484 float[] holder_i = new float[N];
485
486 // Create a test ramp signal.
487 for (int i = 0; i < N; i++) {
488 in_r[i] = (float) i / N;
489 in_i[i] = (float) 0;
490 }
491
492 forwardFFT(in_r, in_i);
493
494 // Copy to new array because IFFT will destroy the FFT results.
495 for (int i = 0; i < N; i++) {
496 holder_r[i] = in_r[i];
497 holder_i[i] = in_i[i];
498 }
499
500 for (int i = 0; i < N; i++) {
501 in_r[i] = in_r[i];
502 in_i[i] = in_i[i];
503 }
504
505 reverseFFT(in_r, in_i);
506 float[] mag = magnitudeSpectrum(in_r, in_i);
507
508 System.out.println("i x1[i] re[i] im[i] tv[i]");
509 for (int i = 0; i < N; i++) {
510 System.out.println(i + "\t" + i + "\t" + holder_r[i] + "\t\t" +
511 holder_i[i] + "\t\t" + mag[i]);
512 }
513
514 // Time a 256K FFT
515 in_r = new float[262144];
516 in_i = new float[262144];
517
518 // Create a test ramp signal.
519 N = 262144;
520 numBits = 18;
521
522 for (int i = 0; i < N; i++) {
523 in_r[i] = (float) i / N;
524 in_i[i] = (float) 0;
525 }
526
527 Timer t1 = new Timer();
528 t1.clear();
529 t1.start();
530
531 forwardFFT(in_r, in_i);
532
533 System.out.println("256K point 1D FFT (using float):");
534 // Stop the timer and report.
535 t1.stop();
536 t1.report();
537
538 }
539
540 /**
541 * 1D FFT utility functions.
542 */
543 public void swap(int i) {
544 float tempr;
545 int j = bitr(i);
546 String js = Integer.toBinaryString(j);
547 String is = Integer.toBinaryString(i);
548
549 // System.out.println("swap "+is+","+js);
550 // System.out.println(Integer.toBinaryString(i)+"bitr="+Integer.toBinaryString(bitr(i)));
551 tempr = r_data[j];
552 r_data[j] = r_data[i];
553 r_data[i] = tempr;
554 tempr = i_data[j];
555 i_data[j] = i_data[i];
556 i_data[i] = tempr;
557 }
558
559 // swap Zi with Zj
560 public static void swapInt(float[] i_data1, float[] r_data1, int i, int j) {
561 float tempr;
562 int ti;
563 int tj;
564 ti = i - 1;
565 tj = j - 1;
566 tempr = r_data1[tj];
567 r_data1[tj] = r_data1[ti];
568 r_data1[ti] = tempr;
569 tempr = i_data1[tj];
570 i_data1[tj] = i_data1[ti];
571 i_data1[ti] = tempr;
572 }
573
574 double getMaxValue(double in[]) {
575 double max;
576 max = -0.99e30;
577
578 for (int i = 0; i < in.length; i++)
579 if (in[i] > max)
580 max = in[i];
581 return max;
582
583 }
584
585 void bitReverse2() {
586 /* bit reversal */
587 int n = r_data.length;
588 int j = 1;
589 int k;
590
591 for (int i = 1; i < n; i++) {
592 if (i < j) {
593 swapInt(i_data, r_data, i, j);
594 } //if
595
596 k = n / 2;
597 while (k >= 1 && k < j) {
598 j = j - k;
599 k = k / 2;
600 }
601 j = j + k;
602
603 } // for
604 }
605
606 void bitReverse() {
607 /* bit reversal */
608 int n = r_data.length;
609 int j = 1;
610 int k;
611
612 for (int i = 1; i < n; i++) {
613 if (i < j) {
614 // this does not work...
615 // why?
616 swap(i);
617 } //if
618 j = bitr(i);
619
620 } // for
621 }
622
623
624 int bitr(int j) {
625 int ans = 0;
626 for (int i = 0; i < numBits; i++) {
627 ans = (ans << 1) + (j & 1);
628 j = j >> 1;
629 }
630 return ans;
631 }
632
633 public void forwardFFT(float[] in_r, float[] in_i) {
634 direction = FORWARD_FFT;
635 fft(in_r, in_i);
636 }
637
638 public void reverseFFT(float[] in_r, float[] in_i) {
639 direction = REVERSE_FFT;
640 fft(in_r, in_i);
641 }
642
643 /**
644 * FFT engine.
645 * 1D Complex fft with floats.
646 * It is Radix 2.
647 * Therefore, the size must be an integral
648 * power of 2, 2*4, 32,64, 128,256,512,1024
649 */
650 public void fft(float in_r[], float in_i[]) {
651 int id;
652 // radix 2 number if sample, 1D of course.
653 int localN;
654 float wtemp, Wjk_r, Wjk_i, Wj_r, Wj_i;
655 float theta, tempr, tempi;
656 int ti, tj;
657
658 // Truncate input data to a power of two
659 int length = 1 << numBits; // length = 2**nu
660 int n = length;
661
662 // Copy passed references to variables to be used within
663 // transforms.fft routines & utilities
664 r_data = in_r;
665 i_data = in_i;
666
667 bitReverse2();
668 for (int m = 1; m <= numBits; m++) {
669 // localN = 2^m;
670 localN = 1 << m;
671 Wjk_r = 1;
672 Wjk_i = 0;
673
674 theta = twoPI / localN;
675
676
677 Wj_r = (float) Math.cos(theta);
678 Wj_i = (float) (direction * Math.sin(theta));
679
680 int nby2 = localN / 2;
681 for (int j = 0; j < nby2; j++) {
682 // This is the FFT innermost loop
683 // Any optimizations that can be made here will yield
684 // great rewards.
685 for (int k = j; k < n; k += localN) {
686 id = k + nby2;
687 tempr = Wjk_r * r_data[id] - Wjk_i * i_data[id];
688 tempi = Wjk_r * i_data[id] + Wjk_i * r_data[id];
689
690 // Zid = Zi -C
691 r_data[id] = r_data[k] - tempr;
692 i_data[id] = i_data[k] - tempi;
693 r_data[k] += tempr;
694 i_data[k] += tempi;
695 }
696
697 // (eq 6.23) and (eq 6.24)
698
699 wtemp = Wjk_r;
700
701 Wjk_r = Wj_r * Wjk_r - Wj_i * Wjk_i;
702 Wjk_i = Wj_r * Wjk_i + Wj_i * wtemp;
703 }
704 }
705 }
706 /**
707 * FFT engine.
708 * 1D Complex fft with floats.
709 * It is Radix 2.
710 * Therefore, the size must be an integral
711 * power of 2, 2*4, 32,64, 128,256,512,1024
712 */
713 public static void fftStatic(FFT fft, float in_r[], float in_i[], int numBits) {
714 int id;
715 // radix 2 number if sample, 1D of course.
716 int localN;
717 float wtemp, Wjk_r, Wjk_i, Wj_r, Wj_i;
718 float theta, tempr, tempi;
719 int ti, tj;
720
721 // Truncate input data to a power of two
722 int length = 1 << numBits; // length = 2**nu
723 int n = length;
724
725 // Copy passed references to variables to be used within
726 // transforms.fft routines & utilities
727 fft.r_data = in_r;
728 fft.i_data = in_i;
729
730 fft.bitReverse2();
731 for (int m = 1; m <= numBits; m++) {
732 // localN = 2^m;
733 localN = 1 << m;
734 Wjk_r = 1;
735 Wjk_i = 0;
736
737 theta = twoPI / localN;
738
739
740 Wj_r = (float) Math.cos(theta);
741 Wj_i = (float) (fft.direction * Math.sin(theta));
742
743 int nby2 = localN / 2;
744 for (int j = 0; j < nby2; j++) {
745 // This is the FFT innermost loop
746 // Any optimizations that can be made here will yield
747 // great rewards.
748 for (int k = j; k < n; k += localN) {
749 id = k + nby2;
750 tempr = Wjk_r * fft.r_data[id] - Wjk_i * fft.i_data[id];
751 tempi = Wjk_r * fft.i_data[id] + Wjk_i * fft.r_data[id];
752
753 // Zid = Zi -C
754 fft.r_data[id] = fft.r_data[k] - tempr;
755 fft.i_data[id] = fft.i_data[k] - tempi;
756 fft.r_data[k] += tempr;
757 fft.i_data[k] += tempi;
758 }
759
760 // (eq 6.23) and (eq 6.24)
761
762 wtemp = Wjk_r;
763
764 Wjk_r = Wj_r * Wjk_r - Wj_i * Wjk_i;
765 Wjk_i = Wj_r * Wjk_i + Wj_i * wtemp;
766 }
767 }
768 }
769
770 public float[][] getRedReal() {
771 return (cR_r);
772 }
773
774 public float[][] getRedImaginary() {
775 return (cR_i);
776 }
777
778 public float[][] getGreenReal() {
779 return (cG_r);
780 }
781
782 public float[][] getGreenImaginary() {
783 return (cG_i);
784 }
785
786 public float[][] getBlueReal() {
787 return (cB_r);
788 }
789
790 public float[][] getBlueImaginary() {
791 return (cB_i);
792 }
793
794 /**
795 * The two arrays must be the same size.
796 */
797 private void copy2dArray(float[][] dst, float[][] src) {
798 for (int i = 0; i < height; i++) {
799 for (int j = 0; j < width; j++) {
800 dst[i][j] = src[i][j];
801 }
802 }
803 }
804
805 private void copyFloatToComplex(float[][] dst_r, float[][] dst_i,
806 float[] imageData) {
807 int k = 0;
808
809 float alternateSign = 1;
810
811 for (int i = 0; i < height; i++) {
812 for (int j = 0; j < width; j++) {
813 // Calculate (-1)**(i+j)
814 alternateSign = ((i + j) % 2 == 0) ? -1 : 1;
815
816 // 1. Put short image array into a complex pair,
817 // (-1)**(i+j) is to put the zero frequency in the
818 // center of the image when it is displayed.
819 // 2. We also take this opportunity to scale the input
820 // by N (width * height).
821 dst_r[i][j] = (float) (imageData[k++] * alternateSign / N);
822 dst_i[i][j] = (float) 0.0;
823 }
824 }
825 }
826
827 private void copyShortToComplex(float[][] dst_r, float[][] dst_i,
828 short[] imageData) {
829 int k = 0;
830
831 float alternateSign = 1;
832
833 for (int i = 0; i < height; i++)
834 for (int j = 0; j < width; j++) {
835 // Calculate (-1)**(i+j)
836 alternateSign = ((i + j) % 2 == 0) ? -1 : 1;
837
838 // 1. Put short image array into a complex pair,
839 // (-1)**(i+j) is to put the zero frequency in the
840 // center of the image when it is displayed.
841 // 2. We also take this opportunity to scale the input
842 // by N (width * height).
843 dst_r[i][j] = (float) (imageData[k++] * alternateSign / N);
844 dst_i[i][j] = (float) 0.0;
845 }
846 }
847
848 private boolean isNotInRange(int x, int low, int high) {
849 if (x < low) return true;
850 if (x >= high) return true;
851 return false;
852 }
853
854 private float[][] Rotate90(float[][] in) {
855 float[][] out = new float[height][width];
856
857 for (int i = 0; i < height; i++) {
858 for (int j = 0; j < width; j++) {
859 int x = height - j - 1;
860 int y = i;
861 if (isNotInRange(x, 0, in.length)) continue;
862 if (isNotInRange(y, 0, in[0].length)) continue;
863 out[i][j] = in[x][y];
864 }
865 }
866
867 return (out);
868 }
869
870 private void rotateInPlace180(float[][] in) {
871 float temp;
872
873 for (int i = 0; i < height / 2; i++) {
874 for (int j = 0; j < width; j++) {
875 temp = in[i][j];
876 in[i][j] = in[height - i - 1][width - j - 1];
877 in[height - i - 1][width - j - 1] = temp;
878 }
879 }
880 }
881
882 private float[] copyRealToFloat(float[][] in_r) {
883 float[] f_data = new float[N];
884 int k = 0;
885
886 for (int i = 0; i < height; i++) {
887 for (int j = 0; j < width; j++) {
888 f_data[k++] = in_r[i][j];
889 }
890 }
891
892 return (f_data);
893 }
894
895 private float[] magnitudeSpectrum(float[][] in_r, float[][] in_i) {
896 float[] mag = new float[N];
897 int k = 0;
898
899 for (int i = 0; i < height; i++) {
900 for (int j = 0; j < width; j++) {
901 mag[k] = (float) Math.sqrt(in_r[i][j] * in_r[i][j] +
902 in_i[i][j] * in_i[i][j]);
903
904 // Since we're iterating through the loop anyway, see what min is.
905 if (minPSD > mag[k])
906 minPSD = mag[k];
907 if (maxPSD < mag[k])
908 maxPSD = mag[k];
909
910 k++;
911 }
912 }
913
914 return (mag);
915 }
916
917 private float[] getPhaseImage(float[][] in_r, float[][] in_i) {
918 float phase[] = new float[N];
919 int k = 0;
920
921 for (int i = 0; i < height; i++)
922 for (int j = 0; j < width; j++) {
923 phase[k] = (float) (in_r[i][j] / in_i[i][j]);
924
925 // Since we're iterating through the loop anyway, see what min is.
926 if (minPSD > phase[k])
927 minPSD = phase[k];
928 if (maxPSD < phase[k])
929 maxPSD = phase[k];
930
931 k++;
932 }
933
934 return (phase);
935 }
936
937 public float[] magnitudeSpectrum(float[] in_r, float[] in_i) {
938 N = in_r.length;
939 float[] mag = new float[N];
940
941 for (int i = 0; i < N; i++) {
942 mag[i] = (float) Math.sqrt(in_r[i] * in_r[i] +
943 in_i[i] * in_i[i]);
944
945 // Since we're iterating through the loop anyway,
946 // find what min is.
947 if (minPSD > mag[i]) {
948 minPSD = mag[i];
949 }
950 }
951
952 return (mag);
953 }
954
955 }
956