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@56: #include "xtract_macros_private.h"
jamie@1: #include <math.h>
jamie@43: #include <string.h>
jamie@43: #include <stdlib.h>
jamie@30: 
jamie@30: #ifdef XTRACT_FFT
jamie@30: 
jamie@1: #include <fftw3.h>
jamie@1: 
jamie@54: int xtract_spectrum(const float *data, const int N, const void *argv, float *result){
jamie@1: 
jamie@56:     float *input, *rfft, q, temp;
jamie@43:     size_t bytes;
jamie@54:     int n , NxN, M, vector;
jamie@1:     fftwf_plan plan;
jamie@1: 
jamie@54:     M = N >> 1;
jamie@56:     NxN = XTRACT_SQ(N);
jamie@54: 
jamie@54:     rfft = (float *)fftwf_malloc(N * sizeof(float));
jamie@43:     input = (float *)malloc(bytes = N * sizeof(float));
jamie@43:     input = memcpy(input, data, bytes);
jamie@1: 
jamie@56:     q = *(float *)argv;
jamie@54:     vector = (int)*((float *)argv+1);
jamie@46: 
jamie@56:     XTRACT_CHECK_q;
jamie@46: 
jamie@54:     plan = fftwf_plan_r2r_1d(N, input, rfft, FFTW_R2HC, FFTW_ESTIMATE);
jamie@1:     
jamie@1:     fftwf_execute(plan);
jamie@54: 
jamie@54:     switch(vector){
jamie@67: 
jamie@67: /*	case XTRACT_MAGNITUDE_SPECTRUM:
jamie@67: 	    for(n = 1; n < M; n++){
jamie@67: 		result[n] = sqrt(XTRACT_SQ(rfft[n]) + 
jamie@67: 			XTRACT_SQ(rfft[N - n - 1])) / N; 
jamie@56: 		result[M + n] = n * q;
jamie@54: 	    }
jamie@54: 	    break;
jamie@67: */
jamie@56: 	case XTRACT_LOG_MAGNITUDE_SPECTRUM:
jamie@67: 	    for(n = 1; n < M; n++){
jamie@67: 		if ((temp = XTRACT_SQ(rfft[n]) + 
jamie@67: 			    XTRACT_SQ(rfft[N - n - 1])) > XTRACT_LOG_LIMIT)
jamie@54: 		    temp = log(sqrt(temp) / N);
jamie@54: 		else
jamie@56: 		    temp = XTRACT_LOG_LIMIT_DB;
jamie@54: 		/*Normalise*/
jamie@67: 		result[n] = 
jamie@67: 		    (temp + XTRACT_DB_SCALE_OFFSET) / XTRACT_DB_SCALE_OFFSET; 
jamie@56: 		result[M + n] = n * q;
jamie@54: 	    }
jamie@54: 	    break;
jamie@67: 
jamie@56: 	case XTRACT_POWER_SPECTRUM:
jamie@67: 	    for(n = 1; n < M; n++){
jamie@67: 		result[n] = (XTRACT_SQ(rfft[n]) + XTRACT_SQ(rfft[N - n - 1])) 
jamie@67: 		    / NxN;
jamie@56: 		result[M + n] = n * q;
jamie@54: 	    }
jamie@54: 	    break;
jamie@67: 
jamie@56: 	case XTRACT_LOG_POWER_SPECTRUM:
jamie@67: 	    for(n = 1; n < M; n++){
jamie@67: 		if ((temp = XTRACT_SQ(rfft[n]) + XTRACT_SQ(rfft[N - n - 1])) > 
jamie@67: 			XTRACT_LOG_LIMIT)
jamie@54: 		    temp = log(temp / NxN);
jamie@54: 		else
jamie@56: 		    temp = XTRACT_LOG_LIMIT_DB; 		
jamie@67: 		result[n] = (temp + XTRACT_DB_SCALE_OFFSET) / 
jamie@67: 		    XTRACT_DB_SCALE_OFFSET; 
jamie@56: 		result[M + n] = n * q;
jamie@54: 	    }
jamie@54: 	    break;
jamie@67: 
jamie@54: 	default:
jamie@54: 	    /* MAGNITUDE_SPECTRUM */
jamie@67: 	    for(n = 1; n < M; n++){
jamie@67: 		result[n] = sqrt(XTRACT_SQ(rfft[n]) + 
jamie@67: 			XTRACT_SQ(rfft[N - n - 1])) / N; 
jamie@56: 		result[M + n] = n * q;
jamie@54: 	    }
jamie@54: 	    break;
jamie@1:     }
jamie@1:     
jamie@67:     /* Set the DC component to 0 */
jamie@67:     result[0] = result[M] = 0.f;
jamie@67:     /* Set the Nyquist */
jamie@56:     result[N] = q * M;
jamie@1:     
jamie@1:     fftwf_destroy_plan(plan);
jamie@54:     fftwf_free(rfft);
jamie@43:     free(input);
jamie@1:     
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@43: int xtract_autocorrelation_fft(const float *data, const int N, const void *argv, float *result){
jamie@1:     
jamie@43:     float *temp, *input;
jamie@43:     size_t bytes;
jamie@1:     int n;
jamie@1:     fftwf_plan plan;
jamie@1: 
jamie@1:     temp = (float *)fftwf_malloc(N * sizeof(float));
jamie@43:     input = (float *)malloc(bytes = N * sizeof(float));
jamie@43:     input = memcpy(input, data, bytes);
jamie@43: 
jamie@43:     plan = fftwf_plan_r2r_1d(N, input, 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@43:     free(input);
jamie@38: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@43: int xtract_mfcc(const float *data, const int N, const void *argv, float *result){
jamie@30: 
jamie@30:     xtract_mel_filter *f;
jamie@43:     float *input;
jamie@43:     size_t bytes;
jamie@30:     int n, filter;
jamie@30: 
jamie@30:     f = (xtract_mel_filter *)argv;
jamie@39:     
jamie@43:     input = (float *)malloc(bytes = N * sizeof(float));
jamie@43:     input = memcpy(input, data, bytes);
jamie@43:     
jamie@30:     for(filter = 0; filter < f->n_filters; filter++){
danstowell@68:         result[filter] = 0.f;
jamie@30:         for(n = 0; n < N; n++){
jamie@43:             result[filter] += input[n] * f->filters[filter][n];
jamie@30:         }
jamie@56:         if(result[filter] < XTRACT_LOG_LIMIT) result[filter] = XTRACT_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@43:     free(input);
jamie@43:     
jamie@56:     return XTRACT_SUCCESS;
jamie@30: }
jamie@30: 
jamie@43: int xtract_dct(const float *data, const int N, const void *argv, float *result){
jamie@30:     
jamie@30:     fftwf_plan plan;
jamie@43:     float *input;
jamie@43:     size_t bytes;
jamie@30: 
jamie@43:     input = (float *)malloc(bytes = N * sizeof(float));
jamie@43:     input = memcpy(input, data, bytes);
jamie@43:     
jamie@30:     plan = 
jamie@43:         fftwf_plan_r2r_1d(N, input, result, FFTW_REDFT00, FFTW_ESTIMATE);
jamie@30:     
jamie@30:     fftwf_execute(plan);
jamie@30:     fftwf_destroy_plan(plan);
jamie@43:     free(input);
jamie@38: 
jamie@56:     return XTRACT_SUCCESS;
jamie@30: }
jamie@30: 
jamie@30: #else
jamie@30: 
jamie@67: int xtract_spectrum(const float *data, const int N, const void *argv, float *result){
jamie@30: 
danstowell@66:     XTRACT_NEEDS_FFTW;
danstowell@66:     return XTRACT_NO_RESULT;
jamie@30: 
jamie@30: }
jamie@30: 
jamie@43: int xtract_autocorrelation_fft(const float *data, const int N, const void *argv, float *result){
jamie@30: 
danstowell@66:     XTRACT_NEEDS_FFTW;
danstowell@66:     return XTRACT_NO_RESULT;
jamie@30: 
jamie@30: }
jamie@30: 
jamie@43: int xtract_mfcc(const float *data, const int N, const void *argv, float *result){
jamie@30: 
danstowell@66:     XTRACT_NEEDS_FFTW;
danstowell@66:     return XTRACT_NO_RESULT;
jamie@30: 
jamie@30: }
jamie@30: 
jamie@43: int xtract_dct(const float *data, const int N, const void *argv, float *result){
jamie@30: 
danstowell@66:     XTRACT_NEEDS_FFTW;
danstowell@66:     return XTRACT_NO_RESULT;
jamie@30: 
jamie@30: }
jamie@30: 
jamie@30: #endif
jamie@30: 
jamie@43: int xtract_autocorrelation(const float *data, const int N, const 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@38: 
jamie@56:     return XTRACT_SUCCESS;
jamie@30: }
jamie@30: 
jamie@43: int xtract_amdf(const float *data, const int N, const 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@38: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@43: int xtract_asdf(const float *data, const int N, const 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@56:             sd += XTRACT_SQ(data[i] - data[i + n]);
jamie@1:         }
jamie@1:         result[n] = sd / N;
jamie@1:     }
jamie@38: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@43: int xtract_bark_coefficients(const float *data, const int N, const void *argv, float *result){
jamie@1: 
jamie@1:     int *limits, band, n;
jamie@1: 
jamie@1:     limits = (int *)argv;
jamie@1:     
jamie@59:     for(band = 0; band < XTRACT_BARK_BANDS - 1; band++){
jamie@1:         for(n = limits[band]; n < limits[band + 1]; n++)
jamie@1:             result[band] += data[n];
jamie@1:     }
jamie@38: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@52: int xtract_peak_spectrum(const float *data, const int N, const void *argv, float *result){
jamie@1: 
jamie@56:     float threshold, max, y, y2, y3, p, q, *input = NULL;
jamie@43:     size_t bytes;
jamie@59:     int n = N, rv = XTRACT_SUCCESS;
jamie@49: 
jamie@56:     threshold = max = y = y2 = y3 = p = q = 0.f;
jamie@1:     
jamie@1:     if(argv != NULL){
jamie@56:         q = ((float *)argv)[0];
jamie@55:         threshold = ((float *)argv)[1];
jamie@1:     }
jamie@49:     else
jamie@56:         rv = XTRACT_BAD_ARGV;
jamie@49: 
jamie@55:     if(threshold < 0 || threshold > 100){
jamie@55:         threshold = 0;
jamie@56:         rv = XTRACT_BAD_ARGV;
jamie@1:     }
jamie@1: 
jamie@56:     XTRACT_CHECK_q;
jamie@49: 
jamie@43:     input = (float *)malloc(bytes = N * sizeof(float));
jamie@43: 
jamie@43:     if(input != NULL)
jamie@43: 	input = memcpy(input, data, bytes);
jamie@43:     else
jamie@56: 	return XTRACT_MALLOC_FAILED;
jamie@43: 
jamie@45:     while(n--)
jamie@56:         max = XTRACT_MAX(max, input[n]);
jamie@1:     
jamie@55:     threshold *= .01 * max;
jamie@1: 
jamie@1:     result[0] = 0;
jamie@59:     result[N] = 0;
jamie@1: 
jamie@59:     for(n = 1; n < N; n++){
jamie@55:         if(input[n] >= threshold){
jamie@43:             if(input[n] > input[n - 1] && input[n] > input[n + 1]){
jamie@59:                 result[N + n] = q * (n + (p = .5 * (y = input[n-1] - 
jamie@52: 				(y3 = input[n+1])) / (input[n - 1] - 2 * 
jamie@52: 				    (y2 = input[n]) + input[n + 1])));
jamie@52:                 result[n] = y2 - .25 * (y - y3) * p;
jamie@1:             }
jamie@1:             else{
jamie@1:                 result[n] = 0;
jamie@59:                 result[N + n] = 0;
jamie@1:             }
jamie@1:         }
jamie@1:         else{
jamie@1:             result[n] = 0;
jamie@59:             result[N + n] = 0;
jamie@1:         }
jamie@1:     }	  
jamie@1:    
jamie@43:     free(input);
jamie@56:     return (rv ? rv : XTRACT_SUCCESS);
jamie@1: }
jamie@41: 	    
jamie@52: int xtract_harmonic_spectrum(const float *data, const int N, const void *argv, float *result){
jamie@38:     
jamie@38:     int n = (N >> 1), M = n; 
jamie@38: 
jamie@43:     const float *freqs, *amps;
jamie@55:     float f0, threshold, ratio, nearest, distance;
jamie@38: 
jamie@52:     amps = data;
jamie@52:     freqs = data + n;
jamie@38:     f0 = *((float *)argv);
jamie@55:     threshold = *((float *)argv+1);
jamie@38: 
jamie@38:     ratio = nearest = distance = 0.f;
jamie@38: 
jamie@38:     while(n--){
jamie@38: 	if(freqs[n]){
jamie@38: 	    ratio = freqs[n] / f0;
jamie@38: 	    nearest = round(ratio);
jamie@38: 	    distance = fabs(nearest - ratio);
jamie@55: 	    if(distance > threshold)
jamie@38: 		result[n] = result[M + n] = 0.f;
jamie@42: 	    else {
jamie@52: 		result[n] = amps[n];
jamie@52: 		result[M + n] = freqs[n];
jamie@42: 	    }
jamie@38: 	}
jamie@38: 	else
jamie@38: 	    result[n] = result[M + n] = 0.f;
jamie@38:     }
jamie@56:     return XTRACT_SUCCESS;
jamie@38: }
jamie@38: