Changes in src/fft.c [237f632:4b6937b]
Legend:
- Unmodified
- Added
- Removed
-
src/fft.c
r237f632 r4b6937b 19 19 20 20 #include "aubio_priv.h" 21 #include "sample.h" 21 #include "fvec.h" 22 #include "cvec.h" 22 23 #include "mathutils.h" 23 24 #include "fft.h" 24 25 25 26 #if FFTW3F_SUPPORT 26 #define fftw_malloc 27 #define fftw_free 28 #define fftw_execute 29 #define fftw_plan_dft_r2c_1d 30 #define fftw_plan_dft_c2r_1d 31 #define fftw_plan_r2r_1d fftwf_plan_r2r_1d32 #define fftw_plan 33 #define fftw_destroy_plan 27 #define fftw_malloc fftwf_malloc 28 #define fftw_free fftwf_free 29 #define fftw_execute fftwf_execute 30 #define fftw_plan_dft_r2c_1d fftwf_plan_dft_r2c_1d 31 #define fftw_plan_dft_c2r_1d fftwf_plan_dft_c2r_1d 32 #define fftw_plan_r2r_1d fftwf_plan_r2r_1d 33 #define fftw_plan fftwf_plan 34 #define fftw_destroy_plan fftwf_destroy_plan 34 35 #endif 35 36 … … 41 42 42 43 struct _aubio_fft_t { 43 uint_t fft_size; 44 uint_t channels; 45 real_t *in, *out; 46 fft_data_t *specdata; 47 fftw_plan pfw, pbw; 44 uint_t winsize; 45 uint_t channels; 46 uint_t fft_size; 47 real_t *in, *out; 48 fftw_plan pfw, pbw; 49 fft_data_t * specdata; /* complex spectral data */ 50 fvec_t * compspec; 48 51 }; 49 52 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) { 53 aubio_fft_t * s = AUBIO_NEW(aubio_fft_t);54 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 53 aubio_fft_t * new_aubio_fft(uint_t winsize, uint_t channels) { 54 aubio_fft_t * s = AUBIO_NEW(aubio_fft_t); 55 s->winsize = winsize; 56 s->channels = channels; 57 /* allocate memory */ 58 s->in = AUBIO_ARRAY(real_t,winsize); 59 s->out = AUBIO_ARRAY(real_t,winsize); 60 s->compspec = new_fvec(winsize,channels); 61 /* create plans */ 59 62 #ifdef HAVE_COMPLEX_H 60 s->pfw = fftw_plan_dft_r2c_1d(size, s->in, s->specdata, FFTW_ESTIMATE); 61 s->pbw = fftw_plan_dft_c2r_1d(size, s->specdata, s->out, FFTW_ESTIMATE); 63 s->fft_size = winsize/2 + 1; 64 s->specdata = (fft_data_t*)fftw_malloc(sizeof(fft_data_t)*s->fft_size); 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); 62 67 #else 63 s->pfw = fftw_plan_r2r_1d(size, s->in, s->specdata, FFTW_R2HC, FFTW_ESTIMATE); 64 s->pbw = fftw_plan_r2r_1d(size, s->specdata, s->out, FFTW_HC2R, FFTW_ESTIMATE); 68 s->fft_size = winsize; 69 s->specdata = (fft_data_t*)fftw_malloc(sizeof(fft_data_t)*s->fft_size); 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); 65 72 #endif 66 73 return s; 67 74 } 68 75 69 76 void del_aubio_fft(aubio_fft_t * s) { 70 /* destroy data */ 71 fftw_destroy_plan(s->pfw); 72 fftw_destroy_plan(s->pbw); 73 fftw_free(s->specdata); 74 AUBIO_FREE(s->out); 75 AUBIO_FREE(s->in ); 76 AUBIO_FREE(s); 77 /* destroy data */ 78 del_fvec(s->compspec); 79 fftw_destroy_plan(s->pfw); 80 fftw_destroy_plan(s->pbw); 81 fftw_free(s->specdata); 82 AUBIO_FREE(s->out); 83 AUBIO_FREE(s->in ); 84 AUBIO_FREE(s); 77 85 } 78 86 79 void aubio_fft_do(const aubio_fft_t * s, 80 const smpl_t * data, fft_data_t * spectrum, 81 const uint_t size) { 82 uint_t i; 83 for (i=0;i<size;i++) s->in[i] = data[i]; 84 fftw_execute(s->pfw); 85 for (i=0;i<size;i++) spectrum[i] = s->specdata[i]; 87 void 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); 86 90 } 87 91 88 void aubio_fft_rdo(const aubio_fft_t * s, 89 const fft_data_t * spectrum, 90 smpl_t * data, 91 const uint_t size) { 92 uint_t i; 93 const smpl_t renorm = 1./(smpl_t)size; 94 for (i=0;i<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; 92 void 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); 97 95 } 98 96 97 void 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); 99 104 #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 //for (i=0;i<size/2+1;i++) AUBIO_DBG("%f\n", norm[i]); 105 } 106 107 void aubio_fft_getphas(smpl_t * phas, fft_data_t * spectrum, uint_t size) { 108 uint_t i; 109 for (i=0;i<size/2+1;i++) phas[i] = ARGC(spectrum[i]); 110 //for (i=0;i<size/2+1;i++) AUBIO_DBG("%f\n", phas[i]); 111 } 112 113 void aubio_fft_getspectrum(fft_data_t * spectrum, smpl_t *norm, smpl_t * phas, uint_t size) { 114 uint_t j; 115 for (j=0; j<size/2+1; j++) { 116 spectrum[j] = CEXPC(I*phas[j]); 117 spectrum[j] *= norm[j]; 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 118 116 } 119 117 } 120 118 119 void 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 #ifdef 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]; 121 130 #else 122 123 void aubio_fft_getnorm(smpl_t * norm, fft_data_t * spectrum, uint_t size) { 124 uint_t i; 125 norm[0] = -spectrum[0]; 126 for (i=1;i<size/2+1;i++) norm[i] = SQRT(SQR(spectrum[i]) + SQR(spectrum[size-i])); 127 //for (i=0;i<size/2+1;i++) AUBIO_DBG("%f\n", norm[i]); 128 } 129 130 void aubio_fft_getphas(smpl_t * phas, fft_data_t * spectrum, uint_t size) { 131 uint_t i; 132 phas[0] = PI; 133 for (i=1;i<size/2+1;i++) phas[i] = atan2f(spectrum[size-i] , spectrum[i]); 134 //for (i=0;i<size/2+1;i++) AUBIO_DBG("%f\n", phas[i]); 135 } 136 137 void aubio_fft_getspectrum(fft_data_t * spectrum, smpl_t *norm, smpl_t * phas, uint_t size) { 138 uint_t j; 139 for (j=0; j<size/2+1; j++) { 140 spectrum[j] = norm[j]*COS(phas[j]); 141 } 142 for (j=1; j<size/2+1; j++) { 143 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 } 144 139 } 145 140 } 146 141 147 #endif 148 149 /* new interface aubio_mfft */ 150 struct _aubio_mfft_t { 151 aubio_fft_t * fft; /* fftw interface */ 152 fft_data_t ** spec; /* complex spectral data */ 153 uint_t winsize; 154 uint_t channels; 155 }; 156 157 aubio_mfft_t * new_aubio_mfft(uint_t winsize, uint_t channels){ 158 uint_t i; 159 aubio_mfft_t * fft = AUBIO_NEW(aubio_mfft_t); 160 fft->winsize = winsize; 161 fft->channels = channels; 162 fft->fft = new_aubio_fft(winsize); 163 fft->spec = AUBIO_ARRAY(fft_data_t*,channels); 164 for (i=0; i < channels; i++) 165 fft->spec[i] = AUBIO_ARRAY(fft_data_t,winsize); 166 return fft; 142 void aubio_fft_get_spectrum(fvec_t * compspec, cvec_t * spectrum) { 143 aubio_fft_get_phas(compspec, spectrum); 144 aubio_fft_get_norm(compspec, spectrum); 167 145 } 168 146 169 /* execute stft */ 170 void aubio_mfft_do (aubio_mfft_t * fft,fvec_t * in,cvec_t * fftgrain){ 171 uint_t i=0; 172 /* execute stft */ 173 for (i=0; i < fft->channels; i++) { 174 aubio_fft_do (fft->fft,in->data[i],fft->spec[i],fft->winsize); 175 /* put norm and phase into fftgrain */ 176 aubio_fft_getnorm(fftgrain->norm[i], fft->spec[i], fft->winsize); 177 aubio_fft_getphas(fftgrain->phas[i], fft->spec[i], fft->winsize); 178 } 147 void aubio_fft_get_realimag(cvec_t * spectrum, fvec_t * compspec) { 148 aubio_fft_get_imag(spectrum, compspec); 149 aubio_fft_get_real(spectrum, compspec); 179 150 } 180 151 181 /* execute inverse fourier transform */ 182 void aubio_mfft_rdo(aubio_mfft_t * fft,cvec_t * fftgrain, fvec_t * out){ 183 uint_t i=0; 184 for (i=0; i < fft->channels; i++) { 185 aubio_fft_getspectrum(fft->spec[i],fftgrain->norm[i],fftgrain->phas[i],fft->winsize); 186 aubio_fft_rdo(fft->fft,fft->spec[i],out->data[i],fft->winsize); 187 } 152 void 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.; 161 } 188 162 } 189 163 190 void del_aubio_mfft(aubio_mfft_t * fft) { 191 uint_t i; 192 for (i=0; i < fft->channels; i++) 193 AUBIO_FREE(fft->spec[i]); 194 AUBIO_FREE(fft->spec); 195 del_aubio_fft(fft->fft); 196 AUBIO_FREE(fft); 164 void 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]; 173 } 197 174 } 175 176 void 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 } 184 } 185 186 void 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.