From: "chris perkins" <chris.perkins@canrem.com>
Subject: fft/dft progs wanted.
Reply-To: "chris perkins" <chris.perkins@canrem.com>
Organization: Canada Remote Systems
Distribution: comp
Date: 13 Oct 92 03:23:13 EST

TO:pj_pycro@csd.uwe.ac.uk

PP>I would be grateful if you could help me ... where can I find the FFT
PP>as well as DFT programs (128, 256, 512 points and perhaps more!!.)
PP>in 'C' source code and if there is any possibility for me FTP' these
PP>programs.

PP>I'm a profoundly deaf postgraduate student doing M.Sc. in Parallel Computing
PP>Systems.

I have not compiled or debugged the code which follows.  It is based
upon code I've written for some DSP chips which have been debugged.
The algorithm will be fine, but there may be some minor errors.  Hope
it helps.

Regards,
Chris

/*==============================================================================
 * CFFT.H  -  complex FFT and IFFT code
 *
 *         Copyright 1992 Chris Perkins
 *
 *============================================================================*/

void init_fft( void );
void c_fft( float *buffer, int n );
void c_ifft( float *buffer, int n );

/*==============================================================================
 * CFFT.C  -  complex FFT and IFFT code
 *
 *         Copyright 1992 Chris Perkins
 *
 *   This code implements the complex Fast Fourier Transform.
 *   The init_fft() function must be called to initialize the sine table
 *         before any call is made to c_fft() or c_ifft().
 *   Both the c_fft() and c_ifft take the same arguments.  The first argument
 *         is a pointer to an array of floats.  The second argument is an int and
 *         is the size of the COMPLEX array.  The float array contains the real
 *         components of the complex numbers in the even indices and the imaginary
 *         components in the odd indices.  The value of the second argument
 *         must therefore be half the size of the float array.
 *   The #define'd constants MAX_SIZE, SINE_SIZE, and PI_2 are related as
 *         follows:  SINE_SIZE = 3/4*MAX_SIZE; PI_2 = 1/4*MAX_SIZE.  With
 *         MAX_SIZE defined as 1024, the maximum size of the FFT is 1024.
 *
 *   NOTE:  THIS CODE HAS NOT BEEN DEBUGGED OR EVEN COMPILED.  USE AT YOUR
 *             OWN RISK.  (...and feel free to make changes.)
 *============================================================================*/

#include <math.h>
#include "cfft.h"

#define MAX_SIZE 1024
#define SINE_SIZE  768
#define PI_2       256

void bit_reverse( float *, int );
static float sines[SINE_SIZE];

/*         FFT sine table initialization */
void
init_fft( void )
{
 double two_pi;
 int ii;

 two_pi = 8.0*atan( 1.0 );
 for( ii = 0; ii < SINE_SIZE; ii++ )
   sines[ii] = sin((two_pi*ii)/MAX_SIZE);
}

/*         Forward FFT             */
void
c_fft( float *buffer, int n )
{
 int step;
 int step2;
 int index;

 bit_reverse( buffer, n );

 for( step = 2, step2 = 2*PI_2; step < 2*n; step *= 2, step2 /= 2 ) {
   for( kk = 0, index = 0; kk < 2*PI_2; kk += step2, index++ ) {
       for(ii = index, jj = index+step/2; ii < n; ii += step, jj += step) {

               buffer[jj*2] = 1/2.0*
                               (buffer[jj*2]*sines[kk+PI_2]+buffer[jj*2+1]*sines[kk]);
               buffer[jj*2+1] = 1/2.0*
                               (buffer[jj*2+1]*sines[kk+PI_2]-buffer[jj*2]*sines[kk]);
               buffer[ii*2] = 1/2.0*buffer[ii*2] + buffer[jj*2];
               buffer[jj*2] = buffer[ii*2] - 2*buffer[jj*2];
               buffer[ii*2+1] = 1/2.0*buffer[ii*2+1] + buffer[jj*2+1];
               buffer[jj*2+1] = buffer[ii*2+1] - 2*buffer[jj*2+1];
               }

(Continued in the next message)
---
  DeLuxe/386 1.25 #9629  Best wishes,         Chris Perkins
--
Canada Remote Systems  - Toronto, Ontario
World's Largest PCBOARD System - 416-629-7000/629-7044
Xref: sparky comp.dsp:2371 sci.electronics:17155