src/fft.c

   1 /*
   2    Copyright (C) 2003 Paul Brossier
   3
   4    This program is free software; you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation; either version 2 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program; if not, write to the Free Software
  16    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  17
  18 */
  19
  20 #include "aubio_priv.h"
  21 #include "sample.h"
  22 #include "mathutils.h"
  23 #include "fft.h"
  24
  25 #if FFTW3F_SUPPORT
  26 #define fftw_malloc             fftwf_malloc
  27 #define fftw_free               fftwf_free
  28 #define fftw_execute            fftwf_execute
  29 #define fftw_plan_dft_r2c_1d    fftwf_plan_dft_r2c_1d
  30 #define fftw_plan_dft_c2r_1d    fftwf_plan_dft_c2r_1d
  31 #define fftw_plan               fftwf_plan
  32 #define fftw_destroy_plan       fftwf_destroy_plan
  33 #endif
  34
  35 #if FFTW3F_SUPPORT
  36 #define real_t smpl_t
  37 #else
  38 #define real_t lsmp_t
  39 #endif
  40
  41 struct _aubio_fft_t {
  42         uint_t fft_size;
  43         uint_t channels;
  44         real_t          *in, *out;
  45         fft_data_t      *specdata;
  46         fftw_plan       pfw, pbw;
  47 };
  48
  49 static aubio_fft_t * aubio_fft_alloc(uint_t size);
  50 static void aubio_fft_free(aubio_fft_t *s);
  51
  52 static aubio_fft_t * aubio_fft_alloc(uint_t size) {
  53         aubio_fft_t * s = AUBIO_NEW(aubio_fft_t);
  54         /* allocate memory */
  55         s->in       = AUBIO_ARRAY(real_t,size);
  56         s->out      = AUBIO_ARRAY(real_t,size);
  57         s->specdata = (fft_data_t*)fftw_malloc(sizeof(fft_data_t)*size);
  58         return s;
  59 }
  60
  61 static void aubio_fft_free(aubio_fft_t * s) {
  62         /* destroy data */
  63         fftw_destroy_plan(s->pfw);
  64         fftw_destroy_plan(s->pbw);
  65         if (s->specdata)        fftw_free(s->specdata);
  66         if (s->out)             AUBIO_FREE(s->out);
  67         if (s->in )             AUBIO_FREE(s->in );
  68 }
  69
  70 aubio_fft_t * new_aubio_fft(uint_t size) {
  71         aubio_fft_t * s =(aubio_fft_t *)aubio_fft_alloc(size);
  72         /* create plans */
  73         s->pfw = fftw_plan_dft_r2c_1d(size, s->in,  s->specdata, FFTW_ESTIMATE);
  74         s->pbw = fftw_plan_dft_c2r_1d(size, s->specdata, s->out, FFTW_ESTIMATE);
  75         return s;
  76 }
  77
  78 void del_aubio_fft(aubio_fft_t * s) {
  79         aubio_fft_free(s);
  80         AUBIO_FREE(s);
  81 }
  82
  83 void aubio_fft_do(const aubio_fft_t * s,
  84                 const smpl_t * data, fft_data_t * spectrum,
  85                 const uint_t size) {
  86         uint_t i;
  87         for (i=0;i<size;i++) s->in[i] = data[i];
  88         fftw_execute(s->pfw);
  89         for (i=0;i<size;i++) spectrum[i] = s->specdata[i];
  90 }
  91
  92 void aubio_fft_rdo(const aubio_fft_t * s,
  93                 const fft_data_t * spectrum,
  94                 smpl_t * data,
  95                 const uint_t size) {
  96         uint_t i;
  97         const smpl_t renorm = 1./(smpl_t)size;
  98         for (i=0;i<size;i++) s->specdata[i] = spectrum[i];
  99         fftw_execute(s->pbw);
 100         for (i=0;i<size;i++) data[i] = s->out[i]*renorm;
 101 }
 102
 103
 104 void aubio_fft_getnorm(smpl_t * norm, fft_data_t * spectrum, uint_t size) {
 105         uint_t i;
 106         for (i=0;i<size;i++) norm[i] = ABSC(spectrum[i]);
 107 }
 108
 109 void aubio_fft_getphas(smpl_t * phas, fft_data_t * spectrum, uint_t size) {
 110         uint_t i;
 111         for (i=0;i<size;i++) phas[i] = ARGC(spectrum[i]);
 112 }
 113
 114
 115 /* new interface aubio_mfft */
 116 struct _aubio_mfft_t {
 117         aubio_fft_t * fft;      /* fftw interface */
 118         fft_data_t ** spec;     /* complex spectral data */
 119         uint_t winsize;
 120         uint_t channels;
 121 };
 122
 123 aubio_mfft_t * new_aubio_mfft(uint_t winsize, uint_t channels){
 124         uint_t i;
 125         aubio_mfft_t * fft = AUBIO_NEW(aubio_mfft_t);
 126         fft->winsize       = winsize;
 127         fft->channels      = channels;
 128         fft->fft           = new_aubio_fft(winsize);
 129         fft->spec          = AUBIO_ARRAY(fft_data_t*,channels);
 130         for (i=0; i < channels; i++)
 131                 fft->spec[i] = AUBIO_ARRAY(fft_data_t,winsize);
 132         return fft;
 133 }
 134
 135 /* execute stft */
 136 void aubio_mfft_do (aubio_mfft_t * fft,fvec_t * in,cvec_t * fftgrain){
 137         uint_t i=0;
 138         /* execute stft */
 139         for (i=0; i < fft->channels; i++) {
 140                 aubio_fft_do (fft->fft,in->data[i],fft->spec[i],fft->winsize);
 141                 /* put norm and phase into fftgrain */
 142                 aubio_fft_getnorm(fftgrain->norm[i], fft->spec[i], fft->winsize/2+1);
 143                 aubio_fft_getphas(fftgrain->phas[i], fft->spec[i], fft->winsize/2+1);
 144         }
 145 }
 146
 147 /* execute inverse fourier transform */
 148 void aubio_mfft_rdo(aubio_mfft_t * fft,cvec_t * fftgrain, fvec_t * out){
 149         uint_t i=0,j;
 150         for (i=0; i < fft->channels; i++) {
 151                 for (j=0; j<fft->winsize/2+1; j++) {
 152                         fft->spec[i][j]  = CEXPC(I*aubio_unwrap2pi(fftgrain->phas[i][j]));
 153                         fft->spec[i][j] *= fftgrain->norm[i][j];
 154                 }
 155                 aubio_fft_rdo(fft->fft,fft->spec[i],out->data[i],fft->winsize);
 156         }
 157 }
 158
 159 void del_aubio_mfft(aubio_mfft_t * fft) {
 160         uint_t i;
 161         for (i=0; i < fft->channels; i++)
 162                 AUBIO_FREE(fft->spec[i]);
 163         AUBIO_FREE(fft->spec);
 164         aubio_fft_free(fft->fft);
 165         AUBIO_FREE(fft);
 166 }