fft.c

   1 /*
   2   ZynAddSubFX - a software synthesizer
   3
   4   FFTwrapper.c  -  A wrapper for Fast Fourier Transforms
   5   Copyright (C) 2002-2005 Nasca Octavian Paul
   6   Author: Nasca Octavian Paul
   7
   8   This program is free software; you can redistribute it and/or modify
   9   it under the terms of version 2 of the GNU General Public License
  10   as published by the Free Software Foundation.
  11
  12   This program is distributed in the hope that it will be useful,
  13   but WITHOUT ANY WARRANTY; without even the implied warranty of
  14   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15   GNU General Public License (version 2) for more details.
  16
  17   You should have received a copy of the GNU General Public License (version 2)
  18   along with this program; if not, write to the Free Software Foundation,
  19   Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  20
  21 */
  22
  23 #include <math.h>
  24 #include <stdlib.h>
  25
  26 #ifdef FFTW_VERSION_2
  27 # include <fftw.h>
  28 /* If you got error messages about rfftw.h, replace the next include  line with "#include <srfftw.h>"
  29    or with "#include <drfftw.h> (if one doesn't work try the other). It may be necessary to replace
  30    the <fftw.h> with <dfftw.h> or <sfftw.h>. If the neither one doesn't work,
  31    please install latest version of fftw(recomanded from the sources) from www.fftw.org.
  32    If you'll install fftw3 you need to change the Makefile.inc
  33    Hope all goes right." */
  34 # include <rfftw.h>
  35 #else
  36 # include <fftw3.h>
  37 # define fftw_real double
  38 # define rfftw_plan fftw_plan
  39 #endif
  40
  41 #include "fft.h"
  42
  43 struct zyn_fft
  44 {
  45   int size;
  46   fftw_real * tmp_data1;
  47   fftw_real * tmp_data2;
  48   rfftw_plan plan;
  49   rfftw_plan plan_inv;
  50 };
  51
  52 zyn_fft_handle
  53 zyn_fft_create(
  54   int fftsize)
  55 {
  56   struct zyn_fft * fft_ptr;
  57
  58   fft_ptr = malloc(sizeof(struct zyn_fft));
  59
  60   fft_ptr->size = fftsize;
  61
  62   fft_ptr->tmp_data1 = malloc(sizeof(fftw_real) * fftsize);
  63   fft_ptr->tmp_data2 = malloc(sizeof(fftw_real) * fftsize);
  64
  65 #ifdef FFTW_VERSION_2
  66   fft_ptr->plan = rfftw_create_plan(fftsize, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE | FFTW_IN_PLACE);
  67   fft_ptr->plan_inv = rfftw_create_plan(fftsize, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE | FFTW_IN_PLACE);
  68 #else
  69   fft_ptr->plan = fftw_plan_r2r_1d(fftsize, fft_ptr->tmp_data1, fft_ptr->tmp_data1, FFTW_R2HC, FFTW_ESTIMATE);
  70   fft_ptr->plan_inv = fftw_plan_r2r_1d(fftsize, fft_ptr->tmp_data2, fft_ptr->tmp_data2, FFTW_HC2R, FFTW_ESTIMATE);
  71 #endif
  72
  73   return (zyn_fft_handle)fft_ptr;
  74 }
  75
  76 #define fft_ptr ((struct zyn_fft *)handle)
  77
  78 void
  79 zyn_fft_destroy(
  80   zyn_fft_handle handle)
  81 {
  82 #ifdef FFTW_VERSION_2
  83   rfftw_destroy_plan(fft_ptr->plan);
  84   rfftw_destroy_plan(fft_ptr->plan_inv);
  85 #else
  86   fftw_destroy_plan(fft_ptr->plan);
  87   fftw_destroy_plan(fft_ptr->plan_inv);
  88 #endif
  89
  90   free(fft_ptr->tmp_data1);
  91   free(fft_ptr->tmp_data2);
  92 }
  93
  94 /*
  95  * do the Fast Fourier Transform
  96  */
  97 void
  98 zyn_fft_smps2freqs(
  99   zyn_fft_handle handle,
 100   REALTYPE * smps,
 101   struct zyn_fft_freqs * freqs_ptr)
 102 {
 103   int i;
 104   fftw_real * tmp_data_ptr;
 105
 106   for (i = 0 ; i < fft_ptr->size ; i++)
 107   {
 108     fft_ptr->tmp_data1[i] = smps[i];
 109   }
 110
 111 #ifdef FFTW_VERSION_2
 112   rfftw_one(fft_ptr->plan, fft_ptr->tmp_data1, fft_ptr->tmp_data2);
 113   tmp_data_ptr = fft_ptr->tmp_data2;
 114 #else
 115   fftw_execute(fft_ptr->plan);
 116   tmp_data_ptr = fft_ptr->tmp_data1;
 117 #endif
 118
 119   for (i = 0 ; i < fft_ptr->size / 2 ; i++)
 120   {
 121     freqs_ptr->c[i] = tmp_data_ptr[i];
 122     if (i != 0)
 123     {
 124       freqs_ptr->s[i] = tmp_data_ptr[fft_ptr->size - i];
 125     }
 126   }
 127
 128   fft_ptr->tmp_data2[fft_ptr->size / 2] = 0.0;
 129 }
 130
 131 /*
 132  * do the Inverse Fast Fourier Transform
 133  */
 134 void
 135 zyn_fft_freqs2smps(
 136   zyn_fft_handle handle,
 137   struct zyn_fft_freqs * freqs_ptr,
 138   REALTYPE * smps)
 139 {
 140   int i;
 141   fftw_real * tmp_data_ptr;
 142
 143   fft_ptr->tmp_data2[fft_ptr->size / 2] = 0.0;
 144
 145 #ifdef FFTW_VERSION_2
 146   tmp_data_ptr = fft_ptr->tmp_data1;
 147 #else
 148   tmp_data_ptr = fft_ptr->tmp_data2;
 149 #endif
 150
 151   for (i = 0 ; i < fft_ptr->size / 2 ; i++)
 152   {
 153     tmp_data_ptr[i] = freqs_ptr->c[i];
 154     if (i != 0)
 155     {
 156       tmp_data_ptr[fft_ptr->size - i] = freqs_ptr->s[i];
 157     }
 158   }
 159
 160 #ifdef FFTW_VERSION_2
 161   rfftw_one(fft_ptr->plan_inv, fft_ptr->tmp_data1, fft_ptr->tmp_data2);
 162 #else
 163   fftw_execute(fft_ptr->plan_inv);
 164 #endif
 165
 166   for (i = 0 ; i < fft_ptr->size ; i++)
 167   {
 168     smps[i] = fft_ptr->tmp_data2[i];
 169   }
 170 }
 171
 172 void
 173 zyn_fft_freqs_init(
 174   struct zyn_fft_freqs * f,
 175   int size)
 176 {
 177   f->c = malloc(size * sizeof(zyn_sample_type));
 178   f->s = malloc(size * sizeof(zyn_sample_type));
 179
 180   silence_two_buffers(f->c, f->s, size);
 181 }
 182
 183 void
 184 zyn_fft_freqs_uninit(
 185   struct zyn_fft_freqs * f)
 186 {
 187   free(f->c);
 188   free(f->s);
 189
 190   f->c = NULL;
 191   f->s = NULL;
 192 }