Reply by glen herrmannsfeldt●March 23, 20152015-03-23
Evgeny Filatov <e.v.filatov@ieee.org> wrote:
> On 23.03.2015 3:35, Randy Yates wrote:
>> "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?
> Nope. The license which allows that is LGPL (Lesser GPL).
It depends on the meaning of "mere aggregation".
I won't try to define it.
-- glen
Reply by Evgeny Filatov●March 23, 20152015-03-23
On 23.03.2015 3:35, Randy Yates wrote:
> "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?
>
Nope. The license which allows that is LGPL (Lesser GPL).
Evgeny.
Reply by Randy Yates●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
Reply by robert bristow-johnson●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 mnentwig●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 not
want their programs protected by the GPL
_____________________________
Posted through www.DSPRelated.com
Reply by robert bristow-johnson●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 Steve Underwood●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 MarkSci●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 MarkSci●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 Randy Yates●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)?