cfftwtest.c 4.93 KB
Newer Older
1 2 3 4 5 6
/*! \file cfftwtest.c
 * \brief a small test program for fftw (implementation)
 * 
 * ----------------------------------------------------------------------------
 * 
 * ----
7
 * libfourier is free software; you can redistribute it and/or modify
8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version. 
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 * ----
 *
 * \author Thomas Forbriger
 * \date 12/09/2007
 * 
 * a small test program for fftw (implementation)
 * 
 * Copyright (c) 2007 by Thomas Forbriger (BFO Schiltach) 
 * 
 * REVISIONS and CHANGES 
 *  - 12/09/2007   V1.0   Thomas Forbriger
 *
 *  Note:
 *  The library code is migrated to FFTW3. This program however was just meant
 *  to be a testing sandbox to develop the first library code. It is not used
 *  in production nor is it used to test the library. This code will not be
 *  migrated and becomes obsolete.
 * 
 * ============================================================================
 */
#define TF_CFFTWTEST_C_VERSION \
  "TF_CFFTWTEST_C   V1.0   "

#include<drfftw.h>
#include<stdio.h>
#include<math.h>

void fill(fftw_real* a, int n, int m)
{
  int k;
  double f1, f2;
  f1=3.141592653589/n;
  f2=f1*m*2;
  /* fill input array */
  for (k=0; k<n; ++k)
  {
    a[k]=sin(k*f1)*sin(k*f1)*sin(k*f2);
  }
}

/*
 * the code provided in the tutorial for rfftw_one produces segmentation
 * faults if used together with heap memory
 *
 * solution: always use drfftw.h together with double precision transforms!
 */

void process(int n, int m)
{
  // create arrays 
  fftw_real* in=malloc(n*sizeof(fftw_real));
  fftw_real* out=malloc(n*sizeof(fftw_real));
  fftw_real* power_spectrum=malloc((n/2+1)*sizeof(fftw_real));
  if (in == NULL) { printf("could not alloc in\n"); abort(); }
  if (out == NULL) { printf("could not alloc out\n"); abort(); }
  if (power_spectrum == NULL) 
  { printf("could not alloc power_spectrum\n"); abort(); }
  rfftw_plan p;
  int k;
  fill(in, n, m);
  p = rfftw_create_plan(n, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
  rfftw_one(p, in, out);
  power_spectrum[0] = out[0]*out[0];  // DC component 
  for (k = 1; k < (n+1)/2; ++k)  // (k < N/2 rounded up) 
       power_spectrum[k] = out[k]*out[k] + out[n-k]*out[n-k];
  if (n % 2 == 0) // N is even 
       power_spectrum[n/2] = out[n/2]*out[n/2];  // Nyquist freq. 
  rfftw_destroy_plan(p);
  printf("\nPower spectrum of %d point signal with %d periods\n", n, m);
  printf("\ni.e. signal frequency %f\n", (double)m/n);
  for (k=0; k<((n+1)/2); ++k)
  {
    printf("  %3.3d, %5f: %10e\n",k, (double)k/n, power_spectrum[k]);
  }
  free(in);
  free(out);
  free(power_spectrum);
}


#define N 2048
int main()
{
  printf("%s\n", TF_CFFTWTEST_C_VERSION);

/*
 * the code provided in the tutorial for rfftw_one produces segmentation
 * faults if used together with stack memory for N=2000 or larger
 *
 * segfaults can be prevented if in and out array dimensions are set to the
 * double size that is required
 *
 * solution: always use drfftw.h together with double precision transforms!
 */

  fftw_real in[N], out[N], power_spectrum[N/2+1];
  rfftw_plan p;
  int k;
  int m;
  m=5;
  fill(in, N, m);
  p = rfftw_create_plan(N, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
  rfftw_one(p, in, out);
  power_spectrum[0] = out[0]*out[0];  // DC component 
  for (k = 1; k < (N+1)/2; ++k)  // (k < N/2 rounded up) 
       power_spectrum[k] = out[k]*out[k] + out[N-k]*out[N-k];
  if (N % 2 == 0) // N is even 
       power_spectrum[N/2] = out[N/2]*out[N/2];  // Nyquist freq. 
  rfftw_destroy_plan(p);
  printf("\nPower spectrum of %d point signal with %d periods\n", N, m);
  printf("\ni.e. signal frequency %f\n", (double)m/N);
  for (k=0; k<((N+1)/2); ++k)
  {
    printf("  %3.3d, %5f: %10e\n",k, (double)k/N, power_spectrum[k]);
  }
  
  fill(in, N, m);
  p = rfftw_create_plan(N, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
  rfftw_one(p, in, out);
  power_spectrum[0] = out[0]*out[0];  // DC component 
  for (k = 1; k < (N+1)/2; ++k)  // (k < N/2 rounded up) 
       power_spectrum[k] = out[k]*out[k] + out[N-k]*out[N-k];
  if (N % 2 == 0) // N is even 
       power_spectrum[N/2] = out[N/2]*out[N/2];  // Nyquist freq. 
  rfftw_destroy_plan(p);
  printf("\nPower spectrum of %d point signal with %d periods\n", N, m);
  printf("\ni.e. signal frequency %f\n", (double)m/N);
  for (k=0; k<((N+1)/2); ++k)
  {
    printf("  %3.3d, %5f: %10e\n",k, (double)k/N, power_spectrum[k]);
  }
  
  process(200, 20);
  process(50000, 2000);
}

/* ----- END OF cfftwtest.c ----- */