From eb976fa09171036cbaeaabf920708b2d39c49acc Mon Sep 17 00:00:00 2001 From: Miller Puckette Date: Sat, 3 Jun 2006 18:22:09 +0000 Subject: Removing renamed files svn path=/trunk/; revision=5163 --- pd/src/d_fft_mayer.c | 423 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 423 insertions(+) create mode 100644 pd/src/d_fft_mayer.c (limited to 'pd/src/d_fft_mayer.c') diff --git a/pd/src/d_fft_mayer.c b/pd/src/d_fft_mayer.c new file mode 100644 index 00000000..860b3120 --- /dev/null +++ b/pd/src/d_fft_mayer.c @@ -0,0 +1,423 @@ +/* +** FFT and FHT routines +** Copyright 1988, 1993; Ron Mayer +** +** mayer_fht(fz,n); +** Does a hartley transform of "n" points in the array "fz". +** mayer_fft(n,real,imag) +** Does a fourier transform of "n" points of the "real" and +** "imag" arrays. +** mayer_ifft(n,real,imag) +** Does an inverse fourier transform of "n" points of the "real" +** and "imag" arrays. +** mayer_realfft(n,real) +** Does a real-valued fourier transform of "n" points of the +** "real" array. The real part of the transform ends +** up in the first half of the array and the imaginary part of the +** transform ends up in the second half of the array. +** mayer_realifft(n,real) +** The inverse of the realfft() routine above. +** +** +** NOTE: This routine uses at least 2 patented algorithms, and may be +** under the restrictions of a bunch of different organizations. +** Although I wrote it completely myself, it is kind of a derivative +** of a routine I once authored and released under the GPL, so it +** may fall under the free software foundation's restrictions; +** it was worked on as a Stanford Univ project, so they claim +** some rights to it; it was further optimized at work here, so +** I think this company claims parts of it. The patents are +** held by R. Bracewell (the FHT algorithm) and O. Buneman (the +** trig generator), both at Stanford Univ. +** If it were up to me, I'd say go do whatever you want with it; +** but it would be polite to give credit to the following people +** if you use this anywhere: +** Euler - probable inventor of the fourier transform. +** Gauss - probable inventor of the FFT. +** Hartley - probable inventor of the hartley transform. +** Buneman - for a really cool trig generator +** Mayer(me) - for authoring this particular version and +** including all the optimizations in one package. +** Thanks, +** Ron Mayer; mayer@acuson.com +** +*/ + +/* This is a slightly modified version of Mayer's contribution; write +* msp@ucsd.edu for the original code. Kudos to Mayer for a fine piece +* of work. -msp +*/ + +/* These pragmas are only used for MSVC, not MinGW or Cygwin */ +#ifdef _MSC_VER +#pragma warning( disable : 4305 ) /* uncast const double to float */ +#pragma warning( disable : 4244 ) /* uncast double to float */ +#pragma warning( disable : 4101 ) /* unused local variables */ +#endif + +/* the following is needed only to declare pd_fft() as exportable in MSW */ +#include "m_pd.h" + +#define REAL float +#define GOOD_TRIG + +#ifdef GOOD_TRIG +#else +#define FAST_TRIG +#endif + +#if defined(GOOD_TRIG) +#define FHT_SWAP(a,b,t) {(t)=(a);(a)=(b);(b)=(t);} +#define TRIG_VARS \ + int t_lam=0; +#define TRIG_INIT(k,c,s) \ + { \ + int i; \ + for (i=2 ; i<=k ; i++) \ + {coswrk[i]=costab[i];sinwrk[i]=sintab[i];} \ + t_lam = 0; \ + c = 1; \ + s = 0; \ + } +#define TRIG_NEXT(k,c,s) \ + { \ + int i,j; \ + (t_lam)++; \ + for (i=0 ; !((1<1) \ + { \ + for (j=k-i+2 ; (1<>1; (!((k2^=k)&k)); k>>=1); + if (k1>k2) + { + aa=fz[k1];fz[k1]=fz[k2];fz[k2]=aa; + } + } + for ( k=0 ; (1<> 1; + fi = fz; + gi = fi + kx; + fn = fz + n; + do + { + REAL g0,f0,f1,g1,f2,g2,f3,g3; + f1 = fi[0 ] - fi[k1]; + f0 = fi[0 ] + fi[k1]; + f3 = fi[k2] - fi[k3]; + f2 = fi[k2] + fi[k3]; + fi[k2] = f0 - f2; + fi[0 ] = f0 + f2; + fi[k3] = f1 - f3; + fi[k1] = f1 + f3; + g1 = gi[0 ] - gi[k1]; + g0 = gi[0 ] + gi[k1]; + g3 = SQRT2 * gi[k3]; + g2 = SQRT2 * gi[k2]; + gi[k2] = g0 - g2; + gi[0 ] = g0 + g2; + gi[k3] = g1 - g3; + gi[k1] = g1 + g3; + gi += k4; + fi += k4; + } while (fi