I need a FFT algorithm that doesn't use SIN, COS, EXP function.

Hi,

I've a FFT problem. This is the situation:
 1.I'm working with a DSP56824 (Motorola).
 2.I need to do a FFT (maybe 512 or 1024 points).
 3.I'm working with a C compiler (Codewarrior, but I DON'T HAVE SDK).
 4.I was trying to do a FFT with standart FFT C library. Here I have
problems, because, the algorithm calls to sin or cos function, but the
compiler doesn't support it, I think thta I need to work with data
tables. The question is:

Someone has a FFT library or algorithm, that doesn't uses sin or cos
(it would have a data table)?

Or if someone knows where I can download SKD for motorola dsp56824
FREE?

Or if you have another solution or another probed algorithm, PLEASE
send me it.

Thank you.

Diego.

P.S. Excuse my english. Thank you again.

For FFT usually people do not use SIN and COS function in practice.
They use simple table look-up. First you create a table with either
half cycle or 1/4 of cycle, then you look up the table with an index to
an array.

You can find the algorithm from many DSP tutorials talking about this
implementation.

dfalvarado@utpl.edu.ec wrote:
> Hi,
> 
> I've a FFT problem. This is the situation:
>  1.I'm working with a DSP56824 (Motorola).
>  2.I need to do a FFT (maybe 512 or 1024 points).
>  3.I'm working with a C compiler (Codewarrior, but I DON'T HAVE SDK).
>  4.I was trying to do a FFT with standart FFT C library. Here I have
> problems, because, the algorithm calls to sin or cos function, but the
> compiler doesn't support it, I think thta I need to work with data
> tables. The question is:
> 
> Someone has a FFT library or algorithm, that doesn't uses sin or cos
> (it would have a data table)?
> 
> Or if someone knows where I can download SKD for motorola dsp56824
> FREE?
> 
> Or if you have another solution or another probed algorithm, PLEASE
> send me it.
> 
> Thank you.
> 
> Diego.
> 
> P.S. Excuse my english. Thank you again.
> 

Are you using the compiler right?  I'm pretty sure sin/cos/exp are 
defined in the standard C math library, but you may have to explicity 
link to it via a build parameter (e.g. -lm)

If that is still no help, you can avoid the math functions by making a 
hard coded fft structure using kissfft.  Precompute what the structure 
values would be, and then compile them into a header file.

http://www.sourceforge.net/projects/kiss_fft


-- Mark

Mark Borgerding wrote:

> http://www.sourceforge.net/projects/kiss_fft

That should be
http://www.sourceforge.net/projects/kissfft

<dfalvarado@utpl.edu.ec> wrote in message
news:1115853350.319429.148020@o13g2000cwo.googlegroups.com...
> Hi,
>
> I've a FFT problem. This is the situation:
>  1.I'm working with a DSP56824 (Motorola).
>  2.I need to do a FFT (maybe 512 or 1024 points).
>  3.I'm working with a C compiler (Codewarrior, but I DON'T HAVE SDK).
>  4.I was trying to do a FFT with standart FFT C library. Here I have
> problems, because, the algorithm calls to sin or cos function, but the
> compiler doesn't support it, I think thta I need to work with data
> tables. The question is:
>
> Someone has a FFT library or algorithm, that doesn't uses sin or cos
> (it would have a data table)?

I'm fairly sure that Motorola has a canned FFT routine written in assembly
for you in one of it's libraries. You can call it directly from C or write a
wrapper to be more C friendly and call it.


> Or if someone knows where I can download SKD for motorola dsp56824
> FREE?
>
> Or if you have another solution or another probed algorithm, PLEASE
> send me it.
>
> Thank you.
>
> Diego.
>
> P.S. Excuse my english. Thank you again.
>

dfalvarado@utpl.edu.ec wrote:

> I've a FFT problem. This is the situation:
>  1.I'm working with a DSP56824 (Motorola).
>  2.I need to do a FFT (maybe 512 or 1024 points).
>  3.I'm working with a C compiler (Codewarrior, but I DON'T HAVE SDK).
>  4.I was trying to do a FFT with standart FFT C library. Here I have
> problems, because, the algorithm calls to sin or cos function, but the
> compiler doesn't support it, I think thta I need to work with data
> tables. The question is:

> Someone has a FFT library or algorithm, that doesn't uses sin or cos
> (it would have a data table)?

I believe that most FFT routines derive them through the half angle
formula, which uses sqrt() and some simple arithmetic.  Otherwise,
you can supply a table.

-- glen

