Changeset d57d1de for src/fft.c


Ignore:
Timestamp:
Nov 7, 2007, 4:51:42 PM (16 years ago)
Author:
Paul Brossier <piem@piem.org>
Branches:
feature/autosink, feature/cnn, feature/cnn_org, feature/constantq, feature/crepe, feature/crepe_org, feature/pitchshift, feature/pydocstrings, feature/timestretch, fix/ffmpeg5, master, pitchshift, sampler, timestretch, yinfft+
Children:
4368223
Parents:
e00f769 (diff), 038852a (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

update fft.py tests, merge from banane

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/fft.c

    re00f769 rd57d1de  
    1919
    2020#include "aubio_priv.h"
    21 #include "sample.h"
     21#include "fvec.h"
     22#include "cvec.h"
    2223#include "mathutils.h"
    2324#include "fft.h"
    … …  
    4142
    4243struct _aubio_fft_t {
     44  uint_t winsize;
     45  uint_t channels;
    4346  uint_t fft_size;
    44   uint_t channels;
    45   real_t    *in, *out;
    46   fft_data_t   *specdata;
     47  real_t *in, *out;
    4748  fftw_plan   pfw, pbw;
     49  fft_data_t * specdata;     /* complex spectral data */
     50  fvec_t * compspec;
    4851};
    4952
    50 static void aubio_fft_getspectrum(fft_data_t * spectrum, smpl_t *norm, smpl_t * phas, uint_t size);
    51 
    52 aubio_fft_t * new_aubio_fft(uint_t size) {
     53aubio_fft_t * new_aubio_fft(uint_t winsize, uint_t channels) {
    5354  aubio_fft_t * s = AUBIO_NEW(aubio_fft_t);
     55  s->winsize  = winsize;
     56  s->channels = channels;
    5457  /* allocate memory */
    55   s->in       = AUBIO_ARRAY(real_t,size);
    56   s->out      = AUBIO_ARRAY(real_t,size);
     58  s->in       = AUBIO_ARRAY(real_t,winsize);
     59  s->out      = AUBIO_ARRAY(real_t,winsize);
     60  s->compspec = new_fvec(winsize,channels);
    5761  /* create plans */
    5862#ifdef HAVE_COMPLEX_H
    59   s->fft_size = size/2+1;
     63  s->fft_size = winsize/2+1;
    6064  s->specdata = (fft_data_t*)fftw_malloc(sizeof(fft_data_t)*s->fft_size);
    61   s->pfw = fftw_plan_dft_r2c_1d(size, s->in,  s->specdata, FFTW_ESTIMATE);
    62   s->pbw = fftw_plan_dft_c2r_1d(size, s->specdata, s->out, FFTW_ESTIMATE);
     65  s->pfw = fftw_plan_dft_r2c_1d(winsize, s->in,  s->specdata, FFTW_ESTIMATE);
     66  s->pbw = fftw_plan_dft_c2r_1d(winsize, s->specdata, s->out, FFTW_ESTIMATE);
    6367#else
    64   s->fft_size = size;
     68  s->fft_size = winsize;
    6569  s->specdata = (fft_data_t*)fftw_malloc(sizeof(fft_data_t)*s->fft_size);
    66   s->pfw = fftw_plan_r2r_1d(size, s->in,  s->specdata, FFTW_R2HC, FFTW_ESTIMATE);
    67   s->pbw = fftw_plan_r2r_1d(size, s->specdata, s->out, FFTW_HC2R, FFTW_ESTIMATE);
     70  s->pfw = fftw_plan_r2r_1d(winsize, s->in,  s->specdata, FFTW_R2HC, FFTW_ESTIMATE);
     71  s->pbw = fftw_plan_r2r_1d(winsize, s->specdata, s->out, FFTW_HC2R, FFTW_ESTIMATE);
    6872#endif
    6973  return s;
    … …  
    7276void del_aubio_fft(aubio_fft_t * s) {
    7377  /* destroy data */
     78  del_fvec(s->compspec);
    7479  fftw_destroy_plan(s->pfw);
    7580  fftw_destroy_plan(s->pbw);
    … …  
    8085}
    8186
    82 void aubio_fft_do(const aubio_fft_t * s,
    83     const smpl_t * data, fft_data_t * spectrum, const uint_t size) {
    84   uint_t i;
    85   for (i=0;i<size;i++) s->in[i] = data[i];
    86   fftw_execute(s->pfw);
    87   for (i=0; i < s->fft_size; i++) spectrum[i] = s->specdata[i];
     87void aubio_fft_do(aubio_fft_t * s, fvec_t * input, cvec_t * spectrum) {
     88  aubio_fft_do_complex(s, input, s->compspec);
     89  aubio_fft_get_spectrum(s->compspec, spectrum);
    8890}
    8991
    90 void aubio_fft_rdo(const aubio_fft_t * s,
    91     const fft_data_t * spectrum, smpl_t * data, const uint_t size) {
    92   uint_t i;
    93   const smpl_t renorm = 1./(smpl_t)size;
    94   for (i=0; i < s->fft_size; i++) s->specdata[i] = spectrum[i];
    95   fftw_execute(s->pbw);
    96   for (i=0;i<size;i++) data[i] = s->out[i]*renorm;
     92void aubio_fft_rdo(aubio_fft_t * s, cvec_t * spectrum, fvec_t * output) {
     93  aubio_fft_get_realimag(spectrum, s->compspec);
     94  aubio_fft_rdo_complex(s, s->compspec, output);
    9795}
    9896
    99 #ifdef HAVE_COMPLEX_H
    100 
    101 void aubio_fft_getnorm(smpl_t * norm, fft_data_t * spectrum, uint_t size) {
    102   uint_t i;
    103   for (i=0;i<size/2+1;i++) norm[i] = ABSC(spectrum[i]);
    104 }
    105 
    106 void aubio_fft_getphas(smpl_t * phas, fft_data_t * spectrum, uint_t size) {
    107   uint_t i;
    108   for (i=0;i<size/2+1;i++) phas[i] = ARGC(spectrum[i]);
    109 }
    110 
    111 void aubio_fft_getspectrum(fft_data_t * spectrum, smpl_t *norm, smpl_t * phas, uint_t size) {
    112   uint_t j;
    113   for (j=0; j<size/2+1; j++) {
    114     spectrum[j]  = CEXPC(I*phas[j]);
    115     spectrum[j] *= norm[j];
     97void aubio_fft_do_complex(aubio_fft_t * s, fvec_t * input, fvec_t * compspec) {
     98  uint_t i, j;
     99  for (i = 0; i < s->channels; i++) {
     100    for (j=0; j < s->winsize; j++) {
     101      s->in[j] = input->data[i][j];
     102    }
     103    fftw_execute(s->pfw);
     104#if HAVE_COMPLEX_H
     105    compspec->data[i][0] = REAL(s->specdata[0]);
     106    for (j = 1; j < s->fft_size -1 ; j++) {
     107      compspec->data[i][j] = REAL(s->specdata[j]);
     108      compspec->data[i][compspec->length - j] = IMAG(s->specdata[j]);
     109    }
     110    compspec->data[i][s->fft_size-1] = REAL(s->specdata[s->fft_size-1]);
     111#else
     112    for (j = 0; j < s->fft_size; j++) {
     113      compspec->data[i][j] = s->specdata[j];
     114    }
     115#endif
    116116  }
    117117}
    118118
     119void aubio_fft_rdo_complex(aubio_fft_t * s, fvec_t * compspec, fvec_t * output) {
     120  uint_t i, j;
     121  const smpl_t renorm = 1./(smpl_t)s->winsize;
     122  for (i = 0; i < compspec->channels; i++) {
     123#if HAVE_COMPLEX_H
     124    s->specdata[0] = compspec->data[i][0];
     125    for (j=1; j < s->fft_size - 1; j++) {
     126      s->specdata[j] = compspec->data[i][j] +
     127        I * compspec->data[i][compspec->length - j];
     128    }
     129    s->specdata[s->fft_size - 1] = compspec->data[i][s->fft_size - 1];
    119130#else
    120 
    121 void aubio_fft_getnorm(smpl_t * norm, fft_data_t * spectrum, uint_t size) {
    122   uint_t i;
    123   norm[0] = spectrum[0];
    124   for (i=1;i<size/2;i++) norm[i] = SQRT((SQR(spectrum[i]) + SQR(spectrum[size-i])));
    125   norm[size/2] = spectrum[size/2];
    126 }
    127 
    128 void aubio_fft_getphas(smpl_t * phas, fft_data_t * spectrum, uint_t size) {
    129   uint_t i;
    130   phas[0] = 0;
    131   for (i=1;i<size/2+1;i++) phas[i] = atan2f(spectrum[size-i] , spectrum[i]);
    132   phas[size/2] = 0;
    133 }
    134 
    135 void aubio_fft_getspectrum(fft_data_t * spectrum, smpl_t *norm, smpl_t * phas, uint_t size) {
    136   uint_t j;
    137   for (j=0; j<size/2+1; j++) {
    138     spectrum[j]       = norm[j]*COS(phas[j]);
    139   }
    140   for (j=1; j<size/2+1; j++) {
    141     spectrum[size-j]  = norm[j]*SIN(phas[j]);
     131    for (j=0; j < s->fft_size; j++) {
     132      s->specdata[j] = compspec->data[i][j];
     133    }
     134#endif
     135    fftw_execute(s->pbw);
     136    for (j = 0; j < output->length; j++) {
     137      output->data[i][j] = s->out[j]*renorm;
     138    }
    142139  }
    143140}
    144141
    145 #endif
    146 
    147 /* new interface aubio_mfft */
    148 struct _aubio_mfft_t {
    149         aubio_fft_t * fft;      /* fftw interface */
    150         fft_data_t ** spec;     /* complex spectral data */
    151         uint_t winsize;
    152         uint_t channels;
    153 };
    154 
    155 aubio_mfft_t * new_aubio_mfft(uint_t winsize, uint_t channels){
    156   uint_t i;
    157   aubio_mfft_t * fft = AUBIO_NEW(aubio_mfft_t);
    158   fft->winsize       = winsize;
    159   fft->channels      = channels;
    160   fft->fft           = new_aubio_fft(winsize);
    161   fft->spec          = AUBIO_ARRAY(fft_data_t*,channels);
    162   for (i=0; i < channels; i++)
    163     fft->spec[i] = AUBIO_ARRAY(fft_data_t,winsize);
    164   return fft;
     142void aubio_fft_get_spectrum(fvec_t * compspec, cvec_t * spectrum) {
     143  aubio_fft_get_phas(compspec, spectrum);
     144  aubio_fft_get_norm(compspec, spectrum);
    165145}
    166146
    167 /* execute stft */
    168 void aubio_mfft_do (aubio_mfft_t * fft,fvec_t * in,cvec_t * fftgrain){
    169   uint_t i=0;
    170   /* execute stft */
    171   for (i=0; i < fft->channels; i++) {
    172     aubio_fft_do (fft->fft,in->data[i],fft->spec[i],fft->winsize);
    173     /* put norm and phase into fftgrain */
    174     aubio_fft_getnorm(fftgrain->norm[i], fft->spec[i], fft->winsize);
    175     aubio_fft_getphas(fftgrain->phas[i], fft->spec[i], fft->winsize);
     147void aubio_fft_get_realimag(cvec_t * spectrum, fvec_t * compspec) {
     148  aubio_fft_get_imag(spectrum, compspec);
     149  aubio_fft_get_real(spectrum, compspec);
     150}
     151
     152void aubio_fft_get_phas(fvec_t * compspec, cvec_t * spectrum) {
     153  uint_t i, j;
     154  for (i = 0; i < spectrum->channels; i++) {
     155    spectrum->phas[i][0] = 0.;
     156    for (j=1; j < spectrum->length - 1; j++) {
     157      spectrum->phas[i][j] = atan2f(compspec->data[i][compspec->length-j],
     158          compspec->data[i][j]);
     159    }
     160    spectrum->phas[i][spectrum->length-1] = 0.;
    176161  }
    177162}
    178163
    179 /* execute inverse fourier transform */
    180 void aubio_mfft_rdo(aubio_mfft_t * fft,cvec_t * fftgrain, fvec_t * out){
    181   uint_t i=0;
    182   for (i=0; i < fft->channels; i++) {
    183     aubio_fft_getspectrum(fft->spec[i],fftgrain->norm[i],fftgrain->phas[i],fft->winsize);
    184     aubio_fft_rdo(fft->fft,fft->spec[i],out->data[i],fft->winsize);
     164void aubio_fft_get_norm(fvec_t * compspec, cvec_t * spectrum) {
     165  uint_t i, j = 0;
     166  for (i = 0; i < spectrum->channels; i++) {
     167    spectrum->norm[i][0] = compspec->data[i][0];
     168    for (j=1; j < spectrum->length - 1; j++) {
     169      spectrum->norm[i][j] = SQRT(SQR(compspec->data[i][j])
     170          + SQR(compspec->data[i][compspec->length - j]) );
     171    }
     172    spectrum->norm[i][spectrum->length-1] = compspec->data[i][compspec->length/2];
    185173  }
    186174}
    187175
    188 void del_aubio_mfft(aubio_mfft_t * fft) {
    189   uint_t i;
    190   for (i=0; i < fft->channels; i++)
    191     AUBIO_FREE(fft->spec[i]);
    192   AUBIO_FREE(fft->spec);
    193   del_aubio_fft(fft->fft);
    194   AUBIO_FREE(fft);       
     176void aubio_fft_get_imag(cvec_t * spectrum, fvec_t * compspec) {
     177  uint_t i, j;
     178  for (i = 0; i < compspec->channels; i++) {
     179    for (j = 1; j < compspec->length / 2 + 1; j++) {
     180      compspec->data[i][compspec->length - j] =
     181        spectrum->norm[i][j]*SIN(spectrum->phas[i][j]);
     182    }
     183  }
    195184}
     185
     186void aubio_fft_get_real(cvec_t * spectrum, fvec_t * compspec) {
     187  uint_t i, j;
     188  for (i = 0; i < compspec->channels; i++) {
     189    for (j = 0; j< compspec->length / 2 + 1; j++) {
     190      compspec->data[i][j] =
     191        spectrum->norm[i][j]*COS(spectrum->phas[i][j]);
     192    }
     193  }
     194}
Note: See TracChangeset for help on using the changeset viewer.