Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • src/filterbank.c

    r45134c5 rbc4ba75  
    1919*/
    2020
    21 /* part of this mfcc implementation were inspired from LibXtract
    22    http://libxtract.sourceforge.net/
    23 */
    2421
    2522#include "aubio_priv.h"
    2623#include "sample.h"
    2724#include "filterbank.h"
     25
     26#include "stdio.h"
    2827
    2928#define USE_EQUAL_GAIN 1
     
    5352}
    5453
    55 aubio_filterbank_t * new_aubio_filterbank_mfcc(uint_t n_filters, uint_t win_s, smpl_t samplerate, smpl_t freq_min, smpl_t freq_max){
     54aubio_filterbank_t * new_aubio_filterbank_mfcc(uint_t n_filters, uint_t win_s, uint_t samplerate, smpl_t freq_min, smpl_t freq_max){
     55 
    5656  smpl_t nyquist = samplerate/2.;
    5757  uint_t style = 1;
     
    142142    for(k = next_peak + 1; k < fb->win_s; k++)
    143143      fb->filters[n]->data[0][k]=0.f;
     144
     145
    144146  }
    145147
     
    149151  free(fft_peak);
    150152
     153
    151154  return fb;
    152155
    153156}
    154157
     158/*
     159FB initialization based on Slaney's auditory toolbox
     160TODO:
     161  *solve memory leak problems while
     162  *solve quantization issues when constructing signal:
     163    *bug for win_s=512
     164    *corrections for win_s=1024 -> why even filters with smaller amplitude
     165
     166*/
     167
     168aubio_filterbank_t * new_aubio_filterbank_mfcc2(uint_t n_filters, uint_t win_s, uint_t samplerate, smpl_t freq_min, smpl_t freq_max){
     169 
     170  aubio_filterbank_t * fb = new_aubio_filterbank(n_filters, win_s);
     171 
     172 
     173  //slaney params
     174  smpl_t lowestFrequency = 133.3333;
     175  smpl_t linearSpacing = 66.66666666;
     176  smpl_t logSpacing = 1.0711703;
     177
     178  uint_t linearFilters = 13;
     179  uint_t logFilters = 27;
     180  uint_t allFilters = linearFilters + logFilters;
     181 
     182  //buffers for computing filter frequencies
     183  fvec_t * freqs=new_fvec(allFilters+2 , 1);
     184 
     185  fvec_t * lower_freqs=new_fvec( allFilters, 1);
     186  fvec_t * upper_freqs=new_fvec( allFilters, 1);
     187  fvec_t * center_freqs=new_fvec( allFilters, 1);
     188
     189 
     190  fvec_t * triangle_heights=new_fvec( allFilters, 1);
     191  //lookup table of each bin frequency in hz
     192  fvec_t * fft_freqs=new_fvec(win_s, 1);
     193
     194  uint_t filter_cnt, bin_cnt;
     195 
     196  //first step: filling all the linear filter frequencies
     197  for(filter_cnt=0; filter_cnt<linearFilters; filter_cnt++){
     198    freqs->data[0][filter_cnt]=lowestFrequency+ filter_cnt*linearSpacing;
     199  }
     200  smpl_t lastlinearCF=freqs->data[0][filter_cnt-1];
     201 
     202  //second step: filling all the log filter frequencies
     203  for(filter_cnt=0; filter_cnt<logFilters+2; filter_cnt++){
     204    freqs->data[0][filter_cnt+linearFilters]=lastlinearCF*(pow(logSpacing,filter_cnt+1));
     205  }
     206
     207  //Option 1. copying interesting values to lower_freqs, center_freqs and upper freqs arrays
     208  //TODO: would be nicer to have a reference to freqs->data, anyway we do not care in this init step
     209   
     210  for(filter_cnt=0; filter_cnt<allFilters; filter_cnt++){
     211    lower_freqs->data[0][filter_cnt]=freqs->data[0][filter_cnt];
     212    center_freqs->data[0][filter_cnt]=freqs->data[0][filter_cnt+1];
     213    upper_freqs->data[0][filter_cnt]=freqs->data[0][filter_cnt+2];
     214  }
     215
     216
     217  //computing triangle heights so that each triangle has unit area
     218  for(filter_cnt=0; filter_cnt<allFilters; filter_cnt++){
     219    triangle_heights->data[0][filter_cnt]=2./(upper_freqs->data[0][filter_cnt]-lower_freqs->data[0][filter_cnt]);
     220  }
     221 
     222 
     223  //AUBIO_DBG
     224  AUBIO_DBG("filter tables frequencies\n");
     225  for(filter_cnt=0; filter_cnt<allFilters; filter_cnt++)
     226    AUBIO_DBG("filter n. %d %f %f %f %f\n",filter_cnt, lower_freqs->data[0][filter_cnt], center_freqs->data[0][filter_cnt], upper_freqs->data[0][filter_cnt], triangle_heights->data[0][filter_cnt]);
     227
     228  //filling the fft_freqs lookup table, which assigns the frequency in hz to each bin
     229  for(bin_cnt=0; bin_cnt<win_s; bin_cnt++){
     230    //TODO: check the formula!
     231    fft_freqs->data[0][bin_cnt]= (smpl_t)samplerate* (smpl_t)bin_cnt/ (smpl_t)win_s;
     232  }
     233
     234  //building each filter table
     235  for(filter_cnt=0; filter_cnt<allFilters; filter_cnt++){
     236
     237    //TODO:check special case : lower freq =0
     238    //calculating rise increment in mag/Hz
     239    smpl_t riseInc= triangle_heights->data[0][filter_cnt]/(center_freqs->data[0][filter_cnt]-lower_freqs->data[0][filter_cnt]);
     240   
     241
     242    AUBIO_DBG("\nfilter %d",filter_cnt);
     243    //zeroing begining of filter
     244    AUBIO_DBG("\nzero begin\n");
     245    for(bin_cnt=0; bin_cnt<win_s-1; bin_cnt++){
     246      //zeroing beigining of array
     247      fb->filters[filter_cnt]->data[0][bin_cnt]=0.f;
     248      AUBIO_DBG(".");
     249      //AUBIO_DBG("%f %f %f\n", fft_freqs->data[0][bin_cnt], fft_freqs->data[0][bin_cnt+1], lower_freqs->data[0][filter_cnt]);
     250      if(fft_freqs->data[0][bin_cnt]<= lower_freqs->data[0][filter_cnt] && fft_freqs->data[0][bin_cnt+1]> lower_freqs->data[0][filter_cnt]){
     251        break;
     252      }
     253    }
     254    bin_cnt++;
     255   
     256    AUBIO_DBG("\npos slope\n");
     257    //positive slope
     258    for(; bin_cnt<win_s-1; bin_cnt++){
     259      AUBIO_DBG(".");
     260      fb->filters[filter_cnt]->data[0][bin_cnt]=(fft_freqs->data[0][bin_cnt]-lower_freqs->data[0][filter_cnt])*riseInc;
     261      //if(fft_freqs->data[0][bin_cnt]<= center_freqs->data[0][filter_cnt] && fft_freqs->data[0][bin_cnt+1]> center_freqs->data[0][filter_cnt])
     262      if(fft_freqs->data[0][bin_cnt+1]> center_freqs->data[0][filter_cnt])
     263        break;
     264    }
     265    //bin_cnt++;
     266   
     267   
     268    //negative slope
     269    AUBIO_DBG("\nneg slope\n");   
     270    for(; bin_cnt<win_s-1; bin_cnt++){
     271      //AUBIO_DBG(".");
     272     
     273      //checking whether last value is less than 0...
     274      smpl_t val=triangle_heights->data[0][filter_cnt]-(fft_freqs->data[0][bin_cnt]-center_freqs->data[0][filter_cnt])*riseInc;
     275      if(val>=0)
     276        fb->filters[filter_cnt]->data[0][bin_cnt]=val;
     277      else fb->filters[filter_cnt]->data[0][bin_cnt]=0.f;
     278     
     279      //if(fft_freqs->data[0][bin_cnt]<= upper_freqs->data[0][bin_cnt] && fft_freqs->data[0][bin_cnt+1]> upper_freqs->data[0][filter_cnt])
     280      //TODO: CHECK whether bugfix correct
     281      if(fft_freqs->data[0][bin_cnt+1]> upper_freqs->data[0][filter_cnt])
     282        break;
     283    }
     284    //bin_cnt++;
     285   
     286    AUBIO_DBG("\nzero end\n");
     287    //zeroing tail
     288    for(; bin_cnt<win_s; bin_cnt++)
     289      //AUBIO_DBG(".");
     290      fb->filters[filter_cnt]->data[0][bin_cnt]=0.f;
     291
     292  }
     293 
     294 
     295
     296  del_fvec(freqs);
     297  del_fvec(lower_freqs);
     298  del_fvec(upper_freqs);
     299  del_fvec(center_freqs);
     300
     301  del_fvec(triangle_heights);
     302  del_fvec(fft_freqs);
     303
     304  return fb;
     305
     306}
     307
     308void aubio_dump_filterbank(aubio_filterbank_t * fb){
     309
     310  FILE * mlog;
     311  mlog=fopen("filterbank.txt","w");
     312 
     313  int k,n;
     314  //dumping filter values
     315  //smpl_t area_tmp=0.f;
     316  for(n = 0; n < fb->n_filters; n++){
     317    for(k = 0; k < fb->win_s; k++){
     318      fprintf(mlog,"%f ",fb->filters[n]->data[0][k]);
     319    }
     320    fprintf(mlog,"\n");
     321  }
     322 
     323  if(mlog) fclose(mlog);
     324}
    155325
    156326void del_aubio_filterbank(aubio_filterbank_t * fb){
Note: See TracChangeset for help on using the changeset viewer.