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@56: #include "xtract_macros_private.h"
jamie@1: #include "math.h"
jamie@5: #include <stdlib.h>
jamie@43: #include <string.h>
jamie@78: #include <stdio.h>
jamie@1: 
jamie@85: #ifndef powf
jamie@85:     #define powf pow
jamie@85: #endif
jamie@85: 
jamie@85: #ifndef expf
jamie@85:     #define expf exp
jamie@85: #endif
jamie@85: 
jamie@88: int test(void){
jamie@88:     printf("Hello world");
jamie@88: }
jamie@88: 
jamie@43: int xtract_mean(const float *data, const int N, const void *argv, float *result){
jamie@25: 
jamie@1:     int n = N;
jamie@87:     
jamie@87:     *result = 0.f;
jamie@1: 
jamie@1:     while(n--)
jamie@42: 	*result += data[n];
jamie@87:     
jamie@1:     *result /= N;
jamie@38: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@43: int xtract_variance(const float *data, const int N, const void *argv, float *result){
jamie@25: 
jamie@1:     int n = N;
jamie@1: 
jamie@1:     while(n--)
jamie@43: 	*result += pow(data[n] - *(float *)argv, 2);
jamie@25: 
jamie@43:     *result = *result / (N - 1);
jamie@44: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@43: int xtract_standard_deviation(const float *data, const int N, const void *argv, float *result){
jamie@25: 
jamie@61:     *result = sqrt(*(float *)argv);
jamie@25: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@43: int xtract_average_deviation(const float *data, const int N, const void *argv, float *result){
jamie@25: 
jamie@1:     int n = N;
jamie@44: 
jamie@1:     while(n--)
jamie@42: 	*result += fabs(data[n] - *(float *)argv);
jamie@25: 
jamie@1:     *result /= N;
jamie@1: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@43: int xtract_skewness(const float *data, const int N, const void *argv,  float *result){
jamie@25: 
jamie@1:     int n = N;
jamie@1: 
jamie@59:     float temp = 0.f;
jamie@25: 
jamie@42:     while(n--){
jamie@42: 	temp = (data[n] - ((float *)argv)[0]) / ((float *)argv)[1];
jamie@42: 	*result += pow(temp, 3);
jamie@42:     }
jamie@1: 
jamie@42:     *result /= N;
jamie@44: 
jamie@59: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@43: int xtract_kurtosis(const float *data, const int N, const void *argv,  float *result){
jamie@25: 
jamie@1:     int n = N;
jamie@1: 
jamie@42:     float temp;
jamie@25: 
jamie@42:     while(n--){
jamie@42: 	temp = (data[n] - ((float *)argv)[0]) / ((float *)argv)[1];
jamie@42: 	*result += pow(temp, 4);
jamie@42:     }
jamie@25: 
jamie@42:     *result /= N;
jamie@42:     *result -= 3.0f;
jamie@44: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@52: int xtract_spectral_centroid(const float *data, const int N, const void *argv,  float *result){
jamie@25: 
jamie@37:     int n = (N >> 1);
jamie@11: 
jamie@43:     const float *freqs, *amps;
jamie@43:     float FA = 0.f, A = 0.f;
jamie@11: 
jamie@52:     amps = data;
jamie@52:     freqs = data + n;
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@56:     return XTRACT_SUCCESS;
jamie@11: }
jamie@11: 
jamie@52: int xtract_spectral_mean(const float *data, const int N, const void *argv, float *result){
jamie@52: 
jamie@52:     return xtract_spectral_centroid(data, N, argv, result);
jamie@52: 
jamie@52: }
jamie@52: 
jamie@52: int xtract_spectral_variance(const float *data, const int N, const void *argv, float *result){
jamie@52: 
jamie@53:     int m;
jamie@53:     float A = 0.f;
jamie@53:     const float *freqs, *amps;
jamie@52: 
jamie@53:     m = N >> 1;
jamie@52: 
jamie@53:     amps = data;
jamie@53:     freqs = data + m;
jamie@52: 
jamie@53:     while(m--){
jamie@53: 	A += amps[m];
jamie@59: 	*result += powf((freqs[m]  - *(float *)argv) * amps[m], 2);
jamie@53:     }
jamie@53: 
jamie@59:     *result = *result / (A /*- 1*/);
jamie@52: 
jamie@56:     return XTRACT_SUCCESS;
jamie@52: }
jamie@52: 
jamie@52: int xtract_spectral_standard_deviation(const float *data, const int N, const void *argv, float *result){
jamie@52: 
jamie@59:     *result = sqrt(*(float *)argv); 
jamie@52: 
jamie@56:     return XTRACT_SUCCESS;
jamie@52: }
jamie@52: 
jamie@52: int xtract_spectral_average_deviation(const float *data, const int N, const void *argv, float *result){
jamie@52: 
jamie@53:     int m;
jamie@53:     float A = 0.f;
jamie@53:     const float *freqs, *amps;
jamie@52: 
jamie@53:     m = N >> 1;
jamie@52: 
jamie@53:     amps = data;
jamie@53:     freqs = data + m;
jamie@52: 
jamie@53:     while(m--){
jamie@53: 	A += amps[m];
jamie@53: 	*result += fabs((amps[m] * freqs[m]) - *(float *)argv);
jamie@53:     }
jamie@53: 
jamie@53:     *result /= A;
jamie@52: 
jamie@56:     return XTRACT_SUCCESS;
jamie@52: }
jamie@52: 
jamie@52: int xtract_spectral_skewness(const float *data, const int N, const void *argv,  float *result){
jamie@52: 
jamie@53:     int m;
jamie@53:     float temp, A = 0.f;
jamie@53:     const float *freqs, *amps;
jamie@52: 
jamie@53:     m = N >> 1;
jamie@53: 
jamie@53:     amps = data;
jamie@53:     freqs = data + m;
jamie@52: 
jamie@52:     while(m--){
jamie@53: 	A += amps[m];
jamie@53: 	temp = ((amps[m] * freqs[m]) - 
jamie@52: 		((float *)argv)[0]) / ((float *)argv)[1];
jamie@52: 	*result += pow(temp, 3);
jamie@52:     }
jamie@52: 
jamie@53:     *result /= A;
jamie@52: 
jamie@56:     return XTRACT_SUCCESS;
jamie@52: }
jamie@52: 
jamie@52: int xtract_spectral_kurtosis(const float *data, const int N, const void *argv,  float *result){
jamie@52: 
jamie@53:     int m;
jamie@53:     float temp, A = 0.f;
jamie@53:     const float *freqs, *amps;
jamie@52: 
jamie@53:     m = N >> 1;
jamie@53: 
jamie@53:     amps = data;
jamie@53:     freqs = data + m;
jamie@52: 
jamie@52:     while(m--){
jamie@53: 	A += amps[m];
jamie@53: 	temp = ((amps[m] * freqs[m]) - 
jamie@52: 		((float *)argv)[0]) / ((float *)argv)[1];
jamie@52: 	*result += pow(temp, 4);
jamie@52:     }
jamie@52: 
jamie@53:     *result /= A;
jamie@52:     *result -= 3.0f;
jamie@52: 
jamie@56:     return XTRACT_SUCCESS;
jamie@52: }
jamie@52: 
jamie@43: int xtract_irregularity_k(const float *data, const int N, const void *argv, float *result){
jamie@25: 
jamie@1:     int n,
jamie@37: 	M = N - 1;
danstowell@84: 	
danstowell@84: 	*result = 0.f;
danstowell@84: 	
jamie@1:     for(n = 1; n < M; n++)
jamie@42: 	*result += fabs(data[n] - (data[n-1] + data[n] + data[n+1]) / 3);
jamie@1: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@43: int xtract_irregularity_j(const float *data, const int N, const void *argv, float *result){
jamie@25: 
jamie@1:     int n = N;
jamie@1: 
jamie@59:     double num = 0.f, den = 0.f;
jamie@1: 
jamie@1:     while(n--){
jamie@42: 	num += pow(data[n] - data[n+1], 2);
jamie@42: 	den += pow(data[n], 2);
jamie@1:     }
jamie@25: 
jamie@59:     *result = (float)(num / den);
jamie@1: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@43: int xtract_tristimulus_1(const float *data, const int N, const void *argv, float *result){
jamie@1: 
jamie@1:     int n = N;
jamie@1: 
jamie@42:     float den, p1, temp;
jamie@1: 
jamie@42:     den = p1 = temp = 0.f;
jamie@1: 
jamie@42:     for(n = 0; n < N; n++){
jamie@42: 	if((temp = data[n])){
jamie@42: 	    den += temp;
jamie@42: 	    if(!p1)
jamie@42: 		p1 = temp;
jamie@42: 	}
jamie@42:     }
jamie@42: 
jamie@42:     *result = p1 / den;
jamie@1: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@43: int xtract_tristimulus_2(const float *data, const int N, const void *argv, float *result){
jamie@25: 
jamie@1:     int n = N;
jamie@1: 
jamie@42:     float den, p2, p3, p4, temp;
jamie@1: 
jamie@42:     den = p2 = p3 = p4 = temp = 0.f;
jamie@1: 
jamie@42:     for(n = 0; n < N; n++){
jamie@42: 	if((temp = data[n])){
jamie@42: 	    den += temp;
jamie@42: 	    if(!p2)
jamie@42: 		p2 = temp;
jamie@42: 	    else if(!p3)
jamie@42: 		p3 = temp;
jamie@42: 	    else if(!p4)
jamie@42: 		p4 = temp;
jamie@42: 	}
jamie@42:     }
jamie@42: 
jamie@42:     *result = (p2 + p3 + p4)  / den;
jamie@25: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@43: int xtract_tristimulus_3(const float *data, const int N, const void *argv, float *result){
jamie@25: 
jamie@42:     int n = N, count = 0;
jamie@1: 
jamie@42:     float den, num, temp;
jamie@1: 
jamie@42:     den = num = temp = 0.f;
jamie@1: 
jamie@42:     for(n = 0; n < N; n++){
jamie@42: 	if((temp = data[n])){
jamie@42: 	    den += temp;
jamie@42: 	    if(count >= 5)
jamie@42: 		num += temp;
jamie@42: 	    count++;
jamie@42: 	}
jamie@42:     }
jamie@25: 
jamie@1:     *result = num / den;
jamie@25: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@43: int xtract_smoothness(const float *data, const int N, const void *argv, float *result){
jamie@25: 
jamie@59:     int n, M;
jamie@1: 
jamie@43:     float *input;
jamie@43: 
jamie@43:     input = (float *)malloc(N * sizeof(float));
jamie@59:     memcpy(input, data, N * sizeof(float));
jamie@43: 
jamie@43:     if (input[0] <= 0) input[0] = 1;
jamie@25: 
jamie@59:     M = N - 1;
jamie@59: 
jamie@59:     for(n = 1; n < M; n++){ 
jamie@43: 	if(input[n] <= 0) input[n] = 1;
jamie@59: 	*result += fabs(20 * log(input[n]) - (20 * log(input[n-1]) + 
jamie@59: 		    20 * log(input[n]) + 20 * log(input[n+1])) / 3);
jamie@25:     }
jamie@43: 
jamie@43:     free(input);
jamie@44: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@43: int xtract_spread(const float *data, const int N, const void *argv, float *result){
jamie@1: 
jamie@1:     int n = N;
jamie@1: 
jamie@83:     float num = 0.f, den = 0.f, temp = 0.f;
jamie@83: 
jamie@83:     if(argv == NULL)
jamie@83: 	return XTRACT_BAD_ARGV;
jamie@1: 
jamie@1:     while(n--){
jamie@48: 	temp = n - *(float *)argv;
jamie@56: 	num += XTRACT_SQ(temp) * data[n];
jamie@25: 	den += data[n];
jamie@1:     }
jamie@1: 
jamie@1:     *result = sqrt(num / den);
jamie@25: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@43: int xtract_zcr(const float *data, const int N, const 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@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@43: int xtract_rolloff(const float *data, const int N, const void *argv, float *result){
jamie@1: 
jamie@1:     int n = N;
jamie@55:     float pivot, temp, percentile;
jamie@42: 
jamie@42:     pivot = temp = 0.f;
jamie@55:     percentile = ((float *)argv)[1];
jamie@25: 
jamie@1:     while(n--) pivot += data[n];   
jamie@25: 
jamie@55:     pivot *= percentile / 100.f;
jamie@25: 
jamie@42:     for(n = 0; temp < pivot; n++)
jamie@42: 	temp += data[n];
jamie@1: 
jamie@55:     *result = n * ((float *)argv)[0];
jamie@55:     /* *result = (n / (float)N) * (((float *)argv)[1] * .5); */
jamie@25: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@43: int xtract_loudness(const float *data, const int N, const void *argv, float *result){
jamie@25: 
jamie@47:     int n = N, rv;
jamie@25: 
jamie@56:     if(n > XTRACT_BARK_BANDS) 
jamie@56: 	rv = XTRACT_BAD_VECTOR_SIZE; 
jamie@47:     else
jamie@56: 	rv = XTRACT_SUCCESS;
jamie@1: 
jamie@1:     while(n--)
jamie@59: 	*result += powf(data[n], 0.23);
jamie@38: 
jamie@47:     return rv;
jamie@1: }
jamie@1: 
jamie@43: int xtract_flatness(const float *data, const int N, const void *argv, float *result){
jamie@1: 
jamie@42:     int n;
jamie@1: 
jamie@44:     double num, den, temp; 
jamie@25: 
jamie@44:     den = data[0];
jamie@44:     num = (data[0] == 0.f ? 1.f : data[0]);
jamie@43: 
jamie@44:     for(n = 1; n < N; n++){
jamie@44: 	if((temp = data[n]) != 0.f) {
jamie@44: 	    num *= temp;
jamie@44: 	    den += temp;
jamie@25: 	}
jamie@1:     }
jamie@44: 
jamie@59:     num = powf(num, 1.f / N);
jamie@1:     den /= N;
jamie@25: 
jamie@56:     if(num < XTRACT_VERY_SMALL_NUMBER) 
jamie@56: 	num = XTRACT_VERY_SMALL_NUMBER;
jamie@44: 
jamie@56:     if(den < XTRACT_VERY_SMALL_NUMBER) 
jamie@56: 	den = XTRACT_VERY_SMALL_NUMBER;
jamie@44: 
jamie@59:     *result = 10 * log10(num / den);
jamie@25: 
jamie@56:     return XTRACT_SUCCESS;
jamie@44: 
jamie@1: }
jamie@1: 
jamie@43: int xtract_tonality(const float *data, const int N, const void *argv, float *result){
jamie@25: 
jamie@1:     float sfmdb, sfm;
jamie@25: 
jamie@1:     sfm = *(float *)argv;
jamie@1: 
jamie@59:     sfmdb = sfm / -60.f;
jamie@25: 
jamie@56:     *result = XTRACT_MIN(sfmdb, 1);
jamie@25: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@43: int xtract_crest(const float *data, const int N, const void *argv, float *result){
jamie@25: 
jamie@45:     float max, mean; 
jamie@45: 
jamie@45:     max = mean = 0.f;
jamie@45: 
jamie@45:     max = *(float *)argv;
jamie@45:     mean = *((float *)argv+1);
jamie@45: 
jamie@45:     *result = max / mean;
jamie@45: 
jamie@56:     return XTRACT_SUCCESS;
jamie@25: 
jamie@1: }
jamie@1: 
jamie@43: int xtract_noisiness(const float *data, const int N, const void *argv, float *result){
jamie@25: 
jamie@45:     float h, i, p; /*harmonics, inharmonics, partials */
jamie@45: 
jamie@45:     i = p = h = 0.f;
jamie@45: 
jamie@45:     h = *(float *)argv;
jamie@45:     p = *((float *)argv+1);
jamie@45: 
jamie@45:     i = p - h;
jamie@45: 
jamie@45:     *result = i / p;
jamie@45: 
jamie@56:     return XTRACT_SUCCESS;
jamie@25: 
jamie@1: }
jamie@2: 
jamie@43: int xtract_rms_amplitude(const float *data, const int N, const void *argv, float *result){
jamie@1: 
jamie@1:     int n = N;
jamie@1: 
jamie@56:     while(n--) *result += XTRACT_SQ(data[n]);
jamie@1: 
jamie@1:     *result = sqrt(*result / N);
jamie@25: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@52: int xtract_spectral_inharmonicity(const float *data, const int N, const void *argv, float *result){
jamie@1: 
jamie@41:     int n = N >> 1;
jamie@43:     float num = 0.f, den = 0.f, fund; 
jamie@43:     const float *freqs, *amps;
jamie@1: 
jamie@41:     fund = *(float *)argv;
jamie@52:     amps = data;
jamie@52:     freqs = data + n;
jamie@25: 
jamie@1:     while(n--){
jamie@59: 	if(amps[n]){
jamie@59: 	    num += fabs(freqs[n] - n * fund) * XTRACT_SQ(amps[n]);
jamie@59: 	    den += XTRACT_SQ(amps[n]);
jamie@59: 	}
jamie@1:     }
jamie@1: 
jamie@41:     *result = (2 * num) / (fund * den); 
jamie@25: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@1: 
jamie@43: int xtract_power(const float *data, const int N, const void *argv, float *result){
jamie@1: 
jamie@56:     return XTRACT_FEATURE_NOT_IMPLEMENTED;
jamie@25: 
jamie@1: }
jamie@1: 
jamie@43: int xtract_odd_even_ratio(const float *data, const int N, const void *argv, float *result){
jamie@1: 
jamie@43:     int M = (N >> 1), n;
jamie@1: 
jamie@59:     float num = 0.f, den = 0.f,  temp;
jamie@44: 
jamie@43:     for(n = 0; n < M; n++){
jamie@43: 	if((temp = data[n])){
jamie@59: 	    if(XTRACT_IS_ODD(n)){
jamie@59: 		num += temp;
jamie@43: 	    }
jamie@43: 	    else{
jamie@59: 		den += temp;
jamie@43: 	    }
jamie@43: 	}
jamie@1:     }
jamie@1: 
jamie@1:     *result = num / den;
jamie@25: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@43: int xtract_sharpness(const float *data, const int N, const void *argv, float *result){
jamie@1: 
jamie@48:     int n = N, rv;
jamie@48:     float sl, g, temp; /* sl = specific loudness */
jamie@48: 
jamie@48:     sl = g = temp = 0.f;
jamie@48: 
jamie@56:     if(n > XTRACT_BARK_BANDS) 
jamie@56: 	rv = XTRACT_BAD_VECTOR_SIZE; 
jamie@48:     else
jamie@56: 	rv = XTRACT_SUCCESS;
jamie@48: 
jamie@48: 
jamie@48:     while(n--){
jamie@59: 	sl = powf(data[n], 0.23);
jamie@59: 	g = (n < 15 ? 1.f : 0.066 * expf(0.171 * n));
jamie@49: 	temp += n * g * sl;
jamie@48:     }
jamie@48: 
jamie@48:     *result = 0.11 * temp / N;
jamie@48: 
jamie@48:     return rv;
jamie@25: 
jamie@1: }
jamie@1: 
jamie@52: int xtract_spectral_slope(const float *data, const int N, const void *argv, float *result){
jamie@1: 
jamie@48:     const float *freqs, *amps; 
jamie@48:     float f, a,
jamie@56: 	  F, A, FA, FXTRACT_SQ; /* sums of freqs, amps, freq * amps, freq squared */
jamie@48:     int n, M; 
jamie@48:     
jamie@56:     F = A = FA = FXTRACT_SQ = 0.f;
jamie@48:     n = M = N >> 1;
jamie@48: 
jamie@52:     amps = data;
jamie@52:     freqs = data + n;
jamie@48: 
jamie@48:     while(n--){
jamie@48: 	f = freqs[n];
jamie@48: 	a = amps[n];
jamie@48: 	F += f;
jamie@48: 	A += a;
jamie@48: 	FA += f * a;
jamie@56: 	FXTRACT_SQ += f * f;
jamie@48:     }
jamie@48: 
jamie@56:     *result = (1.f / A) * (M * FA - F * A) / (M * FXTRACT_SQ - F * F); 
jamie@48: 
jamie@56:     return XTRACT_SUCCESS;
jamie@25: 
jamie@1: }
jamie@1: 
jamie@45: int xtract_lowest_value(const float *data, const int N, const void *argv, float *result){
jamie@25: 
jamie@45:     int n = N;
jamie@45:     float temp;
jamie@45: 
jamie@46:     *result = data[--n];
jamie@45: 
jamie@45:     while(n--){
jamie@45:        if((temp = data[n]) > *(float *)argv)	
jamie@56: 	    *result = XTRACT_MIN(*result, data[n]);
jamie@45:     }
jamie@45: 
jamie@56:     return XTRACT_SUCCESS;
jamie@45: }
jamie@45: 
jamie@45: int xtract_highest_value(const float *data, const int N, const void *argv, float *result){
jamie@45: 
jamie@1:     int n = N;
jamie@1: 
jamie@46:     *result = data[--n];
jamie@44: 
jamie@45:     while(n--) 
jamie@56: 	*result = XTRACT_MAX(*result, data[n]);
jamie@44: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@45: 
jamie@45: int xtract_sum(const float *data, const int N, const void *argv, float *result){
jamie@45: 
jamie@45:     int n = N;
jamie@45: 
jamie@45:     while(n--)
jamie@45: 	*result += *data++;
jamie@45: 
jamie@56:     return XTRACT_SUCCESS;
jamie@45: 
jamie@45: }
jamie@45: 
jamie@59: int xtract_nonzero_count(const float *data, const int N, const void *argv, float *result){
jamie@59: 
jamie@59:     int n = N;
jamie@59: 
jamie@59:     while(n--)
jamie@59: 	*result += (*data++ ? 1 : 0);
jamie@59: 
jamie@59:     return XTRACT_SUCCESS;
jamie@59: 
jamie@59: }
jamie@59: 
jamie@43: int xtract_hps(const float *data, const int N, const 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@78:     if(sr == 0)
jamie@78: 	sr = 44100.f;
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@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@5: 
jamie@5: 
jamie@43: int xtract_f0(const float *data, const int N, const void *argv, float *result){
jamie@5: 
jamie@78:     int M, tau, n;
jamie@78:     float sr;
jamie@43:     size_t bytes;
jamie@43:     float f0, err_tau_1, err_tau_x, array_max, 
jamie@43: 	  threshold_peak, threshold_centre,
jamie@43: 	  *input;
jamie@22: 
jamie@25:     sr = *(float *)argv;
jamie@78:     if(sr == 0)
jamie@78: 	sr = 44100.f;
jamie@43: 
jamie@43:     input = (float *)malloc(bytes = N * sizeof(float));
jamie@43:     input = memcpy(input, data, bytes);
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@43: 	if (input[n] > array_max)
jamie@43: 	    array_max = input[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@43: 	if(input[n] > threshold_peak)
jamie@43: 	    input[n] = threshold_peak;
jamie@43: 	else if(input[n] < -threshold_peak)
jamie@43: 	    input[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@43: 	if (input[n] < threshold_centre)
jamie@43: 	    input[n] = 0;
jamie@25: 	else 
jamie@43: 	    input[n] -= threshold_centre;
jamie@25:     }
jamie@25: 
jamie@25:     /* Estimate fundamental freq */
jamie@25:     for (n = 1; n < M; n++)
jamie@43: 	err_tau_1 = err_tau_1 + fabs(input[n] - input[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@43: 	    err_tau_x = err_tau_x + fabs(input[n] - input[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@43: 	    free(input);
jamie@56: 	    return XTRACT_SUCCESS;
jamie@25: 	}
jamie@25:     }
jamie@43:     *result = -0;
jamie@43:     free(input);
jamie@56:     return XTRACT_NO_RESULT;
jamie@5: }
jamie@43: 
jamie@43: int xtract_failsafe_f0(const float *data, const int N, const void *argv, float *result){
jamie@44: 
jamie@59:     float *spectrum = NULL, argf[2], *peaks = NULL, return_code, sr;
jamie@44: 
jamie@43:     return_code = xtract_f0(data, N, argv, result);
jamie@44: 
jamie@56:     if(return_code == XTRACT_NO_RESULT){
jamie@44: 
jamie@59: 	sr = *(float *)argv;
jamie@78: 	if(sr == 0)
jamie@78: 	    sr = 44100.f;
jamie@59: 	spectrum = (float *)malloc(N * sizeof(float));
jamie@43: 	peaks = (float *)malloc(N * sizeof(float));
jamie@59: 	argf[0] = sr;
jamie@59: 	argf[1] = XTRACT_MAGNITUDE_SPECTRUM;
jamie@59: 	xtract_spectrum(data, N, argf, spectrum);
jamie@59: 	argf[1] = 10.f;
jamie@59: 	xtract_peak_spectrum(spectrum, N >> 1, argf, peaks);
jamie@43: 	argf[0] = 0.f;
jamie@59: 	xtract_lowest_value(peaks+(N >> 1), N >> 1, argf, result);
jamie@44: 
jamie@59: 	free(spectrum);
jamie@43: 	free(peaks);
jamie@43:     }
jamie@43: 
jamie@56:     return XTRACT_SUCCESS;
jamie@43: 
jamie@43: }
jamie@44: