fourier.cc

#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#define SWAP(a,b) tempr=(a); (a)=(b); (b)=tempr
const long N = 1 << 13;

void four(float*, unsigned long, int);

int main(int argc, char* argv[]) {
  char* buf = new char[256];
  long nn=N;

  if (argc == 2) {
    int n = atoi(argv[1]);
    if (2<n && n<21) nn = 1 << n;
    else {
      perror("number must be 3-20\n");
      exit(1);
    }
  }

  float* data = new float[2*(nn+1)];
  float x, y;
  for (long i=0; i<nn; i++) {
    data[2*i]=0;
    data[2*i+1]=0;
  }
  
  long m=0;
  while (NULL!=fgets(buf, 255, stdin) || m>=nn) {
    if (buf[0] == '#') continue;
    sscanf(buf, "%f %f", &x, &y);
    data[2*m] = y;
    m++;
  }

  four(data, nn, 1);
  for (unsigned long i=(int)nn/m; i<=nn/2; i++) {
    printf("%f\t%f\n", (float)nn/i,
	   data[2*i]*data[2*i]+data[2*i+1]*data[2*i+1]);
  }
  exit(0);
}

//--------------- fourier converter
void four(float data[], unsigned long nn, int isign) {
  unsigned long n, mmax, m, j, istep, i;
  double wtemp, wr, wpr, wpi, wi, theta;
  float tempr, tempi;

  n = nn << 1;
  j = 1;
  for (i=1; i<n; i+=2) {
    if (j > i) {
      SWAP(data[j], data[i]);
      SWAP(data[j+1], data[i+1]);
    }
    m = n >> 1;
    while (m >= 2 && j > m) {
      j -= m;
      m >>= 1;
    }
    j += m;
  }
  
  mmax = 2;
  while (n > mmax) {
    istep = mmax << 1;
    theta = isign * (6.28318530717959/mmax);
    wtemp = sin(0.5*theta);
    wpr = -2.0 * wtemp * wtemp;
    wpi = sin(theta);
    wr = 1.0;
    wi = 0.0;
    for (m=1; m<mmax; m+=2) {
      for (i=m; i<=n; i+=istep) {
	j = i+mmax;
	tempr = wr*data[j]   - wi*data[j+1];
	tempi = wr*data[j+1] + wi*data[j];
	data[j]   = data[i]   - tempr;
	data[j+1] = data[i+1] - tempi;
	data[i]   += tempr;
	data[i+1] += tempi;
      }
      wr += (wtemp=wr)*wpr - wi*wpi;
      wi += wi*wpr + wtemp*wpi;
    }
    mmax = istep;
  }
}