001: package JSci.maths;
002:
003: /**
004: * The Fourier math library.
005: * This class cannot be subclassed or instantiated because all methods are static.
006: * Use <code>sort(transform(sort(...)))</code> for the discrete analogue of the continuous Fourier transform,
007: * and <code>sort(inverseTransform(sort(...)))</code> for the inverse transform.
008: * @jsci.planetmath FourierTransform
009: * @version 0.9
010: * @author Mark Hale
011: */
012: public final class FourierMath extends AbstractMath implements
013: NumericalConstants {
014: private FourierMath() {
015: }
016:
017: /**
018: * Fourier transform (2Pi convention).
019: * @param data an array containing the positive time part of the signal
020: * followed by the negative time part.
021: * @return an array containing positive frequencies in ascending order
022: * followed by negative frequencies in ascending order.
023: * @author Don Cross
024: */
025: public static Complex[] transform(final Complex data[]) {
026: final int N = data.length;
027: if (!isPowerOf2(N))
028: throw new IllegalArgumentException(
029: "The number of samples must be a power of 2.");
030:
031: final double arrayRe[] = new double[N];
032: final double arrayIm[] = new double[N];
033:
034: final int numBits = numberOfBitsNeeded(N);
035: // Simultaneous data copy and bit-reversal ordering into output
036: for (int i = 0; i < N; i++) {
037: final int j = reverseBits(i, numBits);
038: arrayRe[j] = data[i].real();
039: arrayIm[j] = data[i].imag();
040: }
041: // FFT
042: fft(arrayRe, arrayIm, TWO_PI);
043:
044: final Complex answer[] = new Complex[N];
045: for (int i = 0; i < N; i++)
046: answer[i] = new Complex(arrayRe[i], arrayIm[i]);
047: return answer;
048: }
049:
050: /**
051: * Fourier transform (2Pi convention).
052: * @param dataReal an array containing the positive real time part of the signal
053: * followed by the negative real time part.
054: * @param dataImag an array containing the positive imaginary time part of the signal
055: * followed by the negative imaginary time part.
056: * @return an array containing positive frequencies in ascending order
057: * followed by negative frequencies in ascending order.
058: * @author Don Cross
059: */
060: public static Complex[] transform(final double dataReal[],
061: final double dataImag[]) {
062: final int N = dataReal.length;
063: if (!isPowerOf2(N))
064: throw new IllegalArgumentException(
065: "The number of samples must be a power of 2.");
066:
067: final double arrayRe[] = new double[N];
068: final double arrayIm[] = new double[N];
069:
070: final int numBits = numberOfBitsNeeded(N);
071: // Simultaneous data copy and bit-reversal ordering into output
072: for (int i = 0; i < N; i++) {
073: final int j = reverseBits(i, numBits);
074: arrayRe[j] = dataReal[i];
075: arrayIm[j] = dataImag[i];
076: }
077: // FFT
078: fft(arrayRe, arrayIm, TWO_PI);
079:
080: final Complex answer[] = new Complex[N];
081: for (int i = 0; i < N; i++)
082: answer[i] = new Complex(arrayRe[i], arrayIm[i]);
083: return answer;
084: }
085:
086: /**
087: * Fourier transform (2Pi convention).
088: * @param data an array containing the positive time part of the signal
089: * followed by the negative time part.
090: * @return an array containing positive frequencies in ascending order
091: * followed by negative frequencies in ascending order.
092: */
093: public static Complex[] transform(final double data[]) {
094: final int N = data.length;
095: if (!isPowerOf2(N))
096: throw new IllegalArgumentException(
097: "The number of samples must be a power of 2.");
098:
099: final double arrayRe[] = new double[N];
100: final double arrayIm[] = new double[N];
101:
102: final int numBits = numberOfBitsNeeded(N);
103: // Simultaneous data copy and bit-reversal ordering into output
104: for (int i = 0; i < N; i++) {
105: final int j = reverseBits(i, numBits);
106: arrayRe[j] = data[i];
107: }
108: // FFT
109: fft(arrayRe, arrayIm, TWO_PI);
110:
111: final Complex answer[] = new Complex[N];
112: for (int i = 0; i < N; i++)
113: answer[i] = new Complex(arrayRe[i], arrayIm[i]);
114: return answer;
115: }
116:
117: /**
118: * Inverse Fourier transform (-2Pi convention).
119: * @param data an array containing positive frequencies in ascending order
120: * followed by negative frequencies in ascending order.
121: * @return an array containing the positive time part of the signal
122: * followed by the negative time part.
123: * @author Don Cross
124: */
125: public static Complex[] inverseTransform(final Complex data[]) {
126: final int N = data.length;
127: if (!isPowerOf2(N))
128: throw new IllegalArgumentException(
129: "Data length must be a power of 2.");
130:
131: final double arrayRe[] = new double[N];
132: final double arrayIm[] = new double[N];
133:
134: final int numBits = numberOfBitsNeeded(N);
135: // Simultaneous data copy and bit-reversal ordering into output
136: for (int i = 0; i < N; i++) {
137: final int j = reverseBits(i, numBits);
138: arrayRe[j] = data[i].real();
139: arrayIm[j] = data[i].imag();
140: }
141: // inverse FFT
142: fft(arrayRe, arrayIm, -TWO_PI);
143: // Normalize
144: final Complex answer[] = new Complex[N];
145: final double denom = N;
146: for (int i = 0; i < N; i++)
147: answer[i] = new Complex(arrayRe[i] / denom, arrayIm[i]
148: / denom);
149: return answer;
150: }
151:
152: /**
153: * Inverse Fourier transform (-2Pi convention).
154: * @param dataReal an array containing positive real frequencies in ascending order
155: * followed by negative real frequencies in ascending order.
156: * @param dataImag an array containing positive imaginary frequencies in ascending order
157: * followed by negative imaginary frequencies in ascending order.
158: * @return an array containing the positive time part of the signal
159: * followed by the negative time part.
160: * @author Don Cross
161: */
162: public static Complex[] inverseTransform(final double dataReal[],
163: final double dataImag[]) {
164: final int N = dataReal.length;
165: if (!isPowerOf2(N))
166: throw new IllegalArgumentException(
167: "Data length must be a power of 2.");
168:
169: final double arrayRe[] = new double[N];
170: final double arrayIm[] = new double[N];
171:
172: final int numBits = numberOfBitsNeeded(N);
173: // Simultaneous data copy and bit-reversal ordering into output
174: for (int i = 0; i < N; i++) {
175: final int j = reverseBits(i, numBits);
176: arrayRe[j] = dataReal[i];
177: arrayIm[j] = dataImag[i];
178: }
179: // inverse FFT
180: fft(arrayRe, arrayIm, -TWO_PI);
181: // Normalize
182: final Complex answer[] = new Complex[N];
183: final double denom = N;
184: for (int i = 0; i < N; i++)
185: answer[i] = new Complex(arrayRe[i] / denom, arrayIm[i]
186: / denom);
187: return answer;
188: }
189:
190: /**
191: * Inverse Fourier transform (-2Pi convention).
192: * @param data an array containing positive frequencies in ascending order
193: * followed by negative frequencies in ascending order.
194: * @return an array containing the positive time part of the signal
195: * followed by the negative time part.
196: */
197: public static Complex[] inverseTransform(final double data[]) {
198: final int N = data.length;
199: if (!isPowerOf2(N))
200: throw new IllegalArgumentException(
201: "Data length must be a power of 2.");
202:
203: final double arrayRe[] = new double[N];
204: final double arrayIm[] = new double[N];
205:
206: final int numBits = numberOfBitsNeeded(N);
207: // Simultaneous data copy and bit-reversal ordering into output
208: for (int i = 0; i < N; i++) {
209: final int j = reverseBits(i, numBits);
210: arrayRe[j] = data[i];
211: }
212: // inverse FFT
213: fft(arrayRe, arrayIm, -TWO_PI);
214: // Normalize
215: final Complex answer[] = new Complex[N];
216: final double denom = N;
217: for (int i = 0; i < N; i++)
218: answer[i] = new Complex(arrayRe[i] / denom, arrayIm[i]
219: / denom);
220: return answer;
221: }
222:
223: /**
224: * Common FFT code.
225: * @param twoPi TWO_PI for transform, -TWO_PI for inverse transform.
226: */
227: private static void fft(double arrayRe[], double arrayIm[],
228: final double twoPi) {
229: final int N = arrayRe.length;
230: int blockEnd = 1;
231: for (int blockSize = 2; blockSize <= N; blockSize <<= 1) {
232: final double deltaAngle = twoPi / blockSize;
233: double alpha = Math.sin(0.5 * deltaAngle);
234: alpha = 2.0 * alpha * alpha;
235: final double beta = Math.sin(deltaAngle);
236: for (int i = 0; i < N; i += blockSize) {
237: double angleRe = 1.0;
238: double angleIm = 0.0;
239: for (int j = i, n = 0; n < blockEnd; j++, n++) {
240: final int k = j + blockEnd;
241: // tmp = angle*array[k]
242: double tmpRe = angleRe * arrayRe[k] - angleIm
243: * arrayIm[k];
244: double tmpIm = angleRe * arrayIm[k] + angleIm
245: * arrayRe[k];
246: arrayRe[k] = arrayRe[j] - tmpRe;
247: arrayIm[k] = arrayIm[j] - tmpIm;
248: arrayRe[j] += tmpRe;
249: arrayIm[j] += tmpIm;
250: // angle = angle - (a-bi)*angle
251: tmpRe = alpha * angleRe + beta * angleIm;
252: tmpIm = alpha * angleIm - beta * angleRe;
253: angleRe -= tmpRe;
254: angleIm -= tmpIm;
255: }
256: }
257: blockEnd = blockSize;
258: }
259: }
260:
261: /**
262: * Returns true if x is a power of 2.
263: * @author Don Cross
264: */
265: private static boolean isPowerOf2(final int x) {
266: final int BITS_PER_WORD = 32;
267: for (int i = 1, y = 2; i < BITS_PER_WORD; i++, y <<= 1) {
268: if (x == y)
269: return true;
270: }
271: return false;
272: }
273:
274: /**
275: * Number of bits needed.
276: * @author Don Cross
277: */
278: private static int numberOfBitsNeeded(final int pwrOf2) {
279: if (pwrOf2 < 2)
280: throw new IllegalArgumentException();
281: for (int i = 0;; i++) {
282: if ((pwrOf2 & (1 << i)) > 0)
283: return i;
284: }
285: }
286:
287: /**
288: * Reverse bits.
289: * @author Don Cross
290: */
291: private static int reverseBits(int index, final int numBits) {
292: int i, rev;
293: for (i = rev = 0; i < numBits; i++) {
294: rev = (rev << 1) | (index & 1);
295: index >>= 1;
296: }
297: return rev;
298: }
299:
300: /**
301: * Sorts the output from the Fourier transfom methods into
302: * ascending frequency/time order.
303: */
304: public static Complex[] sort(final Complex output[]) {
305: final Complex ret[] = new Complex[output.length];
306: final int Nby2 = output.length / 2;
307: for (int i = 0; i < Nby2; i++) {
308: ret[Nby2 + i] = output[i];
309: ret[i] = output[Nby2 + i];
310: }
311: return ret;
312: }
313:
314: public static double[] sort(final double input[]) {
315: final double ret[] = new double[input.length];
316: final int Nby2 = input.length / 2;
317: for (int i = 0; i < Nby2; i++) {
318: ret[Nby2 + i] = input[i];
319: ret[i] = input[Nby2 + i];
320: }
321: return ret;
322: }
323: }
|