diff options
Diffstat (limited to 'pd/src/d_mayer_fft.c')
| -rw-r--r-- | pd/src/d_mayer_fft.c | 423 |
1 files changed, 0 insertions, 423 deletions
diff --git a/pd/src/d_mayer_fft.c b/pd/src/d_mayer_fft.c deleted file mode 100644 index 860b3120..00000000 --- a/pd/src/d_mayer_fft.c +++ /dev/null @@ -1,423 +0,0 @@ -/* -** 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 <hans@at.or.at> */ -#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<<i)&t_lam) ; i++); \ - i = k-i; \ - s = sinwrk[i]; \ - c = coswrk[i]; \ - if (i>1) \ - { \ - for (j=k-i+2 ; (1<<j)&t_lam ; j++); \ - j = k - j; \ - sinwrk[i] = halsec[i] * (sinwrk[i-1] + sinwrk[j]); \ - coswrk[i] = halsec[i] * (coswrk[i-1] + coswrk[j]); \ - } \ - } -#define TRIG_RESET(k,c,s) -#endif - -#if defined(FAST_TRIG) -#define TRIG_VARS \ - REAL t_c,t_s; -#define TRIG_INIT(k,c,s) \ - { \ - t_c = costab[k]; \ - t_s = sintab[k]; \ - c = 1; \ - s = 0; \ - } -#define TRIG_NEXT(k,c,s) \ - { \ - REAL t = c; \ - c = t*t_c - s*t_s; \ - s = t*t_s + s*t_c; \ - } -#define TRIG_RESET(k,c,s) -#endif - -static REAL halsec[20]= - { - 0, - 0, - .54119610014619698439972320536638942006107206337801, - .50979557910415916894193980398784391368261849190893, - .50241928618815570551167011928012092247859337193963, - .50060299823519630134550410676638239611758632599591, - .50015063602065098821477101271097658495974913010340, - .50003765191554772296778139077905492847503165398345, - .50000941253588775676512870469186533538523133757983, - .50000235310628608051401267171204408939326297376426, - .50000058827484117879868526730916804925780637276181, - .50000014706860214875463798283871198206179118093251, - .50000003676714377807315864400643020315103490883972, - .50000000919178552207366560348853455333939112569380, - .50000000229794635411562887767906868558991922348920, - .50000000057448658687873302235147272458812263401372 - }; -static REAL costab[20]= - { - .00000000000000000000000000000000000000000000000000, - .70710678118654752440084436210484903928483593768847, - .92387953251128675612818318939678828682241662586364, - .98078528040323044912618223613423903697393373089333, - .99518472667219688624483695310947992157547486872985, - .99879545620517239271477160475910069444320361470461, - .99969881869620422011576564966617219685006108125772, - .99992470183914454092164649119638322435060646880221, - .99998117528260114265699043772856771617391725094433, - .99999529380957617151158012570011989955298763362218, - .99999882345170190992902571017152601904826792288976, - .99999970586288221916022821773876567711626389934930, - .99999992646571785114473148070738785694820115568892, - .99999998161642929380834691540290971450507605124278, - .99999999540410731289097193313960614895889430318945, - .99999999885102682756267330779455410840053741619428 - }; -static REAL sintab[20]= - { - 1.0000000000000000000000000000000000000000000000000, - .70710678118654752440084436210484903928483593768846, - .38268343236508977172845998403039886676134456248561, - .19509032201612826784828486847702224092769161775195, - .09801714032956060199419556388864184586113667316749, - .04906767432741801425495497694268265831474536302574, - .02454122852291228803173452945928292506546611923944, - .01227153828571992607940826195100321214037231959176, - .00613588464915447535964023459037258091705788631738, - .00306795676296597627014536549091984251894461021344, - .00153398018628476561230369715026407907995486457522, - .00076699031874270452693856835794857664314091945205, - .00038349518757139558907246168118138126339502603495, - .00019174759731070330743990956198900093346887403385, - .00009587379909597734587051721097647635118706561284, - .00004793689960306688454900399049465887274686668768 - }; -static REAL coswrk[20]= - { - .00000000000000000000000000000000000000000000000000, - .70710678118654752440084436210484903928483593768847, - .92387953251128675612818318939678828682241662586364, - .98078528040323044912618223613423903697393373089333, - .99518472667219688624483695310947992157547486872985, - .99879545620517239271477160475910069444320361470461, - .99969881869620422011576564966617219685006108125772, - .99992470183914454092164649119638322435060646880221, - .99998117528260114265699043772856771617391725094433, - .99999529380957617151158012570011989955298763362218, - .99999882345170190992902571017152601904826792288976, - .99999970586288221916022821773876567711626389934930, - .99999992646571785114473148070738785694820115568892, - .99999998161642929380834691540290971450507605124278, - .99999999540410731289097193313960614895889430318945, - .99999999885102682756267330779455410840053741619428 - }; -static REAL sinwrk[20]= - { - 1.0000000000000000000000000000000000000000000000000, - .70710678118654752440084436210484903928483593768846, - .38268343236508977172845998403039886676134456248561, - .19509032201612826784828486847702224092769161775195, - .09801714032956060199419556388864184586113667316749, - .04906767432741801425495497694268265831474536302574, - .02454122852291228803173452945928292506546611923944, - .01227153828571992607940826195100321214037231959176, - .00613588464915447535964023459037258091705788631738, - .00306795676296597627014536549091984251894461021344, - .00153398018628476561230369715026407907995486457522, - .00076699031874270452693856835794857664314091945205, - .00038349518757139558907246168118138126339502603495, - .00019174759731070330743990956198900093346887403385, - .00009587379909597734587051721097647635118706561284, - .00004793689960306688454900399049465887274686668768 - }; - - -#define SQRT2_2 0.70710678118654752440084436210484 -#define SQRT2 2*0.70710678118654752440084436210484 - -void mayer_fht(REAL *fz, int n) -{ -/* REAL a,b; -REAL c1,s1,s2,c2,s3,c3,s4,c4; - REAL f0,g0,f1,g1,f2,g2,f3,g3; */ - int k,k1,k2,k3,k4,kx; - REAL *fi,*fn,*gi; - TRIG_VARS; - - for (k1=1,k2=0;k1<n;k1++) - { - REAL aa; - for (k=n>>1; (!((k2^=k)&k)); k>>=1); - if (k1>k2) - { - aa=fz[k1];fz[k1]=fz[k2];fz[k2]=aa; - } - } - for ( k=0 ; (1<<k)<n ; k++ ); - k &= 1; - if (k==0) - { - for (fi=fz,fn=fz+n;fi<fn;fi+=4) - { - REAL f0,f1,f2,f3; - f1 = fi[0 ]-fi[1 ]; - f0 = fi[0 ]+fi[1 ]; - f3 = fi[2 ]-fi[3 ]; - f2 = fi[2 ]+fi[3 ]; - fi[2 ] = (f0-f2); - fi[0 ] = (f0+f2); - fi[3 ] = (f1-f3); - fi[1 ] = (f1+f3); - } - } - else - { - for (fi=fz,fn=fz+n,gi=fi+1;fi<fn;fi+=8,gi+=8) - { - REAL bs1,bc1,bs2,bc2,bs3,bc3,bs4,bc4, - bg0,bf0,bf1,bg1,bf2,bg2,bf3,bg3; - bc1 = fi[0 ] - gi[0 ]; - bs1 = fi[0 ] + gi[0 ]; - bc2 = fi[2 ] - gi[2 ]; - bs2 = fi[2 ] + gi[2 ]; - bc3 = fi[4 ] - gi[4 ]; - bs3 = fi[4 ] + gi[4 ]; - bc4 = fi[6 ] - gi[6 ]; - bs4 = fi[6 ] + gi[6 ]; - bf1 = (bs1 - bs2); - bf0 = (bs1 + bs2); - bg1 = (bc1 - bc2); - bg0 = (bc1 + bc2); - bf3 = (bs3 - bs4); - bf2 = (bs3 + bs4); - bg3 = SQRT2*bc4; - bg2 = SQRT2*bc3; - fi[4 ] = bf0 - bf2; - fi[0 ] = bf0 + bf2; - fi[6 ] = bf1 - bf3; - fi[2 ] = bf1 + bf3; - gi[4 ] = bg0 - bg2; - gi[0 ] = bg0 + bg2; - gi[6 ] = bg1 - bg3; - gi[2 ] = bg1 + bg3; - } - } - if (n<16) return; - - do - { - REAL s1,c1; - int ii; - k += 2; - k1 = 1 << k; - k2 = k1 << 1; - k4 = k2 << 1; - k3 = k2 + k1; - kx = k1 >> 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<fn); - TRIG_INIT(k,c1,s1); - for (ii=1;ii<kx;ii++) - { - REAL c2,s2; - TRIG_NEXT(k,c1,s1); - c2 = c1*c1 - s1*s1; - s2 = 2*(c1*s1); - fn = fz + n; - fi = fz +ii; - gi = fz +k1-ii; - do - { - REAL a,b,g0,f0,f1,g1,f2,g2,f3,g3; - b = s2*fi[k1] - c2*gi[k1]; - a = c2*fi[k1] + s2*gi[k1]; - f1 = fi[0 ] - a; - f0 = fi[0 ] + a; - g1 = gi[0 ] - b; - g0 = gi[0 ] + b; - b = s2*fi[k3] - c2*gi[k3]; - a = c2*fi[k3] + s2*gi[k3]; - f3 = fi[k2] - a; - f2 = fi[k2] + a; - g3 = gi[k2] - b; - g2 = gi[k2] + b; - b = s1*f2 - c1*g3; - a = c1*f2 + s1*g3; - fi[k2] = f0 - a; - fi[0 ] = f0 + a; - gi[k3] = g1 - b; - gi[k1] = g1 + b; - b = c1*g2 - s1*f3; - a = s1*g2 + c1*f3; - gi[k2] = g0 - a; - gi[0 ] = g0 + a; - fi[k3] = f1 - b; - fi[k1] = f1 + b; - gi += k4; - fi += k4; - } while (fi<fn); - } - TRIG_RESET(k,c1,s1); - } while (k4<n); -} - -void mayer_fft(int n, REAL *real, REAL *imag) -{ - REAL a,b,c,d; - REAL q,r,s,t; - int i,j,k; - for (i=1,j=n-1,k=n/2;i<k;i++,j--) { - a = real[i]; b = real[j]; q=a+b; r=a-b; - c = imag[i]; d = imag[j]; s=c+d; t=c-d; - real[i] = (q+t)*.5; real[j] = (q-t)*.5; - imag[i] = (s-r)*.5; imag[j] = (s+r)*.5; - } - mayer_fht(real,n); - mayer_fht(imag,n); -} - -void mayer_ifft(int n, REAL *real, REAL *imag) -{ - REAL a,b,c,d; - REAL q,r,s,t; - int i,j,k; - mayer_fht(real,n); - mayer_fht(imag,n); - for (i=1,j=n-1,k=n/2;i<k;i++,j--) { - a = real[i]; b = real[j]; q=a+b; r=a-b; - c = imag[i]; d = imag[j]; s=c+d; t=c-d; - imag[i] = (s+r)*0.5; imag[j] = (s-r)*0.5; - real[i] = (q-t)*0.5; real[j] = (q+t)*0.5; - } -} - -void mayer_realfft(int n, REAL *real) -{ - REAL a,b,c,d; - int i,j,k; - mayer_fht(real,n); - for (i=1,j=n-1,k=n/2;i<k;i++,j--) { - a = real[i]; - b = real[j]; - real[j] = (a-b)*0.5; - real[i] = (a+b)*0.5; - } -} - -void mayer_realifft(int n, REAL *real) -{ - REAL a,b,c,d; - int i,j,k; - for (i=1,j=n-1,k=n/2;i<k;i++,j--) { - a = real[i]; - b = real[j]; - real[j] = (a-b); - real[i] = (a+b); - } - mayer_fht(real,n); -} |