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@1: 
jamie@1: /* xtract_vector.c: defines functions that extract a feature as a single value from an input vector */
jamie@1: 
jamie@1: #include "xtract/libxtract.h"
jamie@1: #include <math.h>
jamie@30: 
jamie@30: #ifdef XTRACT_FFT
jamie@30: 
jamie@1: #include <fftw3.h>
jamie@1: 
jamie@1: int xtract_magnitude_spectrum(float *data, int N, void *argv, float *result){
jamie@1: 
jamie@1:     float *temp;
jamie@1:     int n , M = N >> 1;
jamie@1:     fftwf_plan plan;
jamie@1: 
jamie@1:     temp = (float *)fftwf_malloc(N * sizeof(float));
jamie@1: 
jamie@1:     plan = fftwf_plan_r2r_1d(N, data, temp, FFTW_R2HC, FFTW_ESTIMATE);
jamie@1:     
jamie@1:     fftwf_execute(plan);
jamie@1:     
jamie@1:     for(n = 1; n < M; n++){
jamie@1:         result[n] = sqrt(SQ(temp[n]) + SQ(temp[N - n])) / N;
jamie@1:         result[N-n] = 0.0f;
jamie@1:     }
jamie@1:     
jamie@1:     result[0] = fabs(temp[0]) / N;
jamie@1:     result[M] = fabs(temp[M]) / N;
jamie@1:     
jamie@1:     fftwf_destroy_plan(plan);
jamie@1:     fftwf_free(temp);
jamie@1:     
jamie@1: }
jamie@1: 
jamie@1: int xtract_autocorrelation_fft(float *data, int N, void *argv, float *result){
jamie@1:     
jamie@1:     float *temp;
jamie@1:     int n;
jamie@1:     fftwf_plan plan;
jamie@1: 
jamie@1:     temp = (float *)fftwf_malloc(N * sizeof(float));
jamie@1:     plan = fftwf_plan_r2r_1d(N, data, temp, FFTW_HC2R, FFTW_ESTIMATE);
jamie@1: 
jamie@1:     fftwf_execute(plan);
jamie@1:     
jamie@1:     for(n = 0; n < N - 1; n++)
jamie@1:         result[n] = temp[n+1];
jamie@1:     
jamie@1:     fftwf_destroy_plan(plan);
jamie@1:     fftwf_free(temp);
jamie@1: }
jamie@1: 
jamie@30: int xtract_mfcc(float *data, int N, void *argv, float *result){
jamie@30: 
jamie@30:     xtract_mel_filter *f;
jamie@30:     int n, filter;
jamie@30:     fftwf_plan plan;
jamie@30: 
jamie@30:     f = (xtract_mel_filter *)argv;
jamie@30: 
jamie@30:     for(filter = 0; filter < f->n_filters; filter++){
jamie@30:         for(n = 0; n < N; n++){
jamie@30:             result[filter] += data[n] * f->filters[filter][n];
jamie@30:         }
jamie@30:         if(result[filter] < LOG_LIMIT) result[filter] = LOG_LIMIT;
jamie@30:         result[filter] = log(result[filter]);
jamie@30:     }
jamie@30: 
jamie@30:     for(n = filter + 1; n < N; n++) result[n] = 0; 
jamie@30:     
jamie@30:     xtract_dct(result, f->n_filters, NULL, result);
jamie@30: 
jamie@30: }
jamie@30: 
jamie@30: int xtract_dct(float *data, int N, void *argv, float *result){
jamie@30:     
jamie@30:     fftwf_plan plan;
jamie@30: 
jamie@30:     plan = 
jamie@30:         fftwf_plan_r2r_1d(N, data, result, FFTW_REDFT00, FFTW_ESTIMATE);
jamie@30:     
jamie@30:     fftwf_execute(plan);
jamie@30:     fftwf_destroy_plan(plan);
jamie@30: }
jamie@30: 
jamie@30: #else
jamie@30: 
jamie@30: int xtract_magnitude_spectrum(float *data, int N, void *argv, float *result){
jamie@30: 
jamie@30:     NOT_IMPLEMENTED;
jamie@30: 
jamie@30: }
jamie@30: 
jamie@30: int xtract_autocorrelation_fft(float *data, int N, void *argv, float *result){
jamie@30: 
jamie@30:     NOT_IMPLEMENTED;
jamie@30: 
jamie@30: }
jamie@30: 
jamie@30: int xtract_mfcc(float *data, int N, void *argv, float *result){
jamie@30: 
jamie@30:     NOT_IMPLEMENTED;
jamie@30: 
jamie@30: }
jamie@30: 
jamie@30: int xtract_dct(float *data, int N, void *argv, float *result){
jamie@30: 
jamie@30:     NOT_IMPLEMENTED;
jamie@30: 
jamie@30: }
jamie@30: 
jamie@30: #endif
jamie@30: 
jamie@30: int xtract_autocorrelation(float *data, int N, void *argv, float *result){
jamie@30: 
jamie@30:     /* Naive time domain implementation  */
jamie@30:     
jamie@30:     int n = N, i;
jamie@30:     
jamie@30:     float corr;
jamie@30: 
jamie@30:     while(n--){
jamie@30:        corr = 0;
jamie@30:         for(i = 0; i < N - n; i++){
jamie@30:             corr += data[i] * data[i + n];
jamie@30:         }
jamie@30:         result[n] = corr / N;
jamie@30:     }
jamie@30: }
jamie@30: 
jamie@1: int xtract_amdf(float *data, int N, void *argv, float *result){
jamie@1: 
jamie@1:     int n = N, i;
jamie@1:     
jamie@6:     float md, temp;
jamie@1: 
jamie@1:     while(n--){
jamie@1:        md = 0;
jamie@1:         for(i = 0; i < N - n; i++){
jamie@6:             temp = data[i] - data[i + n];
jamie@6: 			temp = (temp < 0 ? -temp : temp);
jamie@6: 			md += temp;
jamie@1:         }
jamie@1:         result[n] = md / N;
jamie@1:     }
jamie@1: }
jamie@1: 
jamie@1: int xtract_asdf(float *data, int N, void *argv, float *result){
jamie@1:     
jamie@1:     int n = N, i;
jamie@1:     
jamie@1:     float sd;
jamie@1: 
jamie@1:     while(n--){
jamie@1:        sd = 0;
jamie@1:         for(i = 0; i < N - n; i++){
jamie@6:             /*sd = 1;*/
jamie@1:             sd += SQ(data[i] - data[i + n]);
jamie@1:         }
jamie@1:         result[n] = sd / N;
jamie@1:     }
jamie@1: }
jamie@1: 
jamie@1: int xtract_bark_coefficients(float *data, int N, void *argv, float *result){
jamie@1: 
jamie@1:     int *limits, band, n;
jamie@1: 
jamie@1:     limits = (int *)argv;
jamie@1:     
jamie@1:     for(band = 0; band < BARK_BANDS; band++){
jamie@1:         for(n = limits[band]; n < limits[band + 1]; n++)
jamie@1:             result[band] += data[n];
jamie@1:     }
jamie@1: }
jamie@1: 
jamie@1: int xtract_peaks(float *data, int N, void *argv, float *result){
jamie@1: 
jamie@1:     float thresh, max, y, y2, y3, p, width, sr; 
jamie@1:     int n = N, M, return_code;
jamie@1:     
jamie@1:     if(argv != NULL){
jamie@1:         thresh = ((float *)argv)[0];
jamie@1:         sr = ((float *)argv)[1];
jamie@1:         return_code = BAD_ARGV;
jamie@1:     }
jamie@1:     else{
jamie@1:         thresh = 0;
jamie@1:         sr = 44100;
jamie@1:     }
jamie@1: 
jamie@1:     M = N >> 1;
jamie@1:     width = sr / N;
jamie@1:     
jamie@1:     y = y2 = y3 = p = max = 0;
jamie@1: 
jamie@1:     if(thresh < 0 || thresh > 100){
jamie@1:         thresh = 0;
jamie@1:         return_code = BAD_ARGV;
jamie@1:     }
jamie@1: 
jamie@1:     if(!sr){
jamie@1:         sr = 44100;
jamie@1:         return_code = BAD_ARGV;
jamie@1:     }
jamie@1: 
jamie@1:     while(n--){
jamie@1:         max = MAX(max, data[n]);
jamie@1:         /* ensure we never take log10(0) */
jamie@1:         /*data[n] = (data[n] < LOG_LIMIT ? LOG_LIMIT : data[n]);*/
jamie@1:         if ((data[n] * 100000) <= 1)
jamie@1:         /* We get a more stable peak this way */
jamie@1: 		    data[n] = 1;
jamie@1:         
jamie@1:     }
jamie@1:     
jamie@1:     thresh *= .01 * max;
jamie@1: 
jamie@1:     result[0] = 0;
jamie@1:     result[M] = 0;
jamie@1: 
jamie@1:     for(n = 1; n < M; n++){
jamie@1:         if(data[n] >= thresh){
jamie@1:             if(data[n] > data[n - 1] && data[n] > data[n + 1]){
jamie@1:                 result[n] = width * (n + (p = .5 * (y = 20 * log10(data[n-1]) - (y3 = 20 * log10(data[n+1]))) / (20 * log10(data[n - 1]) - 2 * (y2 = 20 * log10(data[n])) + 20 * log10(data[n + 1]))));
jamie@1:                 result[M + n] = y2 - .25 * (y - y3) * p;
jamie@1:             }
jamie@1:             else{
jamie@1:                 result[n] = 0;
jamie@1:                 result[M + n] = 0;
jamie@1:             }
jamie@1:         }
jamie@1:         else{
jamie@1:             result[n] = 0;
jamie@1:             result[M + n] = 0;
jamie@1:         }
jamie@1:     }	  
jamie@1:    
jamie@1:     return (return_code ? return_code : SUCCESS);
jamie@1: }
jamie@1: