Fixed point conversion in FFT

Hello guys,

I need some helps on fixed point conversion in FFT. I knew the fixed point theory and knew how to use it for a simple program. However, when I tried to implement it into FFT, I could not get the correct result and found that the problem was due to the bit growth in the butterfly routines. I tried to solve the problem but it did not work. Hope someone could help me. Many Thanks in advance. Below is my floating point C code:

Code:

---

#include <stdio.h>
#include <math.h>
#include <malloc.h>
#include "fftrm.h"

#define M_PI 3.1415926535897932384626433832795

float *WR;
float *WI;
doublereal *DWR;
doublereal *DWI;
void rmpo( int *rv, int *rvp)
{ int value_h;
  int n;
  n = 1;
  *rvp = -1;
  value_h = 1;
  while ( value_h > 0)
  { value_h = *rv - n;
  (*rvp)++;
    n += n;
  }
}
void zfftrmc (doublecomplex *X, int M, int P, float D)
{
 int MV2,MM1,J,I,K,L,LE,LE1,IP,IQ,IND,IND1,R;
 int I1, J1, IPOTR;
 float A,B;
 float WCOS,WSIN;
 float VR,VI;
 float ARG;
 DWR = (doublereal *)calloc(M,sizeof(doublereal));
 DWI = (doublereal *)calloc(M,sizeof(doublereal));
 LE = 1;
 IND = 0;
 for (L=1; L<=P; L++)          /* Twiddle factor routines */
 {LE1 = LE;
  LE = LE*2;
  DWR[IND] = 1.0;
  DWI[IND] = 0.0;
  ARG = (float)M_PI/(float)LE1;
  WCOS = (float)cos(ARG);
  WSIN = (float)(D*sin(ARG));
  for (R=1; R<=LE1; R++)
  {IND1 = IND+1;
  A = (float)DWR[IND];
  B = (float)DWI[IND];
  DWR[IND1] = A*WCOS - B*WSIN;
  DWI[IND1] = B*WCOS + A*WSIN;
  ++IND;
  }
  }       

  MV2 = M/2;            /* bit reversal  routines*/
  MM1 = M-1;
  J = 1;
  for (I=1; I<=MM1; I++)
  { if(I>=J)
    goto P1;
        J1 = J-1;
        I1 = I-1;
        VR = (float)RE(X[J1]);
        VI = (float)IM(X[J1]);
        RE(X[J1]) = RE(X[I1]);
        IM(X[J1]) = IM(X[I1]);
        RE(X[I1]) = VR;
        IM(X[I1]) = VI;
        P1:K = MV2;
        P2: if(K>=J) goto P3;
        J = J-K;
        K = K/2;
        goto P2;
        P3: J = J+K;
  }
  IND = 0;
  LE = 1;
  for (L=1; L<=P; L++)    /* Number of stage  */
  { LE1 = LE;
    LE = LE*2;
    for (R=0; R<LE1; R++)
        { WCOS = (float)DWR[IND];
          WSIN = (float)DWI[IND];
          IND = IND+1;
          for (IQ=R; IQ<M; IQ+=LE)          /* Butterfly routines */
          {    IP = IQ+LE1;
                A = (float)RE(X[IP]);
                B = (float)IM(X[IP]);
                VR = A*WCOS - B*WSIN;
                VI = B*WCOS + A*WSIN;
                RE(X[IP]) = RE(X[IQ]) - VR;
                IM(X[IP]) = IM(X[IQ]) - VI;
                RE(X[IQ]) = RE(X[IQ]) + VR;
                IM(X[IQ]) = IM(X[IQ]) + VI;
          }
        }
  }
  free(DWR);
  free(DWI);
}

/* fftrm.h */

#ifndef FFTRM_H
#define FFTRM_H
#define RE(z) ((z).r)
#define IM(z) ((z).i)
typedef float doublereal;
typedef struct { doublereal r, i;} doublecomplex;
int zffts (doublecomplex *X, int M);
int ziffts (doublecomplex *X, int M);
void zfftrmc (doublecomplex *X, int M, int P, float D);
void rmpo (int *rv, int *rvp);
#endif

---

Thanks a lot.