From aadd27aa7e4ebb2e420472b537f8231dbd4568c8 Mon Sep 17 00:00:00 2001 From: Paul Brossier Date: Wed, 7 Nov 2007 16:30:29 +0100 Subject: [PATCH] fft.c,fft.h: refactorised, no more mfft, only one fft --- src/fft.c | 211 +++++++++++++++++++++++++++--------------------------- src/fft.h | 91 ++++++++++++----------- 2 files changed, 153 insertions(+), 149 deletions(-) diff --git a/src/fft.c b/src/fft.c index e5c1a037..955af5cd 100644 --- a/src/fft.c +++ b/src/fft.c @@ -18,7 +18,8 @@ */ #include "aubio_priv.h" -#include "sample.h" +#include "fvec.h" +#include "cvec.h" #include "mathutils.h" #include "fft.h" @@ -40,37 +41,41 @@ #endif struct _aubio_fft_t { - uint_t fft_size; + uint_t winsize; uint_t channels; - real_t *in, *out; - fft_data_t *specdata; + uint_t fft_size; + real_t *in, *out; fftw_plan pfw, pbw; + fft_data_t * specdata; /* complex spectral data */ + fvec_t * compspec; }; -static void aubio_fft_getspectrum(fft_data_t * spectrum, smpl_t *norm, smpl_t * phas, uint_t size); - -aubio_fft_t * new_aubio_fft(uint_t size) { +aubio_fft_t * new_aubio_fft(uint_t winsize, uint_t channels) { aubio_fft_t * s = AUBIO_NEW(aubio_fft_t); + s->winsize = winsize; + s->channels = channels; /* allocate memory */ - s->in = AUBIO_ARRAY(real_t,size); - s->out = AUBIO_ARRAY(real_t,size); + s->in = AUBIO_ARRAY(real_t,winsize); + s->out = AUBIO_ARRAY(real_t,winsize); + s->compspec = new_fvec(winsize,channels); /* create plans */ #ifdef HAVE_COMPLEX_H - s->fft_size = size/2+1; + s->fft_size = winsize/2+1; s->specdata = (fft_data_t*)fftw_malloc(sizeof(fft_data_t)*s->fft_size); - s->pfw = fftw_plan_dft_r2c_1d(size, s->in, s->specdata, FFTW_ESTIMATE); - s->pbw = fftw_plan_dft_c2r_1d(size, s->specdata, s->out, FFTW_ESTIMATE); + s->pfw = fftw_plan_dft_r2c_1d(winsize, s->in, s->specdata, FFTW_ESTIMATE); + s->pbw = fftw_plan_dft_c2r_1d(winsize, s->specdata, s->out, FFTW_ESTIMATE); #else - s->fft_size = size; + s->fft_size = winsize; s->specdata = (fft_data_t*)fftw_malloc(sizeof(fft_data_t)*s->fft_size); - s->pfw = fftw_plan_r2r_1d(size, s->in, s->specdata, FFTW_R2HC, FFTW_ESTIMATE); - s->pbw = fftw_plan_r2r_1d(size, s->specdata, s->out, FFTW_HC2R, FFTW_ESTIMATE); + s->pfw = fftw_plan_r2r_1d(winsize, s->in, s->specdata, FFTW_R2HC, FFTW_ESTIMATE); + s->pbw = fftw_plan_r2r_1d(winsize, s->specdata, s->out, FFTW_HC2R, FFTW_ESTIMATE); #endif return s; } void del_aubio_fft(aubio_fft_t * s) { /* destroy data */ + del_fvec(s->compspec); fftw_destroy_plan(s->pfw); fftw_destroy_plan(s->pbw); fftw_free(s->specdata); @@ -79,117 +84,111 @@ void del_aubio_fft(aubio_fft_t * s) { AUBIO_FREE(s); } -void aubio_fft_do(const aubio_fft_t * s, - const smpl_t * data, fft_data_t * spectrum, const uint_t size) { - uint_t i; - for (i=0;iin[i] = data[i]; - fftw_execute(s->pfw); - for (i=0; i < s->fft_size; i++) spectrum[i] = s->specdata[i]; +void aubio_fft_do(aubio_fft_t * s, fvec_t * input, cvec_t * spectrum) { + aubio_fft_do_complex(s, input, s->compspec); + aubio_fft_get_spectrum(s->compspec, spectrum); } -void aubio_fft_rdo(const aubio_fft_t * s, - const fft_data_t * spectrum, smpl_t * data, const uint_t size) { - uint_t i; - const smpl_t renorm = 1./(smpl_t)size; - for (i=0; i < s->fft_size; i++) s->specdata[i] = spectrum[i]; - fftw_execute(s->pbw); - for (i=0;iout[i]*renorm; +void aubio_fft_rdo(aubio_fft_t * s, cvec_t * spectrum, fvec_t * output) { + aubio_fft_get_realimag(spectrum, s->compspec); + aubio_fft_rdo_complex(s, s->compspec, output); } -#ifdef HAVE_COMPLEX_H - -void aubio_fft_getnorm(smpl_t * norm, fft_data_t * spectrum, uint_t size) { - uint_t i; - for (i=0;ichannels; i++) { + for (j=0; j < s->winsize; j++) { + s->in[j] = input->data[i][j]; + } + fftw_execute(s->pfw); +#if HAVE_COMPLEX_H + compspec->data[i][0] = REAL(s->specdata[0]); + for (j = 1; j < s->fft_size -1 ; j++) { + compspec->data[i][j] = REAL(s->specdata[j]); + compspec->data[i][compspec->length - j] = IMAG(s->specdata[j]); + } + compspec->data[i][s->fft_size-1] = REAL(s->specdata[s->fft_size-1]); +#else + for (j = 0; j < s->fft_size; j++) { + compspec->data[i][j] = s->specdata[j]; + } +#endif } } +void aubio_fft_rdo_complex(aubio_fft_t * s, fvec_t * compspec, fvec_t * output) { + uint_t i, j; + const smpl_t renorm = 1./(smpl_t)s->winsize; + for (i = 0; i < compspec->channels; i++) { +#if HAVE_COMPLEX_H + s->specdata[0] = compspec->data[i][0]; + for (j=1; j < s->fft_size - 1; j++) { + s->specdata[j] = compspec->data[i][j] + + I * compspec->data[i][compspec->length - j]; + } + s->specdata[s->fft_size - 1] = compspec->data[i][s->fft_size - 1]; #else - -void aubio_fft_getnorm(smpl_t * norm, fft_data_t * spectrum, uint_t size) { - uint_t i; - norm[0] = spectrum[0]; - for (i=1;ifft_size; j++) { + s->specdata[j] = compspec->data[i][j]; + } +#endif + fftw_execute(s->pbw); + for (j = 0; j < output->length; j++) { + output->data[i][j] = s->out[j]*renorm; + } + } } -void aubio_fft_getphas(smpl_t * phas, fft_data_t * spectrum, uint_t size) { - uint_t i; - phas[0] = 0; - for (i=1;iwinsize = winsize; - fft->channels = channels; - fft->fft = new_aubio_fft(winsize); - fft->spec = AUBIO_ARRAY(fft_data_t*,channels); - for (i=0; i < channels; i++) - fft->spec[i] = AUBIO_ARRAY(fft_data_t,winsize); - return fft; +void aubio_fft_get_phas(fvec_t * compspec, cvec_t * spectrum) { + uint_t i, j; + for (i = 0; i < spectrum->channels; i++) { + spectrum->phas[i][0] = 0.; + for (j=1; j < spectrum->length - 1; j++) { + spectrum->phas[i][j] = atan2f(compspec->data[i][compspec->length-j], + compspec->data[i][j]); + } + spectrum->phas[i][spectrum->length-1] = 0.; + } } -/* execute stft */ -void aubio_mfft_do (aubio_mfft_t * fft,fvec_t * in,cvec_t * fftgrain){ - uint_t i=0; - /* execute stft */ - for (i=0; i < fft->channels; i++) { - aubio_fft_do (fft->fft,in->data[i],fft->spec[i],fft->winsize); - /* put norm and phase into fftgrain */ - aubio_fft_getnorm(fftgrain->norm[i], fft->spec[i], fft->winsize); - aubio_fft_getphas(fftgrain->phas[i], fft->spec[i], fft->winsize); +void aubio_fft_get_norm(fvec_t * compspec, cvec_t * spectrum) { + uint_t i, j = 0; + for (i = 0; i < spectrum->channels; i++) { + spectrum->norm[i][0] = compspec->data[i][0]; + for (j=1; j < spectrum->length - 1; j++) { + spectrum->norm[i][j] = SQRT(SQR(compspec->data[i][j]) + + SQR(compspec->data[i][compspec->length - j]) ); + } + spectrum->norm[i][spectrum->length-1] = compspec->data[i][compspec->length/2]; } } -/* execute inverse fourier transform */ -void aubio_mfft_rdo(aubio_mfft_t * fft,cvec_t * fftgrain, fvec_t * out){ - uint_t i=0; - for (i=0; i < fft->channels; i++) { - aubio_fft_getspectrum(fft->spec[i],fftgrain->norm[i],fftgrain->phas[i],fft->winsize); - aubio_fft_rdo(fft->fft,fft->spec[i],out->data[i],fft->winsize); +void aubio_fft_get_imag(cvec_t * spectrum, fvec_t * compspec) { + uint_t i, j; + for (i = 0; i < compspec->channels; i++) { + for (j = 1; j < compspec->length / 2 + 1; j++) { + compspec->data[i][compspec->length - j] = + spectrum->norm[i][j]*SIN(spectrum->phas[i][j]); + } } } -void del_aubio_mfft(aubio_mfft_t * fft) { - uint_t i; - for (i=0; i < fft->channels; i++) - AUBIO_FREE(fft->spec[i]); - AUBIO_FREE(fft->spec); - del_aubio_fft(fft->fft); - AUBIO_FREE(fft); +void aubio_fft_get_real(cvec_t * spectrum, fvec_t * compspec) { + uint_t i, j; + for (i = 0; i < compspec->channels; i++) { + for (j = 0; j< compspec->length / 2 + 1; j++) { + compspec->data[i][j] = + spectrum->norm[i][j]*COS(spectrum->phas[i][j]); + } + } } diff --git a/src/fft.h b/src/fft.h index 72187eb3..1737bd95 100644 --- a/src/fft.h +++ b/src/fft.h @@ -66,93 +66,98 @@ typedef struct _aubio_fft_t aubio_fft_t; /** create new FFT computation object \param size length of the FFT + \param channels number of channels */ -aubio_fft_t * new_aubio_fft(uint_t size); +aubio_fft_t * new_aubio_fft(uint_t size, uint_t channels); /** delete FFT object \param s fft object as returned by new_aubio_fft */ void del_aubio_fft(aubio_fft_t * s); + /** compute forward FFT \param s fft object as returned by new_aubio_fft - \param data input signal + \param input input signal \param spectrum output spectrum - \param size length of the input vector */ -void aubio_fft_do (const aubio_fft_t *s, const smpl_t * data, - fft_data_t * spectrum, const uint_t size); +void aubio_fft_do (aubio_fft_t *s, fvec_t * input, cvec_t * spectrum); /** compute backward (inverse) FFT \param s fft object as returned by new_aubio_fft \param spectrum input spectrum - \param data output signal - \param size length of the input vector + \param output output signal */ -void aubio_fft_rdo(const aubio_fft_t *s, const fft_data_t * spectrum, - smpl_t * data, const uint_t size); -/** compute norm vector from input spectrum +void aubio_fft_rdo (aubio_fft_t *s, cvec_t * spectrum, fvec_t * output); - \param norm magnitude vector output - \param spectrum spectral data input - \param size size of the vectors +/** compute forward FFT + + \param s fft object as returned by new_aubio_fft + \param input real input signal + \param compspec complex output fft real/imag */ -void aubio_fft_getnorm(smpl_t * norm, fft_data_t * spectrum, uint_t size); -/** compute phase vector from input spectrum - - \param phase phase vector output - \param spectrum spectral data input - \param size size of the vectors +void aubio_fft_do_complex (aubio_fft_t *s, fvec_t * input, fvec_t * compspec); +/** compute backward (inverse) FFT from real/imag + + \param s fft object as returned by new_aubio_fft + \param compspec real/imag input fft array + \param output real output array */ -void aubio_fft_getphas(smpl_t * phase, fft_data_t * spectrum, uint_t size); +void aubio_fft_rdo_complex (aubio_fft_t *s, fvec_t * compspec, fvec_t * output); -/** FFT object (using cvec) +/** convert real/imag spectrum to norm/phas spectrum - This object works similarly as aubio_fft_t, except the spectral data is - stored in a cvec_t as two vectors, magnitude and phase. + \param compspec real/imag input fft array + \param spectrum cvec norm/phas output array */ -typedef struct _aubio_mfft_t aubio_mfft_t; - -/** create new FFT computation object +void aubio_fft_get_spectrum(fvec_t * compspec, cvec_t * spectrum); +/** convert real/imag spectrum to norm/phas spectrum - \param winsize length of the FFT - \param channels number of channels + \param compspec real/imag input fft array + \param spectrum cvec norm/phas output array */ -aubio_mfft_t * new_aubio_mfft(uint_t winsize, uint_t channels); -/** compute forward FFT +void aubio_fft_get_realimag(cvec_t * spectrum, fvec_t * compspec); + +/** compute phas spectrum from real/imag parts - \param fft fft object as returned by new_aubio_mfft - \param in input signal - \param fftgrain output spectrum + \param compspec real/imag input fft array + \param spectrum cvec norm/phas output array */ -void aubio_mfft_do (aubio_mfft_t * fft,fvec_t * in,cvec_t * fftgrain); -/** compute backward (inverse) FFT +void aubio_fft_get_phas(fvec_t * compspec, cvec_t * spectrum); +/** compute imaginary part from the norm/phas cvec - \param fft fft object as returned by new_aubio_mfft - \param fftgrain input spectrum (cvec) - \param out output signal + \param spectrum norm/phas input array + \param compspec real/imag output fft array */ -void aubio_mfft_rdo(aubio_mfft_t * fft,cvec_t * fftgrain, fvec_t * out); -/** delete FFT object +void aubio_fft_get_imag(cvec_t * spectrum, fvec_t * compspec); - \param fft fft object as returned by new_aubio_mfft +/** compute norm component from real/imag parts + + \param compspec real/imag input fft array + \param spectrum cvec norm/phas output array */ -void del_aubio_mfft(aubio_mfft_t * fft); +void aubio_fft_get_norm(fvec_t * compspec, cvec_t * spectrum); +/** compute real part from norm/phas components + + \param spectrum norm/phas input array + \param compspec real/imag output fft array +*/ +void aubio_fft_get_real(cvec_t * spectrum, fvec_t * compspec); #ifdef __cplusplus } #endif -#endif +#endif // FFT_H_ -- 2.26.2