[C11][FFTW] FFT examples

main-fftw-bench.c

// [compile] cc -O3 -std=c11 -pedantic -Wall -Wextra main-fftw-bench.c -lfftw3 -o main-fftw-bench

#include <complex.h>
#include <stdio.h>
#include <stdlib.h>
#include <fftw3.h>

typedef fftw_complex CN;
void fft(size_t N, CN c[N], CN ret[N]) {
  fftw_plan fft = fftw_plan_dft_1d(N, c, ret, FFTW_FORWARD, FFTW_ESTIMATE);
  fftw_execute(fft);
  fftw_destroy_plan(fft);  
}
void ifft(size_t N, CN c[N], CN ret[N]) {
  fftw_plan ifft = fftw_plan_dft_1d(N, c, ret, FFTW_BACKWARD, FFTW_ESTIMATE);
  fftw_execute(ifft);
  fftw_destroy_plan(ifft);
  for (size_t i = 0; i < N; i++) ret[i] /= N;
}

// example main
static void print_cn(CN c) {
  printf("%f%s%fi", creal(c), cimag(c) < 0.0 ? "" : "+", cimag(c));
}
static void print_veccn(size_t N, const CN cs[]) {
  puts("[");
  for (size_t i = 0; i < N; i++) {
    printf("  ");
    print_cn(cs[i]);
    puts(",");
  }
  puts("]");
}
static double drandom() {
  return rand() / (double) RAND_MAX * 2 - 1;
}
int main() {
  srand(0);
  //const size_t N = 1024 * 1024 * 4; // for fft only
  const size_t N = 1024 * 1024 * 8; // for fft only
  //const size_t N = 1024 * 1024 + 1; // for fft only
  //const size_t N = 1024 * 1024 * 16; // for fft only
  //const size_t N = 1024 + 1;
  //const size_t N = 1024;
  CN *raw = fftw_alloc_complex(N);
  for (size_t i = 0; i < N; i++) raw[i] = drandom();

  CN *ffreq = fftw_alloc_complex(N), *fseq = fftw_alloc_complex(N);
  fft(N, raw, ffreq);
  ifft(N, ffreq, fseq);
  
  if (0) {printf("fft: "); print_veccn(N, ffreq);}
  if (0) {printf("ifft: "); print_veccn(N, fseq);}

  double fftd = 0.0;
  for (size_t i = 0; i < N; i++) fftd += cabs(raw[i] - fseq[i]);
  printf("raw distance: %.15f\n", fftd / N);

  fftw_free(raw), fftw_free(ffreq), fftw_free(fseq);
  return 0;
}

main-fftw-example.c

// [compile] cc -std=c11 -pedantic -Wall -Wextra main-fftw-example.c -lfftw3 -o main-fftw-example

#include <complex.h>
#include <stdio.h>
#include <fftw3.h>

typedef fftw_complex CN;
void fft(size_t N, CN c[N], CN ret[N]) {
  fftw_plan fft = fftw_plan_dft_1d(N, c, ret, FFTW_FORWARD, FFTW_ESTIMATE);
  fftw_execute(fft);
  fftw_destroy_plan(fft);
}
void ifft(size_t N, CN c[N], CN ret[N]) {
  fftw_plan ifft = fftw_plan_dft_1d(N, c, ret, FFTW_BACKWARD, FFTW_ESTIMATE);
  fftw_execute(ifft);
  fftw_destroy_plan(ifft);
  for (size_t i = 0; i < N; i++) ret[i] /= N;
}

// example main
static void print_cn(CN c) {
  printf("%f%s%fi", creal(c), cimag(c) < 0.0 ? "" : "+", cimag(c));
}
static void print_veccn(size_t N, const CN cs[]) {
  puts("[");
  for (size_t i = 0; i < N; i++) {
    printf("  ");
    print_cn(cs[i]);
    puts(",");
  }
  puts("]");
}
int main() {
  //const size_t N = 8; CN raw[] = {0, 1, 2, 3, 4, 5, 6, 7};
  const size_t N = 15; CN raw[] = {1,3,4,2, 5,6,2,4, 0,1,3,4, 5,62,2};
  //const size_t N = 16; CN raw[] = {1,3,4,2, 5,6,2,4, 0,1,3,4, 5,62,2, 3};

  CN ffreq[N], fseq[N];
  fft(N, raw, ffreq);
  ifft(N, ffreq, fseq);

  printf("fft: "); print_veccn(N, ffreq);
  printf("ifft: "); print_veccn(N, fseq);
  return 0;
}