fft.cpp
Go to the documentation of this file.
1 
6 
7 
8 // Copyright (C) 1991,2,3,4,8: R B Davies
9 
10 
11 #define WANT_MATH
12 // #define WANT_STREAM
13 
14 #include "include.h"
15 
16 #include "newmatap.h"
17 
18 // #include "newmatio.h"
19 
20 #ifdef use_namespace
21 namespace NEWMAT {
22 #endif
23 
24 #ifdef DO_REPORT
25 #define REPORT { static ExeCounter ExeCount(__LINE__,19); ++ExeCount; }
26 #else
27 #define REPORT {}
28 #endif
29 
30 static void cossin(int n, int d, Real& c, Real& s)
31 // calculate cos(twopi*n/d) and sin(twopi*n/d)
32 // minimise roundoff error
33 {
34  REPORT
35  long n4 = n * 4; int sector = (int)floor( (Real)n4 / (Real)d + 0.5 );
36  n4 -= sector * d;
37  if (sector < 0) { REPORT sector = 3 - (3 - sector) % 4; }
38  else { REPORT sector %= 4; }
39  Real ratio = 1.5707963267948966192 * (Real)n4 / (Real)d;
40 
41  switch (sector)
42  {
43  case 0: REPORT c = cos(ratio); s = sin(ratio); break;
44  case 1: REPORT c = -sin(ratio); s = cos(ratio); break;
45  case 2: REPORT c = -cos(ratio); s = -sin(ratio); break;
46  case 3: REPORT c = sin(ratio); s = -cos(ratio); break;
47  }
48 }
49 
51  ColumnVector& Y, int after, int now, int before)
52 {
53  REPORT
54  Tracer trace("FFT(step)");
55  // const Real twopi = 6.2831853071795864769;
56  const int gamma = after * before; const int delta = now * after;
57  // const Real angle = twopi / delta; Real temp;
58  // Real r_omega = cos(angle); Real i_omega = -sin(angle);
59  Real r_arg = 1.0; Real i_arg = 0.0;
60  Real* x = X.Store(); Real* y = Y.Store(); // pointers to array storage
61  const int m = A.Nrows() - gamma;
62 
63  for (int j = 0; j < now; j++)
64  {
65  Real* a = A.Store(); Real* b = B.Store(); // pointers to array storage
66  Real* x1 = x; Real* y1 = y; x += after; y += after;
67  for (int ia = 0; ia < after; ia++)
68  {
69  // generate sins & cosines explicitly rather than iteratively
70  // for more accuracy; but slower
71  cossin(-(j*after+ia), delta, r_arg, i_arg);
72 
73  Real* a1 = a++; Real* b1 = b++; Real* x2 = x1++; Real* y2 = y1++;
74  if (now==2)
75  {
76  REPORT int ib = before;
77  if (ib) for (;;)
78  {
79  REPORT
80  Real* a2 = m + a1; Real* b2 = m + b1; a1 += after; b1 += after;
81  Real r_value = *a2; Real i_value = *b2;
82  *x2 = r_value * r_arg - i_value * i_arg + *(a2-gamma);
83  *y2 = r_value * i_arg + i_value * r_arg + *(b2-gamma);
84  if (!(--ib)) break;
85  x2 += delta; y2 += delta;
86  }
87  }
88  else
89  {
90  REPORT int ib = before;
91  if (ib) for (;;)
92  {
93  REPORT
94  Real* a2 = m + a1; Real* b2 = m + b1; a1 += after; b1 += after;
95  Real r_value = *a2; Real i_value = *b2;
96  int in = now-1; while (in--)
97  {
98  // it should be possible to make this faster
99  // hand code for now = 2,3,4,5,8
100  // use symmetry to halve number of operations
101  a2 -= gamma; b2 -= gamma; Real temp = r_value;
102  r_value = r_value * r_arg - i_value * i_arg + *a2;
103  i_value = temp * i_arg + i_value * r_arg + *b2;
104  }
105  *x2 = r_value; *y2 = i_value;
106  if (!(--ib)) break;
107  x2 += delta; y2 += delta;
108  }
109  }
110 
111  // temp = r_arg;
112  // r_arg = r_arg * r_omega - i_arg * i_omega;
113  // i_arg = temp * i_omega + i_arg * r_omega;
114 
115  }
116  }
117 }
118 
119 
120 void FFTI(const ColumnVector& U, const ColumnVector& V,
121  ColumnVector& X, ColumnVector& Y)
122 {
123  // Inverse transform
124  Tracer trace("FFTI");
125  REPORT
126  FFT(U,-V,X,Y);
127  const Real n = X.Nrows(); X /= n; Y /= (-n);
128 }
129 
131 {
132  // Fourier transform of a real series
133  Tracer trace("RealFFT");
134  REPORT
135  const int n = U.Nrows(); // length of arrays
136  const int n2 = n / 2;
137  if (n != 2 * n2)
138  Throw(ProgramException("Vector length not multiple of 2", U));
139  ColumnVector A(n2), B(n2);
140  Real* a = A.Store(); Real* b = B.Store(); Real* u = U.Store(); int i = n2;
141  while (i--) { *a++ = *u++; *b++ = *u++; }
142  FFT(A,B,A,B);
143  int n21 = n2 + 1;
144  X.resize(n21); Y.resize(n21);
145  i = n2 - 1;
146  a = A.Store(); b = B.Store(); // first els of A and B
147  Real* an = a + i; Real* bn = b + i; // last els of A and B
148  Real* x = X.Store(); Real* y = Y.Store(); // first els of X and Y
149  Real* xn = x + n2; Real* yn = y + n2; // last els of X and Y
150 
151  *x++ = *a + *b; *y++ = 0.0; // first complex element
152  *xn-- = *a++ - *b++; *yn-- = 0.0; // last complex element
153 
154  int j = -1; i = n2/2;
155  while (i--)
156  {
157  Real c,s; cossin(j--,n,c,s);
158  Real am = *a - *an; Real ap = *a++ + *an--;
159  Real bm = *b - *bn; Real bp = *b++ + *bn--;
160  Real samcbp = s * am + c * bp; Real sbpcam = s * bp - c * am;
161  *x++ = 0.5 * ( ap + samcbp); *y++ = 0.5 * ( bm + sbpcam);
162  *xn-- = 0.5 * ( ap - samcbp); *yn-- = 0.5 * (-bm + sbpcam);
163  }
164 }
165 
166 void RealFFTI(const ColumnVector& A, const ColumnVector& B, ColumnVector& U)
167 {
168  // inverse of a Fourier transform of a real series
169  Tracer trace("RealFFTI");
170  REPORT
171  const int n21 = A.Nrows(); // length of arrays
172  if (n21 != B.Nrows() || n21 == 0)
173  Throw(ProgramException("Vector lengths unequal or zero", A, B));
174  const int n2 = n21 - 1; const int n = 2 * n2; int i = n2 - 1;
175 
176  ColumnVector X(n2), Y(n2);
177  Real* a = A.Store(); Real* b = B.Store(); // first els of A and B
178  Real* an = a + n2; Real* bn = b + n2; // last els of A and B
179  Real* x = X.Store(); Real* y = Y.Store(); // first els of X and Y
180  Real* xn = x + i; Real* yn = y + i; // last els of X and Y
181 
182  Real hn = 0.5 / n2;
183  *x++ = hn * (*a + *an); *y++ = - hn * (*a - *an);
184  a++; an--; b++; bn--;
185  int j = -1; i = n2/2;
186  while (i--)
187  {
188  Real c,s; cossin(j--,n,c,s);
189  Real am = *a - *an; Real ap = *a++ + *an--;
190  Real bm = *b - *bn; Real bp = *b++ + *bn--;
191  Real samcbp = s * am - c * bp; Real sbpcam = s * bp + c * am;
192  *x++ = hn * ( ap + samcbp); *y++ = - hn * ( bm + sbpcam);
193  *xn-- = hn * ( ap - samcbp); *yn-- = - hn * (-bm + sbpcam);
194  }
195  FFT(X,Y,X,Y); // have done inverting elsewhere
196  U.resize(n); i = n2;
197  x = X.Store(); y = Y.Store(); Real* u = U.Store();
198  while (i--) { *u++ = *x++; *u++ = - *y++; }
199 }
200 
201 void FFT(const ColumnVector& U, const ColumnVector& V,
202  ColumnVector& X, ColumnVector& Y)
203 {
204  // from Carl de Boor (1980), Siam J Sci Stat Comput, 1 173-8
205  // but first try Sande and Gentleman
206  Tracer trace("FFT");
207  REPORT
208  const int n = U.Nrows(); // length of arrays
209  if (n != V.Nrows() || n == 0)
210  Throw(ProgramException("Vector lengths unequal or zero", U, V));
211  if (n == 1) { REPORT X = U; Y = V; return; }
212 
213  // see if we can use the newfft routine
215  {
216  REPORT
217  X = U; Y = V;
218  if ( FFT_Controller::ar_1d_ft(n,X.Store(),Y.Store()) ) return;
219  }
220 
221  ColumnVector B = V;
222  ColumnVector A = U;
223  X.resize(n); Y.resize(n);
224  const int nextmx = 8;
225  int prime[8] = { 2,3,5,7,11,13,17,19 };
226  int after = 1; int before = n; int next = 0; bool inzee = true;
227  int now = 0; int b1; // initialised to keep gnu happy
228 
229  do
230  {
231  for (;;)
232  {
233  if (next < nextmx) { REPORT now = prime[next]; }
234  b1 = before / now; if (b1 * now == before) { REPORT break; }
235  next++; now += 2;
236  }
237  before = b1;
238 
239  if (inzee) { REPORT fftstep(A, B, X, Y, after, now, before); }
240  else { REPORT fftstep(X, Y, A, B, after, now, before); }
241 
242  inzee = !inzee; after *= now;
243  }
244  while (before != 1);
245 
246  if (inzee) { REPORT A.release(); X = A; B.release(); Y = B; }
247 }
248 
249 // Trigonometric transforms
250 // see Charles Van Loan (1992) "Computational frameworks for the fast
251 // Fourier transform" published by SIAM; section 4.4.
252 
253 void DCT_II(const ColumnVector& U, ColumnVector& V)
254 {
255  // Discrete cosine transform, type II, of a real series
256  Tracer trace("DCT_II");
257  REPORT
258  const int n = U.Nrows(); // length of arrays
259  const int n2 = n / 2; const int n4 = n * 4;
260  if (n != 2 * n2)
261  Throw(ProgramException("Vector length not multiple of 2", U));
262  ColumnVector A(n);
263  Real* a = A.Store(); Real* b = a + n; Real* u = U.Store();
264  int i = n2;
265  while (i--) { *a++ = *u++; *(--b) = *u++; }
266  ColumnVector X, Y;
267  RealFFT(A, X, Y); A.cleanup();
268  V.resize(n);
269  Real* x = X.Store(); Real* y = Y.Store();
270  Real* v = V.Store(); Real* w = v + n;
271  *v = *x;
272  int k = 0; i = n2;
273  while (i--)
274  {
275  Real c, s; cossin(++k, n4, c, s);
276  Real xi = *(++x); Real yi = *(++y);
277  *(++v) = xi * c + yi * s; *(--w) = xi * s - yi * c;
278  }
279 }
280 
282 {
283  // Inverse of discrete cosine transform, type II
284  Tracer trace("DCT_II_inverse");
285  REPORT
286  const int n = V.Nrows(); // length of array
287  const int n2 = n / 2; const int n4 = n * 4; const int n21 = n2 + 1;
288  if (n != 2 * n2)
289  Throw(ProgramException("Vector length not multiple of 2", V));
290  ColumnVector X(n21), Y(n21);
291  Real* x = X.Store(); Real* y = Y.Store();
292  Real* v = V.Store(); Real* w = v + n;
293  *x = *v; *y = 0.0;
294  int i = n2; int k = 0;
295  while (i--)
296  {
297  Real c, s; cossin(++k, n4, c, s);
298  Real vi = *(++v); Real wi = *(--w);
299  *(++x) = vi * c + wi * s; *(++y) = vi * s - wi * c;
300  }
301  ColumnVector A; RealFFTI(X, Y, A);
302  X.cleanup(); Y.cleanup(); U.resize(n);
303  Real* a = A.Store(); Real* b = a + n; Real* u = U.Store();
304  i = n2;
305  while (i--) { *u++ = *a++; *u++ = *(--b); }
306 }
307 
308 void DST_II(const ColumnVector& U, ColumnVector& V)
309 {
310  // Discrete sine transform, type II, of a real series
311  Tracer trace("DST_II");
312  REPORT
313  const int n = U.Nrows(); // length of arrays
314  const int n2 = n / 2; const int n4 = n * 4;
315  if (n != 2 * n2)
316  Throw(ProgramException("Vector length not multiple of 2", U));
317  ColumnVector A(n);
318  Real* a = A.Store(); Real* b = a + n; Real* u = U.Store();
319  int i = n2;
320  while (i--) { *a++ = *u++; *(--b) = -(*u++); }
321  ColumnVector X, Y;
322  RealFFT(A, X, Y); A.cleanup();
323  V.resize(n);
324  Real* x = X.Store(); Real* y = Y.Store();
325  Real* v = V.Store(); Real* w = v + n;
326  *(--w) = *x;
327  int k = 0; i = n2;
328  while (i--)
329  {
330  Real c, s; cossin(++k, n4, c, s);
331  Real xi = *(++x); Real yi = *(++y);
332  *v++ = xi * s - yi * c; *(--w) = xi * c + yi * s;
333  }
334 }
335 
337 {
338  // Inverse of discrete sine transform, type II
339  Tracer trace("DST_II_inverse");
340  REPORT
341  const int n = V.Nrows(); // length of array
342  const int n2 = n / 2; const int n4 = n * 4; const int n21 = n2 + 1;
343  if (n != 2 * n2)
344  Throw(ProgramException("Vector length not multiple of 2", V));
345  ColumnVector X(n21), Y(n21);
346  Real* x = X.Store(); Real* y = Y.Store();
347  Real* v = V.Store(); Real* w = v + n;
348  *x = *(--w); *y = 0.0;
349  int i = n2; int k = 0;
350  while (i--)
351  {
352  Real c, s; cossin(++k, n4, c, s);
353  Real vi = *v++; Real wi = *(--w);
354  *(++x) = vi * s + wi * c; *(++y) = - vi * c + wi * s;
355  }
356  ColumnVector A; RealFFTI(X, Y, A);
357  X.cleanup(); Y.cleanup(); U.resize(n);
358  Real* a = A.Store(); Real* b = a + n; Real* u = U.Store();
359  i = n2;
360  while (i--) { *u++ = *a++; *u++ = -(*(--b)); }
361 }
362 
364 {
365  // Inverse of discrete cosine transform, type I
366  Tracer trace("DCT_inverse");
367  REPORT
368  const int n = V.Nrows()-1; // length of transform
369  const int n2 = n / 2; const int n21 = n2 + 1;
370  if (n != 2 * n2)
371  Throw(ProgramException("Vector length not multiple of 2", V));
372  ColumnVector X(n21), Y(n21);
373  Real* x = X.Store(); Real* y = Y.Store(); Real* v = V.Store();
374  Real vi = *v++; *x++ = vi; *y++ = 0.0;
375  Real sum1 = vi / 2.0; Real sum2 = sum1; vi = *v++;
376  int i = n2-1;
377  while (i--)
378  {
379  Real vi2 = *v++; sum1 += vi2 + vi; sum2 += vi2 - vi;
380  *x++ = vi2; vi2 = *v++; *y++ = vi - vi2; vi = vi2;
381  }
382  sum1 += vi; sum2 -= vi;
383  vi = *v; *x = vi; *y = 0.0; vi /= 2.0; sum1 += vi; sum2 += vi;
384  ColumnVector A; RealFFTI(X, Y, A);
385  X.cleanup(); Y.cleanup(); U.resize(n+1);
386  Real* a = A.Store(); Real* b = a + n; Real* u = U.Store(); v = u + n;
387  i = n2; int k = 0; *u++ = sum1 / n2; *v-- = sum2 / n2;
388  while (i--)
389  {
390  Real s = sin(1.5707963267948966192 * (++k) / n2);
391  Real ai = *(++a); Real bi = *(--b);
392  Real bz = (ai - bi) / 4 / s; Real az = (ai + bi) / 2;
393  *u++ = az - bz; *v-- = az + bz;
394  }
395 }
396 
397 void DCT(const ColumnVector& U, ColumnVector& V)
398 {
399  // Discrete cosine transform, type I
400  Tracer trace("DCT");
401  REPORT
402  DCT_inverse(U, V);
403  V *= (V.Nrows()-1)/2;
404 }
405 
407 {
408  // Inverse of discrete sine transform, type I
409  Tracer trace("DST_inverse");
410  REPORT
411  const int n = V.Nrows()-1; // length of transform
412  const int n2 = n / 2; const int n21 = n2 + 1;
413  if (n != 2 * n2)
414  Throw(ProgramException("Vector length not multiple of 2", V));
415  ColumnVector X(n21), Y(n21);
416  Real* x = X.Store(); Real* y = Y.Store(); Real* v = V.Store();
417  Real vi = *(++v); *x++ = 2 * vi; *y++ = 0.0;
418  int i = n2-1;
419  while (i--) { *y++ = *(++v); Real vi2 = *(++v); *x++ = vi2 - vi; vi = vi2; }
420  *x = -2 * vi; *y = 0.0;
421  ColumnVector A; RealFFTI(X, Y, A);
422  X.cleanup(); Y.cleanup(); U.resize(n+1);
423  Real* a = A.Store(); Real* b = a + n; Real* u = U.Store(); v = u + n;
424  i = n2; int k = 0; *u++ = 0.0; *v-- = 0.0;
425  while (i--)
426  {
427  Real s = sin(1.5707963267948966192 * (++k) / n2);
428  Real ai = *(++a); Real bi = *(--b);
429  Real az = (ai + bi) / 4 / s; Real bz = (ai - bi) / 2;
430  *u++ = az - bz; *v-- = az + bz;
431  }
432 }
433 
434 void DST(const ColumnVector& U, ColumnVector& V)
435 {
436  // Discrete sine transform, type I
437  Tracer trace("DST");
438  REPORT
439  DST_inverse(U, V);
440  V *= (V.Nrows()-1)/2;
441 }
442 
443 // Two dimensional FFT
444 void FFT2(const Matrix& U, const Matrix& V, Matrix& X, Matrix& Y)
445 {
446  Tracer trace("FFT2");
447  REPORT
448  int m = U.Nrows(); int n = U.Ncols();
449  if (m != V.Nrows() || n != V.Ncols() || m == 0 || n == 0)
450  Throw(ProgramException("Matrix dimensions unequal or zero", U, V));
451  X = U; Y = V;
452  int i; ColumnVector CVR; ColumnVector CVI;
453  for (i = 1; i <= m; ++i)
454  {
455  FFT(X.Row(i).t(), Y.Row(i).t(), CVR, CVI);
456  X.Row(i) = CVR.t(); Y.Row(i) = CVI.t();
457  }
458  for (i = 1; i <= n; ++i)
459  {
460  FFT(X.Column(i), Y.Column(i), CVR, CVI);
461  X.Column(i) = CVR; Y.Column(i) = CVI;
462  }
463 }
464 
465 void FFT2I(const Matrix& U, const Matrix& V, Matrix& X, Matrix& Y)
466 {
467  // Inverse transform
468  Tracer trace("FFT2I");
469  REPORT
470  FFT2(U,-V,X,Y);
471  const Real n = X.Nrows() * X.Ncols(); X /= n; Y /= (-n);
472 }
473 
474 
475 #ifdef use_namespace
476 }
477 #endif
478 
479 
void DST_inverse(const ColumnVector &V, ColumnVector &U)
Definition: fft.cpp:406
void DST_II(const ColumnVector &U, ColumnVector &V)
Definition: fft.cpp:308
void RealFFT(const ColumnVector &U, ColumnVector &X, ColumnVector &Y)
Definition: fft.cpp:130
void cleanup()
Definition: newmat4.cpp:1125
Miscellaneous exception (details in character string).
Definition: newmat.h:1947
double Real
Definition: include.h:307
int Nrows() const
Definition: newmat.h:494
void DCT_II(const ColumnVector &U, ColumnVector &V)
Definition: fft.cpp:253
void FFT(const ColumnVector &U, const ColumnVector &V, ColumnVector &X, ColumnVector &Y)
Definition: fft.cpp:201
void FFT2(const Matrix &U, const Matrix &V, Matrix &X, Matrix &Y)
Definition: fft.cpp:444
void DCT_II_inverse(const ColumnVector &V, ColumnVector &U)
Definition: fft.cpp:281
void FFT2I(const Matrix &U, const Matrix &V, Matrix &X, Matrix &Y)
Definition: fft.cpp:465
void resize(int)
Definition: newmat4.cpp:318
GetSubMatrix Column(int f) const
Definition: newmat.h:2152
FloatVector * d
void FFTI(const ColumnVector &U, const ColumnVector &V, ColumnVector &X, ColumnVector &Y)
Definition: fft.cpp:120
TransposedMatrix t() const
Definition: newmat6.cpp:320
Real * Store() const
Definition: newmat.h:497
void DST(const ColumnVector &U, ColumnVector &V)
Definition: fft.cpp:434
static void fftstep(ColumnVector &A, ColumnVector &B, ColumnVector &X, ColumnVector &Y, int after, int now, int before)
Definition: fft.cpp:50
#define Throw(E)
Definition: myexcept.h:191
The usual rectangular matrix.
Definition: newmat.h:625
FloatVector FloatVector * a
int Ncols() const
Definition: newmat.h:495
Real trace(const BaseMatrix &B)
Definition: newmat.h:2099
static void cossin(int n, int d, Real &c, Real &s)
Definition: fft.cpp:30
static bool ar_1d_ft(int PTS, Real *X, Real *Y)
Definition: newfft.cpp:150
void RealFFTI(const ColumnVector &A, const ColumnVector &B, ColumnVector &U)
Definition: fft.cpp:166
static bool OnlyOldFFT
Definition: newmatap.h:145
#define REPORT
Definition: fft.cpp:27
GetSubMatrix Row(int f) const
Definition: newmat.h:2150
void DCT(const ColumnVector &U, ColumnVector &V)
Definition: fft.cpp:397
void DCT_inverse(const ColumnVector &V, ColumnVector &U)
Definition: fft.cpp:363
static bool CanFactor(int PTS)
Definition: newfft.cpp:992
Column vector.
Definition: newmat.h:1008
void DST_II_inverse(const ColumnVector &V, ColumnVector &U)
Definition: fft.cpp:336
void release()
Definition: newmat.h:517


kni
Author(s): Martin Günther
autogenerated on Fri Jan 3 2020 04:01:16