jamie@1: /* libxtract feature extraction library jamie@1: * jamie@1: * Copyright (C) 2006 Jamie Bullock jamie@1: * jamie@1: * This program is free software; you can redistribute it and/or modify jamie@1: * it under the terms of the GNU General Public License as published by jamie@1: * the Free Software Foundation; either version 2 of the License, or jamie@1: * (at your option) any later version. jamie@1: * jamie@1: * This program is distributed in the hope that it will be useful, jamie@1: * but WITHOUT ANY WARRANTY; without even the implied warranty of jamie@1: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the jamie@1: * GNU General Public License for more details. jamie@1: * jamie@1: * You should have received a copy of the GNU General Public License jamie@1: * along with this program; if not, write to the Free Software jamie@1: * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, jamie@1: * USA. jamie@1: */ jamie@1: jamie@107: /* init.c: defines initialisation and free functions. Also contains library constructor routine. */ jamie@1: jamie@98: #ifdef HAVE_CONFIG_H jamie@98: # include jamie@98: #endif jamie@98: jamie@1: #include jamie@26: #include jamie@1: jamie@98: #include "xtract/libxtract.h" jamie@107: #include "xtract_window_private.h" jamie@102: #define DEFINE_GLOBALS jamie@98: #include "xtract_globals_private.h" jamie@98: jamie@98: #ifdef XTRACT_FFT jamie@98: #include jamie@98: jamie@98: #ifndef XTRACT_FFT_OPTIMISATION_LEVEL jamie@98: /* This should never happen */ jamie@98: #define XTRACT_FFT_OPTIMISATION_LEVEL 1 jamie@98: #endif jamie@98: jamie@43: int xtract_init_mfcc(int N, float nyquist, int style, float freq_min, float freq_max, int freq_bands, float **fft_tables){ jamie@1: jamie@39: int n, i, k, *fft_peak, M, next_peak; jamie@1: float norm, mel_freq_max, mel_freq_min, norm_fact, height, inc, val, jamie@107: freq_bw_mel, *mel_peak, *height_norm, *lin_peak; jamie@1: jamie@1: mel_peak = height_norm = lin_peak = NULL; jamie@1: fft_peak = NULL; jamie@1: norm = 1; jamie@1: jamie@1: mel_freq_max = 1127 * log(1 + freq_max / 700); jamie@1: mel_freq_min = 1127 * log(1 + freq_min / 700); jamie@1: freq_bw_mel = (mel_freq_max - mel_freq_min) / freq_bands; jamie@1: jamie@1: mel_peak = (float *)malloc((freq_bands + 2) * sizeof(float)); jamie@1: /* +2 for zeros at start and end */ jamie@1: lin_peak = (float *)malloc((freq_bands + 2) * sizeof(float)); jamie@1: fft_peak = (int *)malloc((freq_bands + 2) * sizeof(int)); jamie@1: height_norm = (float *)malloc(freq_bands * sizeof(float)); jamie@1: jamie@1: if(mel_peak == NULL || height_norm == NULL || jamie@107: lin_peak == NULL || fft_peak == NULL) jamie@107: return XTRACT_MALLOC_FAILED; jamie@107: jamie@1: M = N >> 1; jamie@1: jamie@1: mel_peak[0] = mel_freq_min; danstowell@95: lin_peak[0] = freq_min; // === 700 * (exp(mel_peak[0] / 1127) - 1); jamie@1: fft_peak[0] = lin_peak[0] / nyquist * M; jamie@1: jamie@1: danstowell@100: for (n = 1; n < freq_bands + 2; n++){ danstowell@100: //roll out peak locations - mel, linear and linear on fft window scale jamie@1: mel_peak[n] = mel_peak[n - 1] + freq_bw_mel; jamie@1: lin_peak[n] = 700 * (exp(mel_peak[n] / 1127) -1); jamie@1: fft_peak[n] = lin_peak[n] / nyquist * M; jamie@1: } jamie@1: jamie@1: for (n = 0; n < freq_bands; n++){ danstowell@100: //roll out normalised gain of each peak jamie@56: if (style == XTRACT_EQUAL_GAIN){ jamie@1: height = 1; jamie@1: norm_fact = norm; jamie@1: } jamie@1: else{ jamie@1: height = 2 / (lin_peak[n + 2] - lin_peak[n]); jamie@1: norm_fact = norm / (2 / (lin_peak[2] - lin_peak[0])); jamie@1: } jamie@1: height_norm[n] = height * norm_fact; jamie@1: } jamie@1: jamie@1: i = 0; jamie@107: jamie@1: for(n = 0; n < freq_bands; n++){ jamie@107: jamie@107: // calculate the rise increment danstowell@95: if(n==0) danstowell@95: inc = height_norm[n] / fft_peak[n]; danstowell@95: else jamie@1: inc = height_norm[n] / (fft_peak[n] - fft_peak[n - 1]); jamie@1: val = 0; jamie@107: jamie@107: // zero the start of the array jamie@107: for(k = 0; k < i; k++) jamie@107: fft_tables[n][k] = 0.f; jamie@107: jamie@107: // fill in the rise jamie@1: for(; i <= fft_peak[n]; i++){ jamie@1: fft_tables[n][i] = val; jamie@1: val += inc; jamie@1: } jamie@107: danstowell@95: // calculate the fall increment jamie@1: inc = height_norm[n] / (fft_peak[n + 1] - fft_peak[n]); jamie@107: jamie@1: val = 0; jamie@107: next_peak = fft_peak[n + 1]; jamie@107: jamie@107: // reverse fill the 'fall' jamie@39: for(i = next_peak; i > fft_peak[n]; i--){ jamie@1: fft_tables[n][i] = val; jamie@1: val += inc; jamie@1: } jamie@39: jamie@107: // zero the rest of the array jamie@107: for(k = next_peak + 1; k < N; k++) jamie@107: fft_tables[n][k] = 0.f; jamie@1: } jamie@1: jamie@98: jamie@98: /* Initialise the fft_plan for the DCT */ jamie@98: xtract_init_fft(freq_bands, XTRACT_MFCC); jamie@98: jamie@1: free(mel_peak); jamie@1: free(lin_peak); jamie@1: free(height_norm); jamie@1: free(fft_peak); jamie@1: jamie@56: return XTRACT_SUCCESS; jamie@1: jamie@1: } jamie@1: jamie@98: int xtract_init_fft(int N, int feature_name){ jamie@98: jamie@98: float *input, *output; jamie@98: int optimisation; jamie@107: jamie@98: input = output = NULL; jamie@107: jamie@113: //fprintf(stderr, "Optimisation level: %d\n", XTRACT_FFT_OPTIMISATION_LEVEL); jamie@98: jamie@98: if(XTRACT_FFT_OPTIMISATION_LEVEL == 0) jamie@98: optimisation = FFTW_ESTIMATE; jamie@98: else if(XTRACT_FFT_OPTIMISATION_LEVEL == 1) jamie@98: optimisation = FFTW_MEASURE; jamie@98: else if(XTRACT_FFT_OPTIMISATION_LEVEL == 2) jamie@98: optimisation = FFTW_PATIENT; jamie@98: else jamie@98: optimisation = FFTW_MEASURE; /* The default */ jamie@98: jamie@98: if(feature_name == XTRACT_AUTOCORRELATION_FFT) jamie@98: N <<= 1; jamie@98: jamie@98: input = malloc(N * sizeof(float)); jamie@98: output = malloc(N * sizeof(float)); jamie@98: jamie@98: switch(feature_name){ jamie@98: case XTRACT_SPECTRUM: jamie@102: if(fft_plans.spectrum_plan != NULL) jamie@102: fftwf_destroy_plan(fft_plans.spectrum_plan); jamie@102: fft_plans.spectrum_plan = jamie@98: fftwf_plan_r2r_1d(N, input, output, FFTW_R2HC, optimisation); jamie@98: break; jamie@98: case XTRACT_AUTOCORRELATION_FFT: jamie@102: if(fft_plans.autocorrelation_fft_plan_1 != NULL) jamie@102: fftwf_destroy_plan(fft_plans.autocorrelation_fft_plan_1); jamie@102: if(fft_plans.autocorrelation_fft_plan_2 != NULL) jamie@102: fftwf_destroy_plan(fft_plans.autocorrelation_fft_plan_2); jamie@102: fft_plans.autocorrelation_fft_plan_1 = jamie@98: fftwf_plan_r2r_1d(N, input, output, FFTW_R2HC, optimisation); jamie@102: fft_plans.autocorrelation_fft_plan_2 = jamie@98: fftwf_plan_r2r_1d(N, input, output, FFTW_HC2R, optimisation); jamie@98: break; jamie@98: case XTRACT_DCT: jamie@102: if(fft_plans.dct_plan != NULL) jamie@102: fftwf_destroy_plan(fft_plans.dct_plan); jamie@102: fft_plans.dct_plan = jamie@118: fftwf_plan_r2r_1d(N, input, output, FFTW_REDFT10, optimisation); jamie@98: case XTRACT_MFCC: jamie@102: if(fft_plans.dct_plan != NULL) jamie@102: fftwf_destroy_plan(fft_plans.dct_plan); jamie@102: fft_plans.dct_plan = jamie@98: fftwf_plan_r2r_1d(N, output, output, FFTW_REDFT00, optimisation); jamie@98: break; jamie@98: } jamie@98: jamie@98: free(input); jamie@98: free(output); jamie@98: jamie@98: return XTRACT_SUCCESS; jamie@98: jamie@98: } jamie@98: jamie@110: void xtract_free_fft(void){ jamie@110: if(fft_plans.spectrum_plan != NULL) jamie@110: fftwf_destroy_plan(fft_plans.spectrum_plan); jamie@110: if(fft_plans.autocorrelation_fft_plan_1 != NULL) jamie@110: fftwf_destroy_plan(fft_plans.autocorrelation_fft_plan_1); jamie@110: if(fft_plans.autocorrelation_fft_plan_2 != NULL) jamie@110: fftwf_destroy_plan(fft_plans.autocorrelation_fft_plan_2); jamie@110: if(fft_plans.dct_plan != NULL) jamie@110: fftwf_destroy_plan(fft_plans.dct_plan); jamie@110: // fftwf_cleanup(); jamie@110: } jamie@110: jamie@98: #endif jamie@98: jamie@59: int xtract_init_bark(int N, float sr, int *band_limits){ jamie@1: jamie@38: float edges[] = {0, 100, 200, 300, 400, 510, 630, 770, 920, 1080, 1270, 1480, 1720, 2000, 2320, 2700, 3150, 3700, 4400, 5300, 6400, 7700, 9500, 12000, 15500, 20500, 27000}; /* Takes us up to sr = 54kHz (CCRMA: JOS)*/ jamie@1: jamie@59: int bands = XTRACT_BARK_BANDS; jamie@107: jamie@1: while(bands--) jamie@59: band_limits[bands] = edges[bands] / sr * N; jamie@107: /*FIX shohuld use rounding, but couldn't get it to work */ jamie@38: jamie@56: return XTRACT_SUCCESS; jamie@1: } jamie@98: jamie@107: float *xtract_init_window(const int N, const int type){ jamie@107: jamie@107: float *window; jamie@107: jamie@107: window = malloc(N * sizeof(float)); jamie@107: jamie@107: switch (type) { jamie@107: case XTRACT_GAUSS: jamie@107: gauss(window, N, 0.4); jamie@107: break; jamie@107: case XTRACT_HAMMING: jamie@107: hamming(window, N); jamie@107: break; jamie@107: case XTRACT_HANN: jamie@107: hann(window, N); jamie@107: break; jamie@107: case XTRACT_BARTLETT: jamie@107: bartlett(window, N); jamie@107: break; jamie@107: case XTRACT_TRIANGULAR: jamie@107: triangular(window, N); jamie@107: break; jamie@107: case XTRACT_BARTLETT_HANN: jamie@107: bartlett_hann(window, N); jamie@107: break; jamie@107: case XTRACT_BLACKMAN: jamie@107: blackman(window, N); jamie@107: break; jamie@107: case XTRACT_KAISER: jamie@107: kaiser(window, N, 3 * PI); jamie@107: break; jamie@107: case XTRACT_BLACKMAN_HARRIS: jamie@107: blackman_harris(window, N); jamie@107: break; jamie@107: default: jamie@107: hann(window, N); jamie@107: break; jamie@107: } jamie@107: jamie@107: return window; jamie@107: jamie@107: } jamie@107: jamie@107: void xtract_free_window(float *window){ jamie@107: jamie@107: free(window); jamie@107: jamie@107: } jamie@107: jamie@102: #ifdef __GNUC__ jamie@102: __attribute__((constructor)) void init() jamie@102: #else jamie@102: void _init()ยท jamie@102: #endif jamie@102: { jamie@116: #ifdef XTRACT_FFT jamie@102: fft_plans.spectrum_plan = NULL; jamie@102: fft_plans.autocorrelation_fft_plan_1 = NULL; jamie@102: fft_plans.autocorrelation_fft_plan_2 = NULL; jamie@102: fft_plans.dct_plan = NULL; jamie@116: #endif jamie@102: }