Having spent much time looking for a decent FFT and IFFT routine written in c, for use with real data, I decided to write my own. Writing a routine that worked directly on real numbers was quite a challenge, but the result was a lightning fast transform. The inverse transform took more effort, due to the fact that the lack of symmetry in the *real* FFT meant that doing the inverse required writing a whole new routine, and in fact took far longer than the effort to write the forward transform. Having put in all the effort, I have decided to put it in the public domain, without any of the restrictions that a GNU licence involves. You can find it here: http://ravellescientific.co.uk/sureal-fft/ _____________________________ Posted through www.DSPRelated.com
Public Doamin FFT software for real numbers
Started by ●March 17, 2015
Reply by ●March 17, 20152015-03-17
"MarkSci" <104213@dsprelated> writes:> Having spent much time looking for a decent FFT and IFFT routine written in > c, for use with real data, I decided to write my own. Writing a routine > that worked directly on real numbers was quite a challenge, but the result > was a lightning fast transform. The inverse transform took more effort, due > to the fact that the lack of symmetry in the *real* FFT meant that doing > the inverse required writing a whole new routine, and in fact took far > longer than the effort to write the forward transform. > > Having put in all the effort, I have decided to put it in the public > domain, without any of the restrictions that a GNU licence involves. You > can find it here: > http://ravellescientific.co.uk/sureal-fft/Didn't FFTW already provide this (and more)? -- Randy Yates Digital Signal Labs http://www.digitalsignallabs.com
Reply by ●March 17, 20152015-03-17
Randy Yates <yates@digitalsignallabs.com> writes:> "MarkSci" <104213@dsprelated> writes: > >> Having spent much time looking for a decent FFT and IFFT routine written in >> c, for use with real data, I decided to write my own. Writing a routine >> that worked directly on real numbers was quite a challenge, but the result >> was a lightning fast transform. The inverse transform took more effort, due >> to the fact that the lack of symmetry in the *real* FFT meant that doing >> the inverse required writing a whole new routine, and in fact took far >> longer than the effort to write the forward transform. >> >> Having put in all the effort, I have decided to put it in the public >> domain, without any of the restrictions that a GNU licence involves. You >> can find it here: >> http://ravellescientific.co.uk/sureal-fft/ > > Didn't FFTW already provide this (and more)?Look at the fourth bullet here: http://www.fftw.org/ -- Randy Yates Digital Signal Labs http://www.digitalsignallabs.com
Reply by ●March 17, 20152015-03-17
>Randy Yates <yates@digitalsignallabs.com> writes: > >> >> Didn't FFTW already provide this (and more)? > >Look at the fourth bullet here: http://www.fftw.org/ >-- >Randy Yates >Digital Signal Labs >http://www.digitalsignallabs.com >IIRC, FFTW is C++, my code is C, so it's easier to port to embedded systems _____________________________ Posted through www.DSPRelated.com
Reply by ●March 17, 20152015-03-17
>>Randy Yates <yates@digitalsignallabs.com> writes: >> >>> >>> Didn't FFTW already provide this (and more)? >> >>Look at the fourth bullet here: http://www.fftw.org/ >>-- >>Randy Yates >>Digital Signal Labs >>http://www.digitalsignallabs.com >> > >IIRC, FFTW is C++, my code is C, so it's easier to port to embedded >systems > >_____________________________ >Posted through www.DSPRelated.com >Ooops, I was wrong, that is C. I was getting it mixed up with another C++ version form a French guy. Yes, I vaguely at FFTW, and think it was a bit overkill for what I wanted. _____________________________ Posted through www.DSPRelated.com
Reply by ●March 18, 20152015-03-18
On 03/18/2015 06:21 AM, MarkSci wrote:>>> Randy Yates <yates@digitalsignallabs.com> writes: >>> >>>> >>>> Didn't FFTW already provide this (and more)? >>> >>> Look at the fourth bullet here: http://www.fftw.org/ >>> -- >>> Randy Yates >>> Digital Signal Labs >>> http://www.digitalsignallabs.com >>> >> >> IIRC, FFTW is C++, my code is C, so it's easier to port to embedded >> systems >> > Ooops, I was wrong, that is C. I was getting it mixed up with another C++ > version form a French guy. Yes, I vaguely at FFTW, and think it was a bit > overkill for what I wanted. >If fftw is overkill use kissfft.
Reply by ●March 19, 20152015-03-19
On 3/17/15 6:21 PM, MarkSci wrote:>>> Randy Yates<yates@digitalsignallabs.com> writes: >>> >>>> >>>> Didn't FFTW already provide this (and more)? >>> >>> Look at the fourth bullet here: http://www.fftw.org/ >>> -- >>> Randy Yates >>> Digital Signal Labs >>> http://www.digitalsignallabs.com >>> >> >> IIRC, FFTW is C++, my code is C, so it's easier to port to embedded >> systems >> >> > Ooops, I was wrong, that is C. I was getting it mixed up with another C++ > version form a French guy. Yes, I vaguely at FFTW, and think it was a bit > overkill for what I wanted. >listen, i would love to help you, but i can't seem to find this old program.... YES I DID!!!! this sucker's *OLD*. like from the 1980's. basic Cooley-Tukey Radix-2. i think i avoid multiplying by 1 and j (because in them olden daze, multiplication was more expensive). you have to unwrap lines that my email client wrapped, but otherwise this should just work. and it's simple and clean. (but it's for complex input. for real input, put the even-indexed into the real part, the odd-index into the imag part, do the FFT, and untangle the results on the other side. ___________________________________________________________________________ /* A set of utility programs to compute the Fast Fourier Transform (FFT): N-1 X[k] = SUM { x[n]exp(-j2 pi nk/N) } n=0 and inverse Fast Fourier Transform (iFFT): N-1 x[n] = 1/N SUM { X[k]exp(+j2 pi nk/N) } k=0 To speed things up, multiplication by 1 and j are avoided. The input, x[], is an array of complex numbers (pairs of doubles) of length N = 2^n. The calling program supplies n = log2(N) not the array length, N. The output is placed in BIT REVERSED order in x[]. Call bit_reverse(x, n) to swap back to normal order, if needed. However, iFFT(X, n, stbl) requires its input, X[], to be in bit reversed order making bit reversing unnecessary in some cases, such as convolution. stbl is an array of doubles of length N/4 containing the sine function from 0 to pi/2 used to compute the FFT. Call sin_table(stbl, n) ONLY ONCE before either FFT(x, n, stbl) or iFFT(X, n, stbl) to set up a sin table for FFT computation. Written in THINK C by Robert Bristow-Johnson. */ #define HALFPI 1.570796326794897 #define PI 3.141592653589793 #define TWOPI 6.283185307179586 // #include "complex.h" typedef struct { double real; double imag; } complex; #define Re(z) (z).real #define Im(z) (z).imag void FFT(complex *x, int n, double *stbl) { long size; register long length, step, stepsize, end; register complex *top, *bottom; /* top & bottom of FFT butterfly */ complex temp; size = 1L<<(n-2); end = (long)x + 4*sizeof(temp)*size; length = size; stepsize = 1; while ( length >= 1) { top = x; while ((long)top < end) { bottom = top + 2*length; Re(temp) = Re(*top) - Re(*bottom); /* butterfly: twiddle = 1 */ Im(temp) = Im(*top) - Im(*bottom); Re(*top) += Re(*bottom); Im(*top) += Im(*bottom); *bottom = temp; top++; bottom++; for (step = stepsize; step < size; step += stepsize) { Re(temp) = Re(*top) - Re(*bottom); /* butterfly: twiddle = exp(-j�) */ Im(temp) = Im(*top) - Im(*bottom); Re(*top) += Re(*bottom); Im(*top) += Im(*bottom); Re(*bottom) = Re(temp)*stbl[size - step] + Im(temp)*stbl[step]; Im(*bottom) = Im(temp)*stbl[size - step] - Re(temp)*stbl[step]; top++; bottom++; } Re(temp) = Im(*top) - Im(*bottom); /* butterfly: twiddle = -j */ Im(temp) = Re(*bottom) - Re(*top); Re(*top) += Re(*bottom); Im(*top) += Im(*bottom); *bottom = temp; top++; bottom++; for (step = stepsize; step < size; step += stepsize) { Re(temp) = Im(*top) - Im(*bottom); /* butterfly: twiddle = -j*exp(-j�) */ Im(temp) = Re(*bottom) - Re(*top); Re(*top) += Re(*bottom); Im(*top) += Im(*bottom); Re(*bottom) = Re(temp)*stbl[size - step] + Im(temp)*stbl[step]; Im(*bottom) = Im(temp)*stbl[size - step] - Re(temp)*stbl[step]; top++; bottom++; } top = bottom; } length >>= 1; stepsize <<= 1; } top = x; bottom = x +1; while ((long)top < end) { Re(temp) = Re(*top) - Re(*bottom); /* butterfly: twiddle = 1 */ Im(temp) = Im(*top) - Im(*bottom); Re(*top) += Re(*bottom); Im(*top) += Im(*bottom); *bottom = temp; top += 2; bottom += 2; } } void iFFT(complex *X, int n, double *stbl) { long size; register long length, step, stepsize, end; double scale; register complex *top, *bottom; /* top & bottom of FFT butterfly */ complex temp; size = 1L<<(n-2); end = (long)X + 4*sizeof(temp)*size; scale = 0.25/size; top = X; bottom = X +1; while ((long)top < end) { Re(temp) = (Re(*top) - Re(*bottom))*scale; /* butterfly: twiddle = 1/N */ Im(temp) = (Im(*top) - Im(*bottom))*scale; Re(*top) = (Re(*top) + Re(*bottom))*scale; Im(*top) = (Im(*top) + Im(*bottom))*scale; *bottom = temp; top += 2; bottom += 2; } length = 1; stepsize = size; while ( stepsize >= 1) { top = X; while ((long)top < end) { bottom = top + 2*length; temp = *bottom; /* butterfly: twiddle = 1 */ Re(*bottom) = Re(*top) - Re(temp); Im(*bottom) = Im(*top) - Im(temp); Re(*top) += Re(temp); Im(*top) += Im(temp); top++; bottom++; for (step = stepsize; step < size; step += stepsize) { /* butterfly: twiddle = exp(+j�) */ Re(temp) = Re(*bottom)*stbl[size - step] - Im(*bottom)*stbl[step]; Im(temp) = Im(*bottom)*stbl[size - step] + Re(*bottom)*stbl[step]; Re(*bottom) = Re(*top) - Re(temp); Im(*bottom) = Im(*top) - Im(temp); Re(*top) += Re(temp); Im(*top) += Im(temp); top++; bottom++; } Re(temp) = -Im(*bottom); /* butterfly: twiddle = +j */ Im(temp) = Re(*bottom); Re(*bottom) = Re(*top) - Re(temp); Im(*bottom) = Im(*top) - Im(temp); Re(*top) += Re(temp); Im(*top) += Im(temp); top++; bottom++; for (step = stepsize; step < size; step += stepsize) { /* butterfly: twiddle = +j*exp(+j�) */ Re(temp) = -Im(*bottom)*stbl[size - step] - Re(*bottom)*stbl[step]; Im(temp) = Re(*bottom)*stbl[size - step] - Im(*bottom)*stbl[step]; Re(*bottom) = Re(*top) - Re(temp); Im(*bottom) = Im(*top) - Im(temp); Re(*top) += Re(temp); Im(*top) += Im(temp); top++; bottom++; } top = bottom; } length <<= 1; stepsize >>= 1; } } void sin_table(double *stbl, int n) { register long size, i; double theta; size = 1L<<(n-2); theta = HALFPI/size; for (i = 0; i < size; i++) { stbl[i] = sin(theta*i); } } void bit_reverse(register complex *x, int n) { complex temp; register long k, i, r, size, count; size = (1L<<n) - 1L; for (k = 1L; k < size; k++) { i = k; r = 0; for (count = n; count > 0; count--) { r <<= 1; r += (i & 0x00000001L); i >>= 1; } if (r > k) { temp = x[r]; x[r] = x[k]; x[k] = temp; } } } ___________________________________________________________________________ -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Reply by ●March 21, 20152015-03-21
Thanks. GPL can be a problem in industrial applications. A long time ago I used FFTW in internal code that suddenly became part of a product deal. We had to replace it. The catch is this, and I can't help but grin when I read over the word "protected" :-)>> Non-free licenses may also be purchased from MIT, for users who do notwant their programs protected by the GPL _____________________________ Posted through www.DSPRelated.com
Reply by ●March 22, 20152015-03-22
On 3/21/15 3:03 PM, mnentwig wrote:> Thanks.FWIW. there is **nothing** new or spectacular (like in FFTW) in that 3 decade old FFT i posted. the only optimization i did was: 1. avoid multiplication by 1 or j and 2. do normal-in bit-reversed-out FFT and bit-reversed-in and normal-out for iFFT (that way, if you multiplied in the frequency domain, you could do that to bit-reversed data and not have to run bit_reverse()). i haven't used it in decades, but last i remember, it works. it came straight outa reading O&S. and, again, it's for complex input. for real, either zero the imaginary part or do the untangle thingie.> GPL can be a problem in industrial applications.i found gnu to be mostly useless. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Reply by ●March 22, 20152015-03-22
"mnentwig" <24789@dsprelated> writes:> Thanks. > GPL can be a problem in industrial applications. > A long time ago I used FFTW in internal code that suddenly became part of a > product deal. We had to replace it.Doesn't GPL say that if you keep it as a shared library you're free to use it for-profit? -- Randy Yates Digital Signal Labs http://www.digitalsignallabs.com