On Fri, 13 May 2005 13:15:50 -0700, glen herrmannsfeldt
<gah@ugcs.caltech.edu> wrote:
> dfalvarado@utpl.edu.ec wrote:
> 
>> I've a FFT problem. This is the situation:
>>  1.I'm working with a DSP56824 (Motorola).
>>  2.I need to do a FFT (maybe 512 or 1024 points).
>>  3.I'm working with a C compiler (Codewarrior, but I DON'T HAVE SDK).
>>  4.I was trying to do a FFT with standart FFT C library. Here I have
>> problems, because, the algorithm calls to sin or cos function, but the
>> compiler doesn't support it, I think thta I need to work with data
>> tables. The question is:
> 
>> Someone has a FFT library or algorithm, that doesn't uses sin or cos
>> (it would have a data table)?
> 
> I believe that most FFT routines derive them through the half angle
> formula, which uses sqrt() and some simple arithmetic.  Otherwise,
> you can supply a table.
> 
> -- glen
> 

Motorola probably publishes such a table.  TI publishes one for every
DSP they make (two for the 'c6x generation with selectable endianism)

> Or if you have another solution or another probed algorithm, PLEASE
> send me it.
> Diego.

Here is another one that is even simpler than kissfft.
Not necessarily the fastest fft in town, but it does
the job and doesn't need any sine or cosine tables
other than the constants that are inside the program.

------------------------------------------------------------------------------


/* dft.c        Paul Mennen paul@mennen.org

fftc(*blk,sz)         Complex-Complex DFT   sz = # of complex inputs/outputs
iftc(*blk,sz)         Complex-Complex IFT   sz = # of complex inputs/outputs
fftr(*blk,sz)         Real-Complex DFT      sz = # of real inputs
iftr(*blk,sz)         Complex-Real IFT      sz = # of real outputs

*/

#include "stdio.h"
#include "math.h"

struct blk {float r; float i;};

#define RMS 1.414213562         /* adjust FFT results to read rms */

float tcos[] =
{-1.0,0.0,0.7071068,                  /* cos(pi), cos(pi/2), cos(pi/4) */
 0.9238795,0.9807853,0.9951847,       /* cos(pi/8), ... cos(pi/32) */
 0.9987955,0.9996988,0.9999247,       /* cos(pi/64), ... cos(pi/256) */
 0.9999812,0.9999952938,0.9999988234, /* cos(pi/512), ... cos(pi/2048) */
 0.9999997055,0.9999999264 };         /* cos(pi/4096), cos(pi/8192) */

float tsin[] =
{0.0,1.0,0.7071068,                   /* sin(pi), sin(pi/2), sin(pi/4) */
 0.3826834,0.1950903,0.0980171,       /* sin(pi/8), ... sin(pi/32) */
 4.9067674e-2,2.4541228e-2,1.2271538e-2,     /* sin(pi/64), ... sin(pi/256) 
*/
 6.1358846e-3,3.067956691e-3,1.533980094e-3, /* sin(pi/512), ... 
sin(pi/2048) */
 7.669902711e-4, 3.834950931e-4 };           /* sin(pi/4096), sin(pi/8192) 
*/

void fft(struct blk *x, int sz, int dir)
/* In-place complex to complex transform. From Oppenheim & Schafer, pg 332 
*/
/* dir=-1/1 DFT/IFT, sz = Number of complex elements, must be a power of 2 
*/
{ float Wr,Wi,Ur,Ui,Tr,Ti,Vr,Vi;
  int k,mm,le,le1,j,i,ii,nv2,m;

  m=0; le=j=sz;  while (j >>= 1) m++; /* m = number of folds = log2(sz) */
  for (mm=1; mm<=m; mm++)  /* do m folds */
  { le1 = le >> 1;         /* le1 = le/2 */
    Ur = 1.0;   Ui = 0.0;  /* initial phasor */
    Wr = tcos[m-mm];       /* Wr=cos(PI/le1); (phasor increment) */
    Wi = dir*tsin[m-mm];   /* Wi=dir*sin(PI/le1); (phasor increment) */
    for (j=0; j<le1; j++)  /* do sz/2 butterflys in le/2 groups */
    { i = j;               /* i points to one butterfly leg */
      while (i < sz)       /* do sz/le butterflys */
      { ii = i + le1;      /* ii points to the other butterfly leg */
        Tr = x[i].r + x[ii].r;    Ti = x[i].i + x[ii].i;
        Vr = x[i].r - x[ii].r;    Vi = x[i].i - x[ii].i;
        if (mm == m) { x[ii].r = Vr;  x[ii].i = Vi; }   /* last fold */
        else         { x[ii].r = Vr*Ur - Vi*Ui;   x[ii].i = Vr*Ui + Vi*Ur; }
        x[i].r = Tr;    x[i].i = Ti;
        i += le;
      }
      Tr = Ur;  Ur = Ur*Wr - Ui*Wi;  Ui = Ui*Wr + Tr*Wi; /* update phasor */
      }    /* end for (j=0; j<le1; j++) */
    le = le1;
  }     /* end for (mm=1; mm<=m; mm++) */

  /* bit reversal */
  m = sz >> 1;  j = 0;
  for (i=0; i<sz-1; i++)
  { if (i < j) { Tr = x[j].r;      Ti = x[j].i;     /* swap i & j */
                 x[j].r = x[i].r;  x[j].i = x[i].i;
                 x[i].r = Tr;      x[i].i = Ti;
               }
    k = m;  while (k <= j) { j -= k;  k >>= 1; }
    j += k;
  }
}

void hfold(struct blk *x, int n, int dir)     /* dir=-1/1 for DFT/IFT */
{ float Wr,Wi,Ur,Ui,ang,Zr,Zi,Xr,Xi,Yr,Yi;
  int k,n2,nk;
  k = 0;  n2 = n;  while (n2 >>= 1) k++;  /* k = number of folds */
  /*        Initial phasor (Ur,Ui)           Phasor increment (Wr,Wi) */
  /* DFT:  -cos(PI/n) - i*sin(PI/n)          cos(PI/n) + i*sin(PI/n)  */
  /* IFT:   cos(PI/n) - i*sin(PI/n)          cos(PI/n) - i*sin(PI/n)  */
  Ur = dir * (Wr = tcos[k]);   Wi = dir * (Ui = -tsin[k]);
  n2 = n >> 1;
  for (k=1; k<=n2; k++)
  { Zr = x[k].r + x[nk=n-k].r;    Zi = x[k].i - x[nk].i;
    Xr = x[k].r - x[nk].r;        Xi = x[k].i + x[nk].i;
    Yr = Ui*Xr - Ur*Xi;           Yi = Ui*Xi + Ur*Xr;
    x[k].r  = Zr + Yr;            x[k].i  = Zi + Yi;
    x[nk].r = Zr - Yr;            x[nk].i = Yi - Zi;
    Zr = Ur;  Ur = Ur*Wr - Ui*Wi;  Ui = Ui*Wr + Zr*Wi; /* update phasor */
  }
}

/* fftc(*blk,sz)   Complex-Complex DFT   sz = # of complex inputs/outputs */
void fftc(struct blk *x, int sz)
{ float sc;
  fft(x,sz,-1);
  sc = 1.0/sz;
  while (--sz >= 0) { x[sz].r *= sc;  x[sz].i *= sc; }
}

/* iftc(*blk,sz)   Complex-Complex IFT   sz = # of complex inputs/outputs */
void iftc(struct blk *x, int sz)
{ fft(x,sz,1); }

/* fftr(*blk,sz)   Real-Complex DFT      sz = # of real inputs */
void fftr(struct blk *x, int sz)
{ float sc;
  sc = (0.5*RMS)/sz; /* scale fft results (as rms single sided spectrum) */
  fft(x,sz>>=1,-1);
  hfold(x,sz,-1);       /* hermitian fold */
  x[0].r = RMS * (x[0].r + x[0].i);
  x[0].i = 0.0;         /* don't compute fmax term */
  while (--sz>=0) { x[sz].r *= sc;  x[sz].i *= sc; }
}

/* iftr(*blk,sz)   Complex-Real IFT      sz = # of real outputs */
void iftr(struct blk *x, int sz)
{ float sc;
  hfold(x,sz>>=1,1);  /* hermitian fold */
  x[0].i = x[0].r;
  fft(x,sz,1);
  sc = 1.0/RMS;   /* scale result (assume rms single sided spectrum) */
  while (--sz>=0) { x[sz].r *= sc;  x[sz].i *= sc; } /* scale result */
}

Thank you for your answers.

I have another problem.

I can't use the math.h library, for this reason I can't use the sin cos
function. 

Diego

In article <1115853350.319429.148020@o13g2000cwo.googlegroups.com>,
 dfalvarado@utpl.edu.ec wrote:

>  4.I was trying to do a FFT with standart FFT C library. Here I have
> problems, because, the algorithm calls to sin or cos function, but the
> compiler doesn't support it, I think thta I need to work with data
> tables. The question is:

Use CORDIC to generate sin and cos in fixed-point:
   http://www.worldserver.com/turk/opensource/index.html#CORDIC