1: // Code to implement a decently performing FFT for complex and real valued
2: // signals. See www.lomont.org for a derivation of the relevant algorithms
3: // from first principles. Copyright Chris Lomont 2010.
4: // This code and any ports are free for all to use for any reason as long
5: // as this header is left in place.
6: using System;
7: using System.Collections.Generic;
8: using System.Linq;
9:
10: namespace Lomont
11: {
12: /// <summary>
13: /// Represent a class that performs real or complex valued Fast Fourier
14: /// Transforms. Instantiate it and use the FFT or TableFFT methods to
15: /// compute complex to complex FFTs. Use FFTReal for real to complex
16: /// FFTs which are much faster than standard FFTs.
17: /// </summary>
18: class LomontFFT
19: {
20: /// <summary>
21: /// Compute the forward or inverse Fourier Transform of data, with
22: /// data containing complex valued data as alternating real and
23: /// imaginary parts. The length must be a power of 2.
24: /// </summary>
25: /// <param name="data">The complex data stored as alternating real
26: /// and imaginary parts</param>
27: /// <param name="forward">true for a forward transform, false for
28: /// inverse transform</param>
29: public void FFT(double[] data, bool forward)
30: {
31: int n = data.Length;
32: // checks n is a power of 2 in 2's complement format
33: if ((n & (n - 1)) != 0)
34: throw new ArgumentException(
35: "data length " + n + " in FFT is not a power of 2");
36: n /= 2; // n is the number of samples
37:
38: Reverse(data, n); // bit index data reversal
39:
40: // do transform: so single point transforms, then doubles, etc.
41: double sign = forward ? 1 : -1;
42: int mmax = 1;
43: while (n > mmax)
44: {
45: int istep = 2 * mmax;
46: double theta = sign * Math.PI / mmax;
47: double wr = 1, wi = 0;
48: double wpr = Math.Cos(theta);
49: double wpi = Math.Sin(theta);
50: for (int m = 0; m < istep; m += 2)
51: {
52: for (int k = m; k < 2 * n; k += 2 * istep)
53: {
54: int j = k + istep;
55: double tempr = wr * data[j] - wi * data[j + 1];
56: double tempi = wi * data[j] + wr * data[j + 1];
57: data[j] = data[k] - tempr;
58: data[j + 1] = data[k + 1] - tempi;
59: data[k] = data[k] + tempr;
60: data[k + 1] = data[k + 1] + tempi;
61: }
62: double t = wr; // trig recurrence
63: wr = wr * wpr - wi * wpi;
64: wi = wi * wpr + t * wpi;
65: }
66: mmax = istep;
67: }
68:
69: // inverse scaling in the backward case
70: if (!forward)
71: {
72: double scale = 1.0 / n;
73: for (int i = 0; i < 2 * n; ++i)
74: data[i] *= scale;
75: }
76: }
77:
78:
79: /// <summary>
80: /// Compute the forward or inverse Fourier Transform of data, with data
81: /// containing complex valued data as alternating real and imaginary
82: /// parts. The length must be a power of 2. This method caches values
83: /// and should be slightly faster on repeated uses than then FFT method.
84: /// It is also slightly more accurate.
85: /// </summary>
86: /// <param name="data">The complex data stored as alternating real
87: /// and imaginary parts</param>
88: /// <param name="forward">true for a forward transform, false for
89: /// inverse transform</param>
90: public void TableFFT(double[] data, bool forward)
91: {
92: int n = data.Length;
93: // checks n is a power of 2 in 2's complement format
94: if ((n & (n - 1)) != 0)
95: throw new ArgumentException(
96: "data length " + n + " in FFT is not a power of 2"
97: );
98: n /= 2; // n is the number of samples
99:
100: Reverse(data, n); // bit index data reversal
101:
102: // make table if needed
103: if ((cosTable == null) || (cosTable.Count != n))
104: Initialize(n);
105:
106: // do transform: so single point transforms, then doubles, etc.
107: double sign = forward ? 1 : -1;
108: int mmax = 1;
109: int tptr = 0;
110: while (n > mmax)
111: {
112: int istep = 2 * mmax;
113: double theta = sign * Math.PI / mmax;
114: for (int m = 0; m < istep; m += 2)
115: {
116: double wr = cosTable[tptr];
117: double wi = sign*sinTable[tptr++];
118: for (int k = m; k < 2 * n; k += 2 * istep)
119: {
120: int j = k + istep;
121: double tempr = wr * data[j] - wi * data[j + 1];
122: double tempi = wi * data[j] + wr * data[j + 1];
123: data[j] = data[k] - tempr;
124: data[j + 1] = data[k + 1] - tempi;
125: data[k] = data[k] + tempr;
126: data[k + 1] = data[k + 1] + tempi;
127: }
128: }
129: mmax = istep;
130: }
131:
132: // copy out with optional scaling
133: if (!forward)
134: {
135: double scale = 1.0 / n;
136: for (int i = 0; i < 2 * n; ++i)
137: data[i] *= scale;
138: }
139: }
140:
141: /// <summary>
142: /// Compute the forward or inverse Fourier Transform of data, with
143: /// data containing real valued data only. The output is complex
144: /// valued after the first two entries, stored in alternating real
145: /// and imaginary parts. The first two returned entries are the real
146: /// parts of the first and last value from the conjugate symmetric
147: /// output, which are necessarily real. The length must be a power
148: /// of 2.
149: /// </summary>
150: /// <param name="data">The complex data stored as alternating real
151: /// and imaginary parts</param>
152: /// <param name="forward">true for a forward transform, false for
153: /// inverse transform</param>
154: public void RealFFT(double[] data, bool forward)
155: {
156: int n = data.Length; // # of real inputs, 1/2 the complex length
157: // checks n is a power of 2 in 2's complement format
158: if ((n & (n - 1)) != 0)
159: throw new ArgumentException(
160: "data length " + n + " in FFT is not a power of 2"
161: );
162:
163: double sign = -1;
164: if (forward)
165: { // do packed FFT. This can be changed to FFT to save memory
166: TableFFT(data, forward);
167: sign = 1;
168: }
169:
170: double theta = sign * 2 * Math.PI / n;
171: double wpr = Math.Cos(theta);
172: double wpi = Math.Sin(theta);
173: double wjr = wpr;
174: double wji = wpi;
175: for (int j = 1; j <= n / 4; ++j)
176: {
177: int k = n / 2 - j;
178: double tnr = data[2 * k];
179: double tni = data[2 * k + 1];
180: double tjr = data[2 * j];
181: double tji = data[2 * j + 1];
182:
183: double e = (tjr + tnr);
184: double f = (tji - tni);
185: double a = (tjr - tnr) * wji;
186: double d = (tji + tni) * wji;
187: double b = (tji + tni) * wjr;
188: double c = (tjr - tnr) * wjr;
189:
190: // compute entry y[j]
191: data[2 * j] = 0.5 * (e + sign * (a + b));
192: data[2 * j + 1] = 0.5 * (f - sign * (c - d));
193:
194: // compute entry y[k]
195: data[2 * k] = 0.5 * (e - sign * (a + b));
196: data[2 * k + 1] = 0.5 * (sign * (-c + d) - f);
197:
198: double temp = wjr;
199: // todo - allow more accurate version here? make option?
200: wjr = wjr * wpr - wji * wpi;
201: wji = temp * wpi + wji * wpr;
202: }
203:
204: if (forward)
205: {
206: // compute final y0 and y_{N/2}, store data[0], data[1]
207: double temp = data[0];
208: data[0] += data[1];
209: data[1] = temp - data[1];
210: }
211: else
212: {
213: double temp = data[0]; // unpack the y[j], then invert FFT
214: data[0] = 0.5 * (temp + data[1]);
215: data[1] = 0.5 * (temp - data[1]);
216: // do packed FFT. This can be changed to FFT to save memory
217: TableFFT(data, false);
218: }
219: }
220:
221: #region Internals
222:
223: /// <summary>
224: /// Call this with the size before using the TableFFT version
225: /// Fills in tables for speed. Done automatically in TableFFT
226: /// </summary>
227: /// <param name="size">The size of the FFT in samples</param>
228: void Initialize(int size)
229: {
230: // NOTE: if you do not use garbage collected languages
231: // like C# or Java be sure to free these correctly
232: cosTable = new List<double>();
233: sinTable = new List<double>();
234:
235: // forward pass
236: int n = size;
237: int mmax = 1;
238: while (n > mmax)
239: {
240: int istep = 2 * mmax;
241: double theta = Math.PI / mmax;
242: double wr = 1, wi = 0;
243: double wpi = Math.Sin(theta);
244: // compute in a slightly slower yet more accurate manner
245: double wpr = Math.Sin(theta / 2);
246: wpr = -2 * wpr * wpr;
247: for (int m = 0; m < istep; m += 2)
248: {
249: cosTable.Add(wr);
250: sinTable.Add(wi);
251: double t = wr;
252: wr = wr * wpr - wi * wpi + wr;
253: wi = wi * wpr + t * wpi + wi;
254: }
255: mmax = istep;
256: }
257: }
258:
259: /// <summary>
260: /// Swap data indices whenever index i has binary
261: /// digits reversed from index j, where data is
262: /// two doubles per index.
263: /// </summary>
264: /// <param name="data"></param>
265: /// <param name="n"></param>
266: void Reverse(double [] data, int n)
267: {
268: // bit reverse the indices. This is exercise 5 in section
269: // 7.2.1.1 of Knuth's TAOCP the idea is a binary counter
270: // in k and one with bits reversed in j
271: int j = 0, k = 0; // Knuth R1: initialize
272: int top = n / 2; // this is Knuth's 2^(n-1)
273: while (true)
274: {
275: // Knuth R2: swap - swap j+1 and k+2^(n-1), 2 entries each
276: double t = data[j + 2];
277: data[j + 2] = data[k + n];
278: data[k + n] = t;
279: t = data[j + 3];
280: data[j + 3] = data[k + n + 1];
281: data[k + n + 1] = t;
282: if (j > k)
283: { // swap two more
284: // j and k
285: t = data[j];
286: data[j] = data[k];
287: data[k] = t;
288: t = data[j + 1];
289: data[j + 1] = data[k + 1];
290: data[k + 1] = t;
291: // j + top + 1 and k+top + 1
292: t = data[j + n + 2];
293: data[j + n + 2] = data[k + n + 2];
294: data[k + n + 2] = t;
295: t = data[j + n + 3];
296: data[j + n + 3] = data[k + n + 3];
297: data[k + n + 3] = t;
298: }
299: // Knuth R3: advance k
300: k += 4;
301: if (k >= n)
302: break;
303: // Knuth R4: advance j
304: int h = top;
305: while (j >= h)
306: {
307: j -= h;
308: h /= 2;
309: }
310: j += h;
311: } // bit reverse loop
312: }
313:
314: /// <summary>
315: /// Precomputed sin/cos tables for speed
316: /// </summary>
317: List<double> cosTable;
318: List<double> sinTable;
319:
320: #endregion
321:
322: #region UnitTest
323: /// <summary>
324: /// Return true if unit tests pass
325: /// </summary>
326: /// <returns>true if and only if the tests all passed</returns>
327: public bool UnitTest()
328: {
329: // some tests of various lengths
330: double[] t4 = { 1, 1, 1, 1 }; // input
331: double[] a4r = { 4, 0, 0, 0 }; // real FFT
332: double[] a4c = { 2, 2, 0, 0 }; // complex FFT, ...
333: double[] t4a = { 1, 2, 3, 4 };
334: double[] a4ar = { 10, -2, -2, -2 };
335: double[] a4ac = { 4, 6, -2, -2 };
336: double[] t8 = { 0.100652, -0.442825, -0.457954, -0.00624455, 0.19978, -0.267328, -0.47192, -0.235878 };
337: double[] a8r = { -1.58172, 0.322834, -0.385598, 0.0522465, 1.23031, -0.468031, 0.187343, 0.0243147 };
338: double[] a8c = { -0.629442, -0.952276, -0.328761, -0.161531, 1.23031, -0.46803, 0.130505, -0.189463 };
339: double[] t32 = { -0.333615, 0.468917, 0.884538, 0.0276625, 0.979812, 0.91061, -0.175599, 0.1756, -0.695263, 0.557298, 0.112251, -0.285586, -0.73988, -0.0750604, -0.332421, 0.391004, 0.0588164, -0.18941, -0.416513, -0.596507, 0.659257, -0.654753, -0.472673, 0.875249, -0.00712734, -0.12367, -0.357211, -0.152413, 0.0130609, -0.0342799, 0.818388, 0.671986 };
340: double[] a32r = { 1.96247, -1.97083, 4.71435, 1.34203, 1.41278, 2.2209, -0.301542, 1.30462, 0.717877, -1.42063, -3.19595, -1.52441, -0.474644, -2.90705, 0.747585, 2.44391, -0.125698, -0.247344, -4.4128, -1.07521, -1.28254, 2.42047, -1.30217, -0.450559, -4.49676, -2.19137, 0.193633, 0.848902, 2.05478, -1.91513, 0.417439, 1.79843 };
341: double[] a32c = { -0.00417904, 1.96665, 1.44498, 2.18532, 1.46969, 1.82996, -1.0868, 0.620212, 1.23124, 0.951988, -2.48776, -1.88406, -0.412289, -2.73394, 0.564497, 2.93413, -0.125699, -0.247344, -4.22971, -0.584989, -1.3449, 2.59358, -2.01036, -0.810202, -5.01012, 0.181248, 0.97889, 0.164497, 1.99787, -2.30607, 3.68681, 2.64171 };
342:
343: double[][] tests = { t4, a4r, a4c, t4a, a4ar, a4ac, t8, a8r, a8c, t32, a32r, a32c };//, t32, a32 };
344:
345: bool ret = true;
346: for (int testIndex = 0; testIndex < tests.Length; testIndex += 3)
347: {
348: double[] test = tests[testIndex];
349: double[] answerReal = tests[testIndex + 1];
350: double[] answerComplex = tests[testIndex + 2];
351:
352: ret &= Test(RealFFT, test, answerReal);
353: ret &= Test(FFT, test, answerComplex);
354: ret &= Test(TableFFT, test, answerComplex);
355: }
356: return ret;
357: }
358:
359: /// <summary>
360: /// Test the given function on the given data and see if the result is the given answer.
361: /// </summary>
362: /// <returns>true if matches</returns>
363: bool Test(Action<double[], bool> FFTFunction, double[] test, double[] answer)
364: {
365: bool returnValue = true;
366: var copy = test.ToArray(); // make a copy
367: FFTFunction(copy, true); // forward transform
368: returnValue &= Compare(copy, answer); // check it
369: FFTFunction(copy, false); // backward transform
370: returnValue &= Compare(copy, test); // check it
371: return returnValue;
372: }
373:
374: /// <summary>
375: /// Compare two arrays of doubles for "equality"
376: /// up to a small tolerance
377: /// </summary>
378: /// <param name="arr1"></param>
379: /// <param name="arr2"></param>
380: /// <returns></returns>
381: static bool Compare(double[] arr1, double[] arr2)
382: {
383: if (arr1.Length != arr2.Length)
384: return false;
385: for (int i = 0; i < arr1.Length; ++i)
386: if ((Math.Abs(arr1[i] - arr2[i]) > 0.0001))
387: return false;
388: return true;
389: }
390:
391: #endregion
392: }
393: }
394: // end of file