/* Appended is the original of Ron Mayers FHT code (fft_may.c trigtbl.h) and some of his remarks. */ /* ** FFT and FHT routines ** Copyright 1988, 1993; Ron Mayer ** ** fht(fz,n); ** Does a hartley transform of "n" points in the array "fz". ** fft(n,real,imag) ** Does a fourier transform of "n" points of the "real" and ** "imag" arrays. ** ifft(n,real,imag) ** Does an inverse fourier transform of "n" points of the "real" ** and "imag" arrays. ** realfft(n,real) ** Does a real-valued fourier transform of "n" points of the ** "real" and "imag" arrays. 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. ** 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 ** */ #define GOOD_TRIG #include "trigtbl.h" char fht_version[] = "Brcwl-Hrtly-Ron-dbld"; #define SQRT2_2 0.70710678118654752440084436210484 #define SQRT2 2*0.70710678118654752440084436210484 fht(fz,n) int n; REAL *fz; { REAL a,b; REAL c1,s1,s2,c2,s3,c3,s4,c4; REAL f0,g0,f1,g1,f2,g2,f3,g3; int i,k,k1,k2,k3,k4,kx; REAL *fi,*fn,*gi; TRIG_VARS; for (k1=1,k2=0;k1>1; (!((k2^=k)&k)); k>>=1); if (k1>k2) { a=fz[k1];fz[k1]=fz[k2];fz[k2]=a; } } 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 (fi1) \ { \ for (j=k-i+2 ; (1< real_valued_FFT is faster than the transformation from 1/2pointFFT-> real_valued_FFT so my real-valued fft is even better when compared to most published real valued ffts. 2) Avoid multiplications by 1 and zero (and sometimes sqrt2). Many published routines such as Numerical Recipes seem to spend a lot of time multiplying by cos(0), cos(pi), etc. and almost all seem to do 1/sqrt_2*x+1/sqrt_2*y instead of 1/sqrt_2*(x+y). 3) Faster trig generation. Most algorithms use 1 "sin" library call for each level of the transform; and 2 real multiplications and 2 real additions for each needed trig value within it's loop. I use a stable algorithm to generate each trig value using 1 real multiplication and 1 real addition for each value using a small (log(n)) table of trig values. The tradeoff is that I require much more integer arithmetic for this calculation, including a (n*log(n)) loop; but for multiples of pi/16 or so, my routine still seems faster. By taking advantage of the fact that values required for FFTs are for evenly spaced angles, I avoid all calls to slow trig library functions which are unnecessarily complex because they need to work for arbitrary values. 4) Generate less trig values I use the identities sin(x)=sin(pi-x)=-sin(pi+x)=-sin(-x),etc. to avoid excessive trig calculations, and sin(2x) = 2*cos(x)*sin(x) to allow simpler trig calculations when accuracy permits. A more stable than average trig generator mentioned in (3) above allows me to use the unstable sin(2x) = 2*cos(x)*sin(x) hack for every other "level" in the FFT without the usual loss of accuracy. 5) Mixed 2,4-radix inner loop. By doing two levels in the inner loop, I gain all the advantages of a radix-4 algorithm; primarily reducing integer arithmetic and memory access. This has a great affect on large arrays when paging occurs. 6) Unrolling loops and variable storage to optimize register allocation. I try not to require storing too many values in variables at any one time to ease a compilers register allocation. 7) Remove low levels of the transform out of the loop. It's significantly faster to do 8 point transforms explicitly; rather than using the general loop. One catch to this routine is that at least two of the algorithms used by it are patented(!) (just about any FHT is patented by R. Bracewell; and the stable trig generator is patented by O. Buneman; both at Stanford Univ.) Who owns the copyright rights to it is also probably being debated; since it is a derivative work of a GNU-licensed routine, so subject to their restrictions; it was worked on for a Stanford project, so they have a claim on it; and I optimized it further working for this company, so they probably claim parts of it. Considering Gauss apparently used the equivalent of real valued FFTs in 1805; and Euler did fourier transforms it in the mid 1700s; I'm amazed that people still want to claim this math. Ron Mayer mayer@acuson.com */