I'm having problems with FFT

[secondrate]

Hi guys,

Would someone be nice enough to convert this C code to C++ code for me? I don't know how to deal with Python so I had difficulty doing the conversion.

#include <stdlib.h>
#include <math.h>

/* macros */
#define TWO_PI (6.2831853071795864769252867665590057683943L)

/* function prototypes */
void fft(int N, double (*x)[2], double (*X)[2]);
void fft_rec(int N, int offset, int delta,
             double (*x)[2], double (*X)[2], double (*XX)[2]);
void ifft(int N, double (*x)[2], double (*X)[2]);

/* FFT */
void fft(int N, double (*x)[2], double (*X)[2])
{
  /* Declare a pointer to scratch space. */
  double (*XX)[2] = malloc(2 * N * sizeof(double));

  /* Calculate FFT by a recursion. */
  fft_rec(N, 0, 1, x, X, XX);

  /* Free memory. */
  free(XX);
}

/* FFT recursion */
void fft_rec(int N, int offset, int delta,
             double (*x)[2], double (*X)[2], double (*XX)[2])
{
  int N2 = N/2;            /* half the number of points in FFT */
  int k;                   /* generic index */
  double cs, sn;           /* cosine and sine */
  int k00, k01, k10, k11;  /* indices for butterflies */
  double tmp0, tmp1;       /* temporary storage */

  if(N != 2)  /* Perform recursive step. */
    {
      /* Calculate two (N/2)-point DFT's. */
      fft_rec(N2, offset, 2*delta, x, XX, X);
      fft_rec(N2, offset+delta, 2*delta, x, XX, X);

      /* Combine the two (N/2)-point DFT's into one N-point DFT. */
      for(k=0; k<N2; k++)
        {
          k00 = offset + k*delta;    k01 = k00 + N2*delta;
          k10 = offset + 2*k*delta;  k11 = k10 + delta;
          cs = cos(TWO_PI*k/(double)N); sn = sin(TWO_PI*k/(double)N);
          tmp0 = cs * XX[k11][0] + sn * XX[k11][1];
          tmp1 = cs * XX[k11][1] - sn * XX[k11][0];
          X[k01][0] = XX[k10][0] - tmp0;
          X[k01][1] = XX[k10][1] - tmp1;
          X[k00][0] = XX[k10][0] + tmp0;
          X[k00][1] = XX[k10][1] + tmp1;
        }
    }
  else  /* Perform 2-point DFT. */
    {
      k00 = offset; k01 = k00 + delta;
      X[k01][0] = x[k00][0] - x[k01][0];
      X[k01][1] = x[k00][1] - x[k01][1];
      X[k00][0] = x[k00][0] + x[k01][0];
      X[k00][1] = x[k00][1] + x[k01][1];
    }
}

/* IFFT */
void ifft(int N, double (*x)[2], double (*X)[2])
{
  int N2 = N/2;       /* half the number of points in IFFT */
  int i;              /* generic index */
  double tmp0, tmp1;  /* temporary storage */

  /* Calculate IFFT via reciprocity property of DFT. */
  fft(N, X, x);
  x[0][0] = x[0][0]/N;    x[0][1] = x[0][1]/N;
  x[N2][0] = x[N2][0]/N;  x[N2][1] = x[N2][1]/N;
  for(i=1; i<N2; i++)
    {
      tmp0 = x[i][0]/N;       tmp1 = x[i][1]/N;
      x[i][0] = x[N-i][0]/N;  x[i][1] = x[N-i][1]/N;
      x[N-i][0] = tmp0;       x[N-i][1] = tmp1;
    }
}

and

#include <stdio.h>
#include <stdlib.h>
#include "fft.c"

int main()
{
  int i;                    /* generic index */
  char file[FILENAME_MAX];  /* name of data file */
  int N;                    /* number of points in FFT */
  double (*x)[2];           /* pointer to time-domain samples */
  double (*X)[2];           /* pointer to frequency-domain samples */
  double dummy;             /* scratch variable */
  FILE *fp;                 /* file pointer */

  /* Get name of input file of time-domain samples x. */
  printf("Input file for time-domain samples x(n)? ");
  scanf("%s", file);

  /* Read through entire file to get number N of points in FFT. */
  if(!(fp = fopen(file, "r")))
    {
      printf("   File \'%s\' could not be opened!", file);
      exit(EXIT_FAILURE);
    }
  N=0;
  while(fscanf(fp, "%lg%lg", &dummy, &dummy) == 2) N++;
  printf("N = %d", N);

  /* Check that N = 2^n for some integer n >= 1. */
  if(N >= 2)
    {
      i = N;
      while(i==2*(i/2)) i = i/2;  /* While i is even, factor out a 2. */
    }  /* For N >=2, we now have N = 2^n iff i = 1. */
  if(N < 2 || i != 1)
    {
      printf(", which does not equal 2^n for an integer n >= 1.");
      exit(EXIT_FAILURE);
    }

  /* Allocate time- and frequency-domain memory. */
  x = malloc(2 * N * sizeof(double));
  X = malloc(2 * N * sizeof(double));

  /* Get time-domain samples. */
  rewind(fp);
  for(i=0; i<N; i++) fscanf(fp, "%lg%lg", &x[i][0], &x[i][1]);
  fclose(fp);

  /* Calculate FFT. */
  fft(N, x, X);

  /* Print time-domain samples and resulting frequency-domain samples. */
  printf("\nx(n):");
  for(i=0; i<N; i++) printf("\n   n=%d: %12f %12f", i, x[i][0], x[i][1]);
  printf("\nX(k):");
  for(i=0; i<N; i++) printf("\n   k=%d: %12f %12f", i, X[i][0], X[i][1]);

  /* Clear time-domain samples and calculate IFFT. */
  for(i=0; i<N; i++) x[i][0] = x[i][1] = 0;
  ifft(N, x, X);

  /* Print recovered time-domain samples. */
  printf("\nx(n):");
  for(i=0; i<N; i++) printf("\n   n=%d: %12f %12f", i, x[i][0], x[i][1]);

  /* Write frequency-domain samples X to a file, if desired. */
  printf("\nOutput file for frequency-domain samples X(k)?"
         "\n   (if none, abort program): ");
  scanf("%s", file);
  if(!(fp = fopen(file, "w")))
    {
      printf("   File \'%s\' could not be opened!", file);
      exit(EXIT_FAILURE);
    }
  for(i=0; i<N; i++) fprintf(fp, "%23.15e  %23.15e\n", X[i][0], X[i][1]);
  fclose(fp);
  printf("Samples X(k) were written to file %s.", file);

  /* Free memory. */
  free(x);
  free(X);

  return 0;
}

I'd really really appreciate your help. Thanks so much!