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_scalar.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@5: #include <stdlib.h>
jamie@1: 
jamie@1: int xtract_mean(float *data, int N, void *argv, float *result){
jamie@25: 
jamie@1:     int n = N;
jamie@1: 
jamie@1:     while(n--)
jamie@25: 	*result += *data++;
jamie@25: 
jamie@1:     *result /= N;
jamie@1: }
jamie@1: 
jamie@1: int xtract_variance(float *data, int N, void *argv, float *result){
jamie@25: 
jamie@1:     int n = N;
jamie@1: 
jamie@1:     while(n--)
jamie@25: 	*result += *data++ - *(float *)argv;
jamie@25: 
jamie@1:     *result = SQ(*result) / (N - 1);
jamie@1: }
jamie@1: 
jamie@1: int xtract_standard_deviation(float *data, int N, void *argv, float *result){
jamie@25: 
jamie@1:     *result = sqrt(*(float *)argv);
jamie@25: 
jamie@1: }
jamie@1: 
jamie@1: int xtract_average_deviation(float *data, int N, void *argv, float *result){
jamie@25: 
jamie@1:     int n = N;
jamie@1: 
jamie@1:     while(n--)
jamie@25: 	*result += fabs(*data++ - *(float *)argv);
jamie@25: 
jamie@1:     *result /= N;
jamie@1: 
jamie@1: }
jamie@1: 
jamie@1: int xtract_skewness(float *data, int N, void *argv,  float *result){
jamie@25: 
jamie@1:     int n = N;
jamie@1: 
jamie@1:     while(n--)
jamie@25: 	*result += (*data++ - ((float *)argv)[0]) / ((float *)argv)[1];
jamie@25: 
jamie@1:     *result = pow(*result, 3) / N;
jamie@1: 
jamie@1: }
jamie@1: 
jamie@1: int xtract_kurtosis(float *data, int N, void *argv,  float *result){
jamie@25: 
jamie@1:     int n = N;
jamie@1: 
jamie@1:     while(n--)
jamie@25: 	*result += (*data++ - ((float *)argv)[0]) / ((float *)argv)[1];
jamie@25: 
jamie@1:     *result = pow(*result, 4) / N - 3;
jamie@25: 
jamie@1: }
jamie@1: 
jamie@11: 
jamie@11: int xtract_centroid(float *data, int N, void *argv,  float *result){
jamie@25: 
jamie@37:     int n = (N >> 1);
jamie@11: 
jamie@37:     float *freqs, *amps, FA = 0.f, A = 0.f;
jamie@11: 
jamie@25:     freqs = data;
jamie@25:     amps = data + (N  >>  1);
jamie@25: 
jamie@11:     while(n--){
jamie@25: 	FA += freqs[n] * amps[n];
jamie@25: 	A += amps[n];
jamie@25:     }
jamie@25: 
jamie@25:     *result = FA / A;
jamie@11: 
jamie@11: }
jamie@11: 
jamie@1: int xtract_irregularity_k(float *data, int N, void *argv, float *result){
jamie@25: 
jamie@1:     int n,
jamie@37: 	M = N - 1;
jamie@1: 
jamie@1:     for(n = 1; n < M; n++)
jamie@25: 	*result += abs(data[n] - (data[n-1] + data[n] + data[n+1]) / 3);
jamie@1: 
jamie@1: }
jamie@1: 
jamie@1: int xtract_irregularity_j(float *data, int N, void *argv, float *result){
jamie@25: 
jamie@1:     int n = N;
jamie@1: 
jamie@37:     float num = 0.f, den = 0.f;
jamie@1: 
jamie@1:     while(n--){
jamie@25: 	num += data[n] - data[n+1];
jamie@25: 	den += data[n] * data[n];
jamie@1:     }
jamie@25: 
jamie@1:     *result = num / den;
jamie@1: 
jamie@1: }
jamie@1: 
jamie@1: int xtract_tristimulus_1(float *data, int N, void *argv, float *result){
jamie@1: 
jamie@1:     int n = N;
jamie@1: 
jamie@37:     float den = 0.f;
jamie@1: 
jamie@1:     while(n--)
jamie@25: 	den += data[n];
jamie@1: 
jamie@1:     *result = data[0] / den;
jamie@1: 
jamie@1: }
jamie@1: 
jamie@1: int xtract_tristimulus_2(float *data, int N, void *argv, float *result){
jamie@25: 
jamie@1:     int n = N;
jamie@1: 
jamie@37:     float den = 0.f;
jamie@1: 
jamie@1:     while(n--)
jamie@25: 	den += data[n];
jamie@1: 
jamie@1:     *result = (data[1] + data[2] + data[3])  / den;
jamie@25: 
jamie@1: }
jamie@1: 
jamie@1: int xtract_tristimulus_3(float *data, int N, void *argv, float *result){
jamie@25: 
jamie@1:     int n = N;
jamie@1: 
jamie@37:     float den = 0.f, num = 0.f;
jamie@1: 
jamie@1:     while(n--)
jamie@25: 	den += data[n];
jamie@1: 
jamie@1:     num = den - data[0] + data[1] + data[2] + data[3];
jamie@25: 
jamie@1:     *result = num / den;
jamie@25: 
jamie@1: }
jamie@1: 
jamie@1: int xtract_smoothness(float *data, int N, void *argv, float *result){
jamie@25: 
jamie@1:     int n = N;
jamie@1: 
jamie@1:     if (data[0] <= 0) data[0] = 1;
jamie@1:     if (data[1] <= 0) data[1] = 1;
jamie@25: 
jamie@1:     for(n = 2; n < N; n++){ 
jamie@25: 	if(data[n] <= 0) data[n] = 1;
jamie@25: 	*result += abs(20 * log(data[n-1]) - (20 * log(data[n-2]) + 
jamie@25: 		    20 * log(data[n-1]) + 20 * log(data[n])) / 3);
jamie@25:     }
jamie@1: }
jamie@1: 
jamie@1: int xtract_spread(float *data, int N, void *argv, float *result){
jamie@1: 
jamie@1:     int n = N;
jamie@1: 
jamie@37:     float num = 0.f, den = 0.f, tmp;
jamie@1: 
jamie@1:     while(n--){
jamie@25: 	tmp = n - *(float *)argv;
jamie@25: 	num += SQ(tmp) * data[n];
jamie@25: 	den += data[n];
jamie@1:     }
jamie@1: 
jamie@1:     *result = sqrt(num / den);
jamie@25: 
jamie@1: }
jamie@1: 
jamie@1: int xtract_zcr(float *data, int N, void *argv, float *result){
jamie@1: 
jamie@1:     int n = N;
jamie@25: 
jamie@1:     for(n = 1; n < N; n++)
jamie@25: 	if(data[n] * data[n-1] < 0) (*result)++;
jamie@25: 
jamie@1:     *result /= N;
jamie@25: 
jamie@1: }
jamie@1: 
jamie@1: int xtract_rolloff(float *data, int N, void *argv, float *result){
jamie@1: 
jamie@1:     int n = N;
jamie@37:     float pivot = 0.f, temp = 0.f;
jamie@25: 
jamie@1:     while(n--) pivot += data[n];   
jamie@25: 
jamie@1:     pivot *= *(float *)argv;
jamie@25: 
jamie@1:     for(n = 0; temp < pivot; temp += data[n++]);
jamie@1: 
jamie@1:     *result = n;
jamie@25: 
jamie@1: }
jamie@1: 
jamie@1: int xtract_loudness(float *data, int N, void *argv, float *result){
jamie@25: 
jamie@1:     int n = BARK_BANDS;
jamie@25: 
jamie@1:     /*if(n != N) return BAD_VECTOR_SIZE; */
jamie@1: 
jamie@1:     while(n--)
jamie@25: 	*result += pow(data[n], 0.23);
jamie@1: }
jamie@1: 
jamie@1: 
jamie@1: int xtract_flatness(float *data, int N, void *argv, float *result){
jamie@1: 
jamie@1:     int n = N;
jamie@1: 
jamie@37:     float num = 0.f, den = 0.f;
jamie@25: 
jamie@1:     while(n--){
jamie@25: 	if(data[n] !=0){
jamie@25: 	    num *= data[n];
jamie@25: 	    den += data[n];
jamie@25: 	}
jamie@1:     }
jamie@1: 
jamie@1:     num = pow(num, 1 / N); 
jamie@1:     den /= N;
jamie@25: 
jamie@1:     *result = 10 * log10(num / den);
jamie@25: 
jamie@1: }
jamie@1: 
jamie@1: int xtract_tonality(float *data, int N, void *argv, float *result){
jamie@25: 
jamie@1:     float sfmdb, sfm;
jamie@25: 
jamie@1:     sfm = *(float *)argv;
jamie@1: 
jamie@1:     sfmdb = (sfm > 0 ? (10 * log10(sfm)) / -60 : 0);
jamie@25: 
jamie@1:     *result = MIN(sfmdb, 1);
jamie@25: 
jamie@1: }
jamie@1: 
jamie@1: int xtract_crest(float *data, int N, void *argv, float *result){
jamie@25: 
jamie@25:     NOT_IMPLEMENTED;
jamie@25: 
jamie@1: }
jamie@1: 
jamie@1: int xtract_noisiness(float *data, int N, void *argv, float *result){
jamie@25: 
jamie@25:     NOT_IMPLEMENTED;
jamie@25: 
jamie@1: }
jamie@2: 
jamie@1: int xtract_rms_amplitude(float *data, int N, void *argv, float *result){
jamie@1: 
jamie@1:     int n = N;
jamie@1: 
jamie@1:     while(n--) *result += SQ(data[n]);
jamie@1: 
jamie@1:     *result = sqrt(*result / N);
jamie@25: 
jamie@1: }
jamie@1: 
jamie@1: int xtract_inharmonicity(float *data, int N, void *argv, float *result){
jamie@1: 
jamie@1:     int n = N;
jamie@37:     float num = 0.f, den = 0.f,
jamie@25: 	  *fund, *freq;
jamie@1: 
jamie@1:     fund = *(float **)argv;
jamie@1:     freq = fund+1;
jamie@25: 
jamie@1:     while(n--){
jamie@25: 	num += abs(freq[n] - n * *fund) * SQ(data[n]);
jamie@25: 	den += SQ(data[n]);
jamie@1:     }
jamie@1: 
jamie@1:     *result = (2 * num) / (*fund * den); 
jamie@25: 
jamie@1: }
jamie@1: 
jamie@1: 
jamie@1: int xtract_power(float *data, int N, void *argv, float *result){
jamie@1: 
jamie@25:     NOT_IMPLEMENTED;
jamie@25: 
jamie@1: }
jamie@1: 
jamie@1: int xtract_odd_even_ratio(float *data, int N, void *argv, float *result){
jamie@1: 
jamie@1:     int n = N >> 1, j, k;
jamie@1: 
jamie@37:     float num = 0.f, den = 0.f;
jamie@1: 
jamie@1:     while(n--){
jamie@25: 	j = n * 2;
jamie@25: 	k = j - 1;
jamie@25: 	num += data[k];
jamie@25: 	den += data[j];
jamie@1:     }
jamie@1: 
jamie@1:     *result = num / den;
jamie@25: 
jamie@1: }
jamie@1: 
jamie@1: int xtract_sharpness(float *data, int N, void *argv, float *result){
jamie@1: 
jamie@1:     NOT_IMPLEMENTED;
jamie@25: 
jamie@1: }
jamie@1: 
jamie@1: int xtract_slope(float *data, int N, void *argv, float *result){
jamie@1: 
jamie@1:     NOT_IMPLEMENTED;
jamie@25: 
jamie@1: }
jamie@1: 
jamie@5: int xtract_lowest_match(float *data, int N, void *argv, float *result){
jamie@25: 
jamie@22:     float lowest_match = SR_LIMIT;
jamie@1:     int n = N;
jamie@1: 
jamie@1:     while(n--) {
jamie@25: 	if(data[n] > 0)
jamie@25: 	    lowest_match = MIN(lowest_match, data[n]);
jamie@1:     }
jamie@1: 
jamie@22:     *result = (lowest_match == SR_LIMIT ? 0 : lowest_match);
jamie@25: 
jamie@1: }
jamie@1: 
jamie@1: int xtract_hps(float *data, int N, void *argv, float *result){
jamie@1: 
jamie@1:     int n = N, M, m, l, peak_index, position1_lwr;
jamie@1:     float *coeffs2, *coeffs3, *product, L, 
jamie@25: 	  largest1_lwr, peak, ratio1, sr;
jamie@1: 
jamie@25:     sr = *(float*)argv;
jamie@25: 
jamie@1:     coeffs2 = (float *)malloc(N * sizeof(float));
jamie@1:     coeffs3 = (float *)malloc(N * sizeof(float));
jamie@1:     product = (float *)malloc(N * sizeof(float));
jamie@25: 
jamie@1:     while(n--) coeffs2[n] = coeffs3[n] = 1;
jamie@1: 
jamie@1:     M = N >> 1;
jamie@1:     L = N / 3;
jamie@1: 
jamie@1:     while(M--){
jamie@25: 	m = M << 1;
jamie@25: 	coeffs2[M] = (data[m] + data[m+1]) * 0.5f;
jamie@1: 
jamie@25: 	if(M < L){
jamie@25: 	    l = M * 3;
jamie@25: 	    coeffs3[M] = (data[l] + data[l+1] + data[l+2]) / 3;
jamie@25: 	}
jamie@1:     }
jamie@25: 
jamie@1:     peak_index = peak = 0;
jamie@25: 
jamie@1:     for(n = 1; n < N; n++){
jamie@25: 	product[n] = data[n] * coeffs2[n] * coeffs3[n];
jamie@25: 	if(product[n] > peak){
jamie@25: 	    peak_index = n;
jamie@25: 	    peak = product[n];
jamie@25: 	}
jamie@1:     }
jamie@1: 
jamie@1:     largest1_lwr = position1_lwr = 0;
jamie@1: 
jamie@1:     for(n = 0; n < N; n++){
jamie@25: 	if(data[n] > largest1_lwr && n != peak_index){
jamie@25: 	    largest1_lwr = data[n];
jamie@25: 	    position1_lwr = n;
jamie@25: 	}
jamie@1:     }
jamie@1: 
jamie@1:     ratio1 = data[position1_lwr] / data[peak_index];
jamie@1: 
jamie@1:     if(position1_lwr > peak_index * 0.4 && position1_lwr < 
jamie@25: 	    peak_index * 0.6 && ratio1 > 0.1)
jamie@25: 	peak_index = position1_lwr;
jamie@1: 
jamie@22:     *result = sr / (float)peak_index; 
jamie@25: 
jamie@1:     free(coeffs2);
jamie@1:     free(coeffs3);
jamie@1:     free(product);
jamie@25: 
jamie@1: }
jamie@5: 
jamie@5: 
jamie@5: int xtract_f0(float *data, int N, void *argv, float *result){
jamie@5: 
jamie@25:     int M, sr, tau, n;
jamie@25:     float f0, err_tau_1, err_tau_x, array_max, threshold_peak, threshold_centre;
jamie@22: 
jamie@25:     sr = *(float *)argv;
jamie@25:     /*  threshold_peak = *((float *)argv+1);
jamie@25: 	threshold_centre = *((float *)argv+2);
jamie@25: 	printf("peak: %.2f\tcentre: %.2f\n", threshold_peak, threshold_centre);*/
jamie@25:     /* add temporary dynamic control over thresholds to test clipping effects */
jamie@22: 
jamie@25:     /* FIX: tweak and  make into macros */
jamie@25:     threshold_peak = .8;
jamie@25:     threshold_centre = .3;
jamie@25:     M = N >> 1;
jamie@25:     err_tau_1 = 0;
jamie@25:     array_max = 0;
jamie@25: 
jamie@25:     /* Find the array max */
jamie@25:     for(n = 0; n < N; n++){
jamie@25: 	if (data[n] > array_max)
jamie@25: 	    array_max = data[n];
jamie@12:     }
jamie@25: 
jamie@25:     threshold_peak *= array_max;
jamie@25: 
jamie@25:     /* peak clip */
jamie@25:     for(n = 0; n < N; n++){
jamie@25: 	if(data[n] > threshold_peak)
jamie@25: 	    data[n] = threshold_peak;
jamie@25: 	else if(data[n] < -threshold_peak)
jamie@25: 	    data[n] = -threshold_peak;
jamie@25:     }
jamie@25: 
jamie@25:     threshold_centre *= array_max;
jamie@25: 
jamie@25:     /* Centre clip */
jamie@25:     for(n = 0; n < N; n++){
jamie@25: 	if (data[n] < threshold_centre)
jamie@25: 	    data[n] = 0;
jamie@25: 	else 
jamie@25: 	    data[n] -= threshold_centre;
jamie@25:     }
jamie@25: 
jamie@25:     /* Estimate fundamental freq */
jamie@25:     for (n = 1; n < M; n++)
jamie@25: 	err_tau_1 = err_tau_1 + fabs(data[n] - data[n+1]);
jamie@25:     /* FIX: this doesn't pose too much load if it returns 'early', but if it can't find f0, load can be significant for larger block sizes M^2 iterations! */  
jamie@25:     for (tau = 2; tau < M; tau++){
jamie@25: 	err_tau_x = 0;
jamie@25: 	for (n = 1; n < M; n++){
jamie@25: 	    err_tau_x = err_tau_x + fabs(data[n] - data[n+tau]);
jamie@25: 	}
jamie@25: 	if (err_tau_x < err_tau_1) {
jamie@25: 	    f0 = sr / (tau + (err_tau_x / err_tau_1));
jamie@25: 	    *result = f0;
jamie@25: 	    return SUCCESS;
jamie@25: 	}
jamie@25:     }
jamie@25:     return NO_RESULT;
jamie@5: }