Jamoma API  0.6.0.a19
mayer_fft.cpp
1 // This is the FFT from PureData by Miller Puckette
2 // http://crca.ucsd.edu/~msp/software.html
3 
4 /*
5 ** FFT and FHT routines
6 ** Copyright 1988, 1993; Ron Mayer
7 **
8 ** mayer_fht(fz,n);
9 ** Does a hartley transform of "n" points in the array "fz".
10 ** mayer_fft(n,real,imag)
11 ** Does a fourier transform of "n" points of the "real" and
12 ** "imag" arrays.
13 ** mayer_ifft(n,real,imag)
14 ** Does an inverse fourier transform of "n" points of the "real"
15 ** and "imag" arrays.
16 ** mayer_realfft(n,real)
17 ** Does a real-valued fourier transform of "n" points of the
18 ** "real" array. The real part of the transform ends
19 ** up in the first half of the array and the imaginary part of the
20 ** transform ends up in the second half of the array.
21 ** mayer_realifft(n,real)
22 ** The inverse of the realfft() routine above.
23 **
24 **
25 ** NOTE: This routine uses at least 2 patented algorithms, and may be
26 ** under the restrictions of a bunch of different organizations.
27 ** Although I wrote it completely myself, it is kind of a derivative
28 ** of a routine I once authored and released under the GPL, so it
29 ** may fall under the free software foundation's restrictions;
30 ** it was worked on as a Stanford Univ project, so they claim
31 ** some rights to it; it was further optimized at work here, so
32 ** I think this company claims parts of it. The patents are
33 ** held by R. Bracewell (the FHT algorithm) and O. Buneman (the
34 ** trig generator), both at Stanford Univ.
35 ** If it were up to me, I'd say go do whatever you want with it;
36 ** but it would be polite to give credit to the following people
37 ** if you use this anywhere:
38 ** Euler - probable inventor of the fourier transform.
39 ** Gauss - probable inventor of the FFT.
40 ** Hartley - probable inventor of the hartley transform.
41 ** Buneman - for a really cool trig generator
42 ** Mayer(me) - for authoring this particular version and
43 ** including all the optimizations in one package.
44 ** Thanks,
45 ** Ron Mayer; mayer@acuson.com
46 **
47 */
48 
49 /* This is a slightly modified version of Mayer's contribution; write
50 * msp@ucsd.edu for the original code. Kudos to Mayer for a fine piece
51 * of work. -msp
52 */
53 
54 #include "mayer_fft.h"
55 
56 #define REAL double
57 #define GOOD_TRIG
58 
59 #ifdef GOOD_TRIG
60 #else
61 #define FAST_TRIG
62 #endif
63 
64 #if defined(GOOD_TRIG)
65 #define FHT_SWAP(a,b,t) {(t)=(a);(a)=(b);(b)=(t);}
66 #define TRIG_VARS \
67  int t_lam=0;
68 #define TRIG_INIT(k,c,s) \
69  { \
70  int i; \
71  for (i=2 ; i<=k ; i++) \
72  {coswrk[i]=costab[i];sinwrk[i]=sintab[i];} \
73  t_lam = 0; \
74  c = 1; \
75  s = 0; \
76  }
77 #define TRIG_NEXT(k,c,s) \
78  { \
79  int i,j; \
80  (t_lam)++; \
81  for (i=0 ; !((1<<i)&t_lam) ; i++); \
82  i = k-i; \
83  s = sinwrk[i]; \
84  c = coswrk[i]; \
85  if (i>1) \
86  { \
87  for (j=k-i+2 ; (1<<j)&t_lam ; j++); \
88  j = k - j; \
89  sinwrk[i] = halsec[i] * (sinwrk[i-1] + sinwrk[j]); \
90  coswrk[i] = halsec[i] * (coswrk[i-1] + coswrk[j]); \
91  } \
92  }
93 #define TRIG_RESET(k,c,s)
94 #endif
95 
96 #if defined(FAST_TRIG)
97 #define TRIG_VARS \
98  REAL t_c,t_s;
99 #define TRIG_INIT(k,c,s) \
100  { \
101  t_c = costab[k]; \
102  t_s = sintab[k]; \
103  c = 1; \
104  s = 0; \
105  }
106 #define TRIG_NEXT(k,c,s) \
107  { \
108  REAL t = c; \
109  c = t*t_c - s*t_s; \
110  s = t*t_s + s*t_c; \
111  }
112 #define TRIG_RESET(k,c,s)
113 #endif
114 
115 static REAL halsec[20]=
116  {
117  0,
118  0,
119  .54119610014619698439972320536638942006107206337801,
120  .50979557910415916894193980398784391368261849190893,
121  .50241928618815570551167011928012092247859337193963,
122  .50060299823519630134550410676638239611758632599591,
123  .50015063602065098821477101271097658495974913010340,
124  .50003765191554772296778139077905492847503165398345,
125  .50000941253588775676512870469186533538523133757983,
126  .50000235310628608051401267171204408939326297376426,
127  .50000058827484117879868526730916804925780637276181,
128  .50000014706860214875463798283871198206179118093251,
129  .50000003676714377807315864400643020315103490883972,
130  .50000000919178552207366560348853455333939112569380,
131  .50000000229794635411562887767906868558991922348920,
132  .50000000057448658687873302235147272458812263401372
133  };
134 static REAL costab[20]=
135  {
136  .00000000000000000000000000000000000000000000000000,
137  .70710678118654752440084436210484903928483593768847,
138  .92387953251128675612818318939678828682241662586364,
139  .98078528040323044912618223613423903697393373089333,
140  .99518472667219688624483695310947992157547486872985,
141  .99879545620517239271477160475910069444320361470461,
142  .99969881869620422011576564966617219685006108125772,
143  .99992470183914454092164649119638322435060646880221,
144  .99998117528260114265699043772856771617391725094433,
145  .99999529380957617151158012570011989955298763362218,
146  .99999882345170190992902571017152601904826792288976,
147  .99999970586288221916022821773876567711626389934930,
148  .99999992646571785114473148070738785694820115568892,
149  .99999998161642929380834691540290971450507605124278,
150  .99999999540410731289097193313960614895889430318945,
151  .99999999885102682756267330779455410840053741619428
152  };
153 static REAL sintab[20]=
154  {
155  1.0000000000000000000000000000000000000000000000000,
156  .70710678118654752440084436210484903928483593768846,
157  .38268343236508977172845998403039886676134456248561,
158  .19509032201612826784828486847702224092769161775195,
159  .09801714032956060199419556388864184586113667316749,
160  .04906767432741801425495497694268265831474536302574,
161  .02454122852291228803173452945928292506546611923944,
162  .01227153828571992607940826195100321214037231959176,
163  .00613588464915447535964023459037258091705788631738,
164  .00306795676296597627014536549091984251894461021344,
165  .00153398018628476561230369715026407907995486457522,
166  .00076699031874270452693856835794857664314091945205,
167  .00038349518757139558907246168118138126339502603495,
168  .00019174759731070330743990956198900093346887403385,
169  .00009587379909597734587051721097647635118706561284,
170  .00004793689960306688454900399049465887274686668768
171  };
172 static REAL coswrk[20]=
173  {
174  .00000000000000000000000000000000000000000000000000,
175  .70710678118654752440084436210484903928483593768847,
176  .92387953251128675612818318939678828682241662586364,
177  .98078528040323044912618223613423903697393373089333,
178  .99518472667219688624483695310947992157547486872985,
179  .99879545620517239271477160475910069444320361470461,
180  .99969881869620422011576564966617219685006108125772,
181  .99992470183914454092164649119638322435060646880221,
182  .99998117528260114265699043772856771617391725094433,
183  .99999529380957617151158012570011989955298763362218,
184  .99999882345170190992902571017152601904826792288976,
185  .99999970586288221916022821773876567711626389934930,
186  .99999992646571785114473148070738785694820115568892,
187  .99999998161642929380834691540290971450507605124278,
188  .99999999540410731289097193313960614895889430318945,
189  .99999999885102682756267330779455410840053741619428
190  };
191 static REAL sinwrk[20]=
192  {
193  1.0000000000000000000000000000000000000000000000000,
194  .70710678118654752440084436210484903928483593768846,
195  .38268343236508977172845998403039886676134456248561,
196  .19509032201612826784828486847702224092769161775195,
197  .09801714032956060199419556388864184586113667316749,
198  .04906767432741801425495497694268265831474536302574,
199  .02454122852291228803173452945928292506546611923944,
200  .01227153828571992607940826195100321214037231959176,
201  .00613588464915447535964023459037258091705788631738,
202  .00306795676296597627014536549091984251894461021344,
203  .00153398018628476561230369715026407907995486457522,
204  .00076699031874270452693856835794857664314091945205,
205  .00038349518757139558907246168118138126339502603495,
206  .00019174759731070330743990956198900093346887403385,
207  .00009587379909597734587051721097647635118706561284,
208  .00004793689960306688454900399049465887274686668768
209  };
210 
211 
212 #define SQRT2_2 0.70710678118654752440084436210484
213 #define SQRT2 2*0.70710678118654752440084436210484
214 
215 void mayer_fht(REAL *fz, int n)
216 {
217 /* REAL a,b;
218 REAL c1,s1,s2,c2,s3,c3,s4,c4;
219  REAL f0,g0,f1,g1,f2,g2,f3,g3; */
220  int k,k1,k2,k3,k4,kx;
221  REAL *fi,*fn,*gi;
222  TRIG_VARS;
223 
224  for (k1=1,k2=0;k1<n;k1++)
225  {
226  REAL aa;
227  for (k=n>>1; (!((k2^=k)&k)); k>>=1);
228  if (k1>k2)
229  {
230  aa=fz[k1];fz[k1]=fz[k2];fz[k2]=aa;
231  }
232  }
233  for ( k=0 ; (1<<k)<n ; k++ );
234  k &= 1;
235  if (k==0)
236  {
237  for (fi=fz,fn=fz+n;fi<fn;fi+=4)
238  {
239  REAL f0,f1,f2,f3;
240  f1 = fi[0 ]-fi[1 ];
241  f0 = fi[0 ]+fi[1 ];
242  f3 = fi[2 ]-fi[3 ];
243  f2 = fi[2 ]+fi[3 ];
244  fi[2 ] = (f0-f2);
245  fi[0 ] = (f0+f2);
246  fi[3 ] = (f1-f3);
247  fi[1 ] = (f1+f3);
248  }
249  }
250  else
251  {
252  for (fi=fz,fn=fz+n,gi=fi+1;fi<fn;fi+=8,gi+=8)
253  {
254  REAL bs1,bc1,bs2,bc2,bs3,bc3,bs4,bc4,
255  bg0,bf0,bf1,bg1,bf2,bg2,bf3,bg3;
256  bc1 = fi[0 ] - gi[0 ];
257  bs1 = fi[0 ] + gi[0 ];
258  bc2 = fi[2 ] - gi[2 ];
259  bs2 = fi[2 ] + gi[2 ];
260  bc3 = fi[4 ] - gi[4 ];
261  bs3 = fi[4 ] + gi[4 ];
262  bc4 = fi[6 ] - gi[6 ];
263  bs4 = fi[6 ] + gi[6 ];
264  bf1 = (bs1 - bs2);
265  bf0 = (bs1 + bs2);
266  bg1 = (bc1 - bc2);
267  bg0 = (bc1 + bc2);
268  bf3 = (bs3 - bs4);
269  bf2 = (bs3 + bs4);
270  bg3 = SQRT2*bc4;
271  bg2 = SQRT2*bc3;
272  fi[4 ] = bf0 - bf2;
273  fi[0 ] = bf0 + bf2;
274  fi[6 ] = bf1 - bf3;
275  fi[2 ] = bf1 + bf3;
276  gi[4 ] = bg0 - bg2;
277  gi[0 ] = bg0 + bg2;
278  gi[6 ] = bg1 - bg3;
279  gi[2 ] = bg1 + bg3;
280  }
281  }
282  if (n<16) return;
283 
284  do
285  {
286  REAL s1,c1;
287  int ii;
288  k += 2;
289  k1 = 1 << k;
290  k2 = k1 << 1;
291  k4 = k2 << 1;
292  k3 = k2 + k1;
293  kx = k1 >> 1;
294  fi = fz;
295  gi = fi + kx;
296  fn = fz + n;
297  do
298  {
299  REAL g0,f0,f1,g1,f2,g2,f3,g3;
300  f1 = fi[0 ] - fi[k1];
301  f0 = fi[0 ] + fi[k1];
302  f3 = fi[k2] - fi[k3];
303  f2 = fi[k2] + fi[k3];
304  fi[k2] = f0 - f2;
305  fi[0 ] = f0 + f2;
306  fi[k3] = f1 - f3;
307  fi[k1] = f1 + f3;
308  g1 = gi[0 ] - gi[k1];
309  g0 = gi[0 ] + gi[k1];
310  g3 = SQRT2 * gi[k3];
311  g2 = SQRT2 * gi[k2];
312  gi[k2] = g0 - g2;
313  gi[0 ] = g0 + g2;
314  gi[k3] = g1 - g3;
315  gi[k1] = g1 + g3;
316  gi += k4;
317  fi += k4;
318  } while (fi<fn);
319  TRIG_INIT(k,c1,s1);
320  for (ii=1;ii<kx;ii++)
321  {
322  REAL c2,s2;
323  TRIG_NEXT(k,c1,s1);
324  c2 = c1*c1 - s1*s1;
325  s2 = 2*(c1*s1);
326  fn = fz + n;
327  fi = fz +ii;
328  gi = fz +k1-ii;
329  do
330  {
331  REAL a,b,g0,f0,f1,g1,f2,g2,f3,g3;
332  b = s2*fi[k1] - c2*gi[k1];
333  a = c2*fi[k1] + s2*gi[k1];
334  f1 = fi[0 ] - a;
335  f0 = fi[0 ] + a;
336  g1 = gi[0 ] - b;
337  g0 = gi[0 ] + b;
338  b = s2*fi[k3] - c2*gi[k3];
339  a = c2*fi[k3] + s2*gi[k3];
340  f3 = fi[k2] - a;
341  f2 = fi[k2] + a;
342  g3 = gi[k2] - b;
343  g2 = gi[k2] + b;
344  b = s1*f2 - c1*g3;
345  a = c1*f2 + s1*g3;
346  fi[k2] = f0 - a;
347  fi[0 ] = f0 + a;
348  gi[k3] = g1 - b;
349  gi[k1] = g1 + b;
350  b = c1*g2 - s1*f3;
351  a = s1*g2 + c1*f3;
352  gi[k2] = g0 - a;
353  gi[0 ] = g0 + a;
354  fi[k3] = f1 - b;
355  fi[k1] = f1 + b;
356  gi += k4;
357  fi += k4;
358  } while (fi<fn);
359  }
360  TRIG_RESET(k,c1,s1);
361  } while (k4<n);
362 }
363 
364 void mayer_fft(int n, REAL *real, REAL *imag)
365 {
366  REAL a,b,c,d;
367  REAL q,r,s,t;
368  int i,j,k;
369  for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
370  a = real[i]; b = real[j]; q=a+b; r=a-b;
371  c = imag[i]; d = imag[j]; s=c+d; t=c-d;
372  real[i] = (q+t)*.5; real[j] = (q-t)*.5;
373  imag[i] = (s-r)*.5; imag[j] = (s+r)*.5;
374  }
375  mayer_fht(real,n);
376  mayer_fht(imag,n);
377 }
378 
379 void mayer_ifft(int n, REAL *real, REAL *imag)
380 {
381  REAL a,b,c,d;
382  REAL q,r,s,t;
383  int i,j,k;
384  mayer_fht(real,n);
385  mayer_fht(imag,n);
386  for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
387  a = real[i]; b = real[j]; q=a+b; r=a-b;
388  c = imag[i]; d = imag[j]; s=c+d; t=c-d;
389  imag[i] = (s+r)*0.5; imag[j] = (s-r)*0.5;
390  real[i] = (q-t)*0.5; real[j] = (q+t)*0.5;
391  }
392 }
393 
394 void mayer_realfft(int n, REAL *real)
395 {
396  REAL a,b;
397  int i,j,k;
398 
399  mayer_fht(real,n);
400  for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
401  a = real[i];
402  b = real[j];
403  real[j] = (a-b)*0.5;
404  real[i] = (a+b)*0.5;
405  }
406 }
407 
408 void mayer_realifft(int n, REAL *real)
409 {
410  REAL a,b;
411  int i,j,k;
412 
413  for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
414  a = real[i];
415  b = real[j];
416  real[j] = (a-b);
417  real[i] = (a+b);
418  }
419  mayer_fht(real,n);
420 }