001: /*
002: * $RCSfile: FFT.java,v $
003: *
004: * Copyright (c) 2005 Sun Microsystems, Inc. All rights reserved.
005: *
006: * Use is subject to license terms.
007: *
008: * $Revision: 1.1 $
009: * $Date: 2005/02/11 04:56:26 $
010: * $State: Exp $
011: */
012: package com.sun.media.jai.opimage;
013:
014: import java.awt.image.DataBuffer;
015: import java.text.NumberFormat;
016: import java.util.Arrays;
017: import java.util.Locale;
018: import javax.media.jai.operator.DFTDescriptor;
019: import com.sun.media.jai.util.MathJAI;
020:
021: /**
022: * The Fast Fourier Transform (FFT) class.
023: *
024: * @since EA3
025: */
026: public class FFT {
027: /**
028: * A flag indicating that the transform is not to be scaled.
029: */
030: public static final int SCALING_NONE = DFTDescriptor.SCALING_NONE
031: .getValue();
032:
033: /**
034: * A flag indicating that the transform is to be scaled by the square
035: * root of the product of its dimensions.
036: */
037: public static final int SCALING_UNITARY = DFTDescriptor.SCALING_UNITARY
038: .getValue();
039:
040: /**
041: * A flag indicating that the transform is to be scaled by the product
042: * of its dimensions.
043: */
044: public static final int SCALING_DIMENSIONS = DFTDescriptor.SCALING_DIMENSIONS
045: .getValue();
046:
047: /** Initialization flag. */
048: protected boolean lengthIsSet = false;
049:
050: /** The sign of the exponential. */
051: protected int exponentSign;
052:
053: /** The type of scaling. */
054: protected int scaleType;
055:
056: /** The length of the FFT. */
057: protected int length;
058:
059: /** The number of bits required to represent the length. */
060: private int nbits;
061:
062: /** Indices to map between normal and bit reversed order. */
063: private int[] index;
064:
065: /** The scale factor. */
066: private double scaleFactor;
067:
068: /** Lookup table of cosines. */
069: private double[] wr;
070:
071: /** Lookup table of sines. */
072: private double[] wi;
073:
074: /** Lookup table of cosines for FCT. */
075: private double[] wrFCT;
076:
077: /** Lookup table of sines for FCT. */
078: private double[] wiFCT;
079:
080: /** Work array for real part. */
081: protected double[] real;
082:
083: /** Work array for imaginary part. */
084: protected double[] imag;
085:
086: /**
087: * Construct a new FFT object.
088: *
089: * @param negatedExponent Whether the exponent is negated.
090: * @param scaleType The type of scaling to be applied.
091: * @param length The length of the FFT; must be a positive power of 2.
092: */
093: public FFT(boolean negatedExponent, Integer scaleType, int length) {
094: // Set the exponential sign.
095: exponentSign = negatedExponent ? -1 : 1;
096:
097: // Set the scaling type.
098: this .scaleType = scaleType.intValue();
099:
100: // Set the sequence length and quantities dependent thereon.
101: setLength(length);
102: }
103:
104: /**
105: * Initialize the length-dependent fields.
106: *
107: * @param length The length of the FFT; must be a positive power of 2.
108: */
109: public void setLength(int length) {
110: // Check whether it's necessary to continue.
111: if (lengthIsSet && length == this .length) {
112: return;
113: }
114:
115: // Ensure that the length is a positive power of two.
116: if (!MathJAI.isPositivePowerOf2(length)) {
117: throw new RuntimeException(JaiI18N.getString("FFT0"));
118: }
119:
120: // Cache the length.
121: this .length = length;
122:
123: // Set the scale factor.
124: if (scaleType == SCALING_NONE) {
125: scaleFactor = 1.0;
126: } else if (scaleType == SCALING_UNITARY) {
127: scaleFactor = 1.0 / Math.sqrt(length);
128: } else if (scaleType == SCALING_DIMENSIONS) {
129: scaleFactor = 1.0 / length;
130: } else {
131: // NB: This statement should be unreachable if the scaling
132: // type is properly verified in the operation descriptor.
133: throw new RuntimeException(JaiI18N.getString("FFT1"));
134: }
135:
136: // Calculate the number of bits required to represent the length.
137: int power = 1;
138: nbits = 0;
139: while (power < length) {
140: nbits++;
141: power <<= 1;
142: }
143:
144: // Initialize the bit-reversal LUT.
145: initBitReversalLUT();
146:
147: // Calculate lookup tables of the W values.
148: calculateCoefficientLUTs();
149:
150: // Allocate work buffer memory.
151: if (!lengthIsSet || length > real.length) {
152: real = new double[length];
153: imag = new double[length];
154: }
155:
156: // Set initialization flag.
157: lengthIsSet = true;
158: }
159:
160: /**
161: * Initialize the bit-reversal lookup table.
162: */
163: private void initBitReversalLUT() {
164: // Calculate elements of index[] and fill.
165: index = new int[length];
166: for (int i = 0; i < length; ++i) {
167: int l = i;
168: int power = length >> 1;
169: int irev = 0;
170: for (int k = 0; k < nbits; ++k) {
171: int j = (l & 1);
172: if (j != 0) {
173: irev = irev + power;
174: }
175: l >>= 1;
176: power >>= 1;
177:
178: index[i] = irev;
179: }
180: }
181: }
182:
183: /**
184: * Calculate the sine and cosine lookup tables.
185: */
186: private void calculateCoefficientLUTs() {
187: wr = new double[nbits];
188: wi = new double[nbits];
189:
190: int inode = 1;
191: double cons = exponentSign * Math.PI;
192:
193: for (int bit = 0; bit < nbits; bit++) {
194: wr[bit] = Math.cos(cons / inode);
195: wi[bit] = Math.sin(cons / inode);
196: inode *= 2;
197: }
198: }
199:
200: /**
201: * Calculate the FCT sine and cosine lookup tables.
202: */
203: private void calculateFCTLUTs() {
204: wrFCT = new double[length];
205: wiFCT = new double[length];
206:
207: for (int i = 0; i < length; i++) {
208: double factor = ((i == 0) ? Math.sqrt(1.0 / length) : Math
209: .sqrt(2.0 / length));
210: double freq = Math.PI * i / (2.0 * length);
211: wrFCT[i] = factor * Math.cos(freq);
212: wiFCT[i] = factor * Math.sin(freq);
213: }
214: }
215:
216: /**
217: * Set the internal work data arrays of the FFT object.
218: *
219: * @param dataType The data type of the source data according to
220: * one of the DataBuffer TYPE_* flags. This should be either
221: * DataBuffer.TYPE_FLOAT or DataBuffer.TYPE_DOUBLE.
222: * @param realArg Float or double array of real parts.
223: * @param offsetReal Offset into the array of real parts.
224: * @param strideReal The real array stride value.
225: * @param imagArg Float or double array of imaginary parts.
226: * @param offsetImag Offset into the array of imaginary parts.
227: * @param strideImag The imaginary array stride value.
228: * @param count The number of values to copy.
229: */
230: public void setData(int dataType, Object realArg, int offsetReal,
231: int strideReal, Object imagArg, int offsetImag,
232: int strideImag, int count) {
233: // Copy the parameter arrays.
234: switch (dataType) {
235: case DataBuffer.TYPE_FLOAT: {
236: float[] realFloat = (float[]) realArg;
237: if (imagArg != null) {
238: float[] imagFloat = (float[]) imagArg;
239: if (offsetReal == offsetImag
240: && strideReal == strideImag) {
241: for (int i = 0; i < count; i++) {
242: real[i] = realFloat[offsetReal];
243: imag[i] = imagFloat[offsetReal];
244: offsetReal += strideReal;
245: }
246: } else {
247: for (int i = 0; i < count; i++) {
248: real[i] = realFloat[offsetReal];
249: imag[i] = imagFloat[offsetImag];
250: offsetReal += strideReal;
251: offsetImag += strideImag;
252: }
253: }
254: } else {
255: for (int i = 0; i < count; i++) {
256: real[i] = realFloat[offsetReal];
257: offsetReal += strideReal;
258: }
259: }
260: }
261: break;
262: case DataBuffer.TYPE_DOUBLE: {
263: double[] realDouble = (double[]) realArg;
264: if (strideReal == 1 && strideImag == 1) {
265: System
266: .arraycopy(realDouble, offsetReal, real, 0,
267: count);
268: if (imagArg != null) {
269: System.arraycopy((double[]) imagArg, offsetImag,
270: imag, 0, count);
271: }
272: } else if (imagArg != null) {
273: double[] imagDouble = (double[]) imagArg;
274: if (offsetReal == offsetImag
275: && strideReal == strideImag) {
276: for (int i = 0; i < count; i++) {
277: real[i] = realDouble[offsetReal];
278: imag[i] = imagDouble[offsetReal];
279: offsetReal += strideReal;
280: }
281: } else {
282: for (int i = 0; i < count; i++) {
283: real[i] = realDouble[offsetReal];
284: imag[i] = imagDouble[offsetImag];
285: offsetReal += strideReal;
286: offsetImag += strideImag;
287: }
288: }
289: } else {
290: for (int i = 0; i < count; i++) {
291: real[i] = realDouble[offsetReal];
292: offsetReal += strideReal;
293: }
294: }
295: }
296: break;
297: default:
298: // NB: This statement should be unreachable as the destination
299: // image is required to be a floating point type and the
300: // RasterAccessor is supposed to promote the data type of
301: // all rasters to the "minimum" data type of all source
302: // and destination rasters involved.
303: throw new RuntimeException(dataType
304: + JaiI18N.getString("FFT2"));
305: }
306:
307: // If fewer input than target points fill target with zeros.
308: if (count < length) {
309: Arrays.fill(real, count, length, 0.0);
310: if (imagArg != null) {
311: Arrays.fill(imag, count, length, 0.0);
312: }
313: }
314:
315: if (imagArg == null) {
316: Arrays.fill(imag, 0, length, 0.0);
317: }
318: }
319:
320: /**
321: * Get data from the internal work data arrays of the FFT object.
322: *
323: * @param dataType The data type of the source data according to
324: * one of the DataBuffer TYPE_* flags. This should be either
325: * DataBuffer.TYPE_FLOAT or DataBuffer.TYPE_DOUBLE.
326: * @param realArg Float or double array of real parts.
327: * @param offsetReal Offset into the array of real parts.
328: * @param strideReal The real array stride value.
329: * @param imagArg Float or double array of imaginary parts.
330: * @param offsetImag Offset into the array of imaginary parts.
331: * @param strideImag The imaginary array stride value.
332: * @param count The number of values to copy.
333: */
334: public void getData(int dataType, Object realArg, int offsetReal,
335: int strideReal, Object imagArg, int offsetImag,
336: int strideImag) {
337: switch (dataType) {
338: case DataBuffer.TYPE_FLOAT: {
339: float[] realFloat = (float[]) realArg;
340: if (imagArg != null) {
341: float[] imagFloat = (float[]) imagArg;
342: if (offsetReal == offsetImag
343: && strideReal == strideImag) {
344: for (int i = 0; i < length; i++) {
345: int idx = index[i];
346: realFloat[offsetReal] = (float) this .real[idx];
347: imagFloat[offsetReal] = (float) this .imag[idx];
348: offsetReal += strideReal;
349: }
350: } else {
351: for (int i = 0; i < length; i++) {
352: int idx = index[i];
353: realFloat[offsetReal] = (float) this .real[idx];
354: imagFloat[offsetImag] = (float) this .imag[idx];
355: offsetReal += strideReal;
356: offsetImag += strideImag;
357: }
358: }
359: } else { // imagArg == null
360: for (int i = 0; i < length; i++) {
361: realFloat[offsetReal] = (float) this .real[index[i]];
362: offsetReal += strideReal;
363: }
364: }
365: }
366: break;
367: case DataBuffer.TYPE_DOUBLE: {
368: double[] realDouble = (double[]) realArg;
369: if (imagArg != null) {
370: double[] imagDouble = (double[]) imagArg;
371: if (offsetReal == offsetImag
372: && strideReal == strideImag) {
373: for (int i = 0; i < length; i++) {
374: int idx = index[i];
375: realDouble[offsetReal] = this .real[idx];
376: imagDouble[offsetReal] = this .imag[idx];
377: offsetReal += strideReal;
378: }
379: } else {
380: for (int i = 0; i < length; i++) {
381: int idx = index[i];
382: realDouble[offsetReal] = this .real[idx];
383: imagDouble[offsetImag] = this .imag[idx];
384: offsetReal += strideReal;
385: offsetImag += strideImag;
386: }
387: }
388: } else { // imagArg == null
389: for (int i = 0; i < length; i++) {
390: realDouble[offsetReal] = this .real[index[i]];
391: offsetReal += strideReal;
392: }
393: }
394: }
395: break;
396: default:
397: // NB: This statement should be unreachable as the destination
398: // image is required to be a floating point type and the
399: // RasterAccessor is supposed to promote the data type of
400: // all rasters to the "minimum" data type of all source
401: // and destination rasters involved.
402: throw new RuntimeException(dataType
403: + JaiI18N.getString("FFT2"));
404: }
405: }
406:
407: /**
408: * Set the internal work data arrays of the FFT object for use with
409: * an FCT operation.
410: *
411: * @param dataType The data type of the source data according to
412: * one of the DataBuffer TYPE_* flags. This should be either
413: * DataBuffer.TYPE_FLOAT or DataBuffer.TYPE_DOUBLE.
414: * @param data The data as a float[] or double[].
415: * @param offset Offset into the data array.
416: * @param stride The array stride value.
417: * @param count The number of values to copy.
418: */
419: public void setFCTData(int dataType, Object data, int offset,
420: int stride, int count) {
421: // Copy the parameter arrays.
422: switch (dataType) {
423: case DataBuffer.TYPE_FLOAT: {
424: float[] realFloat = (float[]) data;
425: for (int i = 0; i < count; i++) {
426: imag[i] = realFloat[offset];
427: offset += stride;
428: }
429: for (int i = count; i < length; i++) {
430: imag[i] = 0.0;
431: }
432: int k = length - 1;
433: int j = 0;
434: for (int i = 0; i < k; i++) {
435: real[i] = imag[j++];
436: real[k--] = imag[j++];
437: }
438: }
439: break;
440: case DataBuffer.TYPE_DOUBLE: {
441: double[] realDouble = (double[]) data;
442: for (int i = 0; i < count; i++) {
443: imag[i] = realDouble[offset];
444: offset += stride;
445: }
446: for (int i = count; i < length; i++) {
447: imag[i] = 0.0;
448: }
449: int k = length - 1;
450: int j = 0;
451: for (int i = 0; i < k; i++) {
452: real[i] = imag[j++];
453: real[k--] = imag[j++];
454: }
455: }
456: break;
457: default:
458: // NB: This statement should be unreachable as the destination
459: // image is required to be a floating point type and the
460: // RasterAccessor is supposed to promote the data type of
461: // all rasters to the "minimum" data type of all source
462: // and destination rasters involved.
463: throw new RuntimeException(dataType
464: + JaiI18N.getString("FFT2"));
465: }
466:
467: // Always clear imaginary part.
468: Arrays.fill(imag, 0, length, 0.0);
469: }
470:
471: /**
472: * Get data from the internal work data arrays of the FFT object after
473: * an IFCT operation.
474: *
475: * @param dataType The data type of the source data according to
476: * one of the DataBuffer TYPE_* flags. This should be either
477: * DataBuffer.TYPE_FLOAT or DataBuffer.TYPE_DOUBLE.
478: * @param data The data as a float[] or double[].
479: * @param offset Offset into the data array.
480: * @param stride The array stride value.
481: * @param count The number of values to copy.
482: */
483: public void getFCTData(int dataType, Object data, int offset,
484: int stride) {
485: if (wrFCT == null || wrFCT.length != length) {
486: calculateFCTLUTs();
487: }
488:
489: switch (dataType) {
490: case DataBuffer.TYPE_FLOAT: {
491: float[] realFloat = (float[]) data;
492: for (int i = 0; i < length; i++) {
493: int idx = index[i];
494: realFloat[offset] = (float) (wrFCT[i] * real[idx] + wiFCT[i]
495: * imag[idx]);
496: offset += stride;
497: }
498: }
499: break;
500: case DataBuffer.TYPE_DOUBLE: {
501: double[] realDouble = (double[]) data;
502: for (int i = 0; i < length; i++) {
503: int idx = index[i];
504: realDouble[offset] = wrFCT[i] * real[idx] + wiFCT[i]
505: * imag[idx];
506: offset += stride;
507: }
508: }
509: break;
510: default:
511: // NB: This statement should be unreachable as the destination
512: // image is required to be a floating point type and the
513: // RasterAccessor is supposed to promote the data type of
514: // all rasters to the "minimum" data type of all source
515: // and destination rasters involved.
516: throw new RuntimeException(dataType
517: + JaiI18N.getString("FFT2"));
518: }
519: }
520:
521: /**
522: * Set the internal work data arrays of the FFT object for use with
523: * an IFCT operation.
524: *
525: * @param dataType The data type of the source data according to
526: * one of the DataBuffer TYPE_* flags. This should be either
527: * DataBuffer.TYPE_FLOAT or DataBuffer.TYPE_DOUBLE.
528: * @param data The data as a float[] or double[].
529: * @param offset Offset into the data array.
530: * @param stride The array stride value.
531: * @param count The number of values to copy.
532: */
533: public void setIFCTData(int dataType, Object data, int offset,
534: int stride, int count) {
535: if (wrFCT == null || wrFCT.length != length) {
536: calculateFCTLUTs();
537: }
538:
539: // Copy the parameter arrays.
540: switch (dataType) {
541: case DataBuffer.TYPE_FLOAT: {
542: float[] realFloat = (float[]) data;
543: for (int i = 0; i < count; i++) {
544: float r = realFloat[offset];
545: real[i] = r * wrFCT[i];
546: imag[i] = r * wiFCT[i];
547: offset += stride;
548: }
549: }
550: break;
551: case DataBuffer.TYPE_DOUBLE: {
552: double[] realDouble = (double[]) data;
553: for (int i = 0; i < count; i++) {
554: double r = realDouble[offset];
555: real[i] = r * wrFCT[i];
556: imag[i] = r * wiFCT[i];
557: offset += stride;
558: }
559: }
560: break;
561: default:
562: // NB: This statement should be unreachable as the destination
563: // image is required to be a floating point type and the
564: // RasterAccessor is supposed to promote the data type of
565: // all rasters to the "minimum" data type of all source
566: // and destination rasters involved.
567: throw new RuntimeException(dataType
568: + JaiI18N.getString("FFT2"));
569: }
570:
571: // If fewer input than target points fill target with zeros.
572: if (count < length) {
573: Arrays.fill(real, count, length, 0.0);
574: Arrays.fill(imag, count, length, 0.0);
575: }
576: }
577:
578: /**
579: * Get data from the internal work data arrays of the FFT object after
580: * an IFCT operation.
581: *
582: * @param dataType The data type of the source data according to
583: * one of the DataBuffer TYPE_* flags. This should be either
584: * DataBuffer.TYPE_FLOAT or DataBuffer.TYPE_DOUBLE.
585: * @param data The data as a float[] or double[].
586: * @param offset Offset into the data array.
587: * @param stride The array stride value.
588: * @param count The number of values to copy.
589: */
590: public void getIFCTData(int dataType, Object data, int offset,
591: int stride) {
592: switch (dataType) {
593: case DataBuffer.TYPE_FLOAT: {
594: float[] realFloat = (float[]) data;
595: int k = length - 1;
596: for (int i = 0; i < k; i++) {
597: realFloat[offset] = (float) real[index[i]];
598: offset += stride;
599: realFloat[offset] = (float) real[index[k--]];
600: offset += stride;
601: }
602: }
603: break;
604: case DataBuffer.TYPE_DOUBLE: {
605: double[] realDouble = (double[]) data;
606: int k = length - 1;
607: for (int i = 0; i < k; i++) {
608: realDouble[offset] = (float) real[index[i]];
609: offset += stride;
610: realDouble[offset] = (float) real[index[k--]];
611: offset += stride;
612: }
613: }
614: break;
615: default:
616: // NB: This statement should be unreachable as the destination
617: // image is required to be a floating point type and the
618: // RasterAccessor is supposed to promote the data type of
619: // all rasters to the "minimum" data type of all source
620: // and destination rasters involved.
621: throw new RuntimeException(dataType
622: + JaiI18N.getString("FFT2"));
623: }
624: }
625:
626: /**
627: * Calculate the DFT of a complex sequence using the FFT algorithm.
628: */
629: public void transform() {
630: int i, k, j, l; // Index variables
631:
632: Integer i18n = new Integer(length);
633: NumberFormat numberFormatter = NumberFormat
634: .getNumberInstance(Locale.getDefault());
635:
636: if (real.length < length || imag.length < length) {
637: throw new RuntimeException(numberFormatter.format(i18n)
638: + JaiI18N.getString("FFT3"));
639: }
640:
641: int inode = 1;
642: int ipair;
643: for (l = 0; l < nbits; ++l) {
644: double cosp = 1.0; // initial w values
645: double sinp = 0.0;
646: ipair = 2 * inode; // calc pair separation
647: for (k = 0; k < inode; ++k) {// sequence through array
648: for (i = k; i < length; i += ipair) {
649: j = i + inode; // calc other node index
650: int iIndex = index[i];
651: int jIndex = index[j];
652: double rtemp = real[jIndex] * cosp
653: - (imag[jIndex] * sinp);
654: double itemp = imag[jIndex] * cosp
655: + (real[jIndex] * sinp);
656: real[jIndex] = real[iIndex] - rtemp; // calc butterfly
657: imag[jIndex] = imag[iIndex] - itemp;
658: real[iIndex] = real[iIndex] + rtemp;
659: imag[iIndex] = imag[iIndex] + itemp;
660: }
661: double costmp = cosp;
662: cosp = cosp * wr[l] - sinp * wi[l]; // update cosp, sinp
663: sinp = costmp * wi[l] + sinp * wr[l];
664: }
665: inode = inode * 2; // new nodal dist
666: }
667:
668: if (scaleFactor != 1.0) { // multiply by non-unity scale factor
669: for (i = 0; i < length; ++i) {
670: real[i] = real[i] * scaleFactor;
671: imag[i] = imag[i] * scaleFactor;
672: }
673: }
674: }
675: }
|