jamie@141: /*
jamie@141:  * Copyright (C) 2012 Jamie Bullock
jamie@140:  *
jamie@141:  * Permission is hereby granted, free of charge, to any person obtaining a copy
jamie@141:  * of this software and associated documentation files (the "Software"), to
jamie@141:  * deal in the Software without restriction, including without limitation the
jamie@141:  * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
jamie@141:  * sell copies of the Software, and to permit persons to whom the Software is
jamie@141:  * furnished to do so, subject to the following conditions:
jamie@1:  *
jamie@141:  * The above copyright notice and this permission notice shall be included in
jamie@141:  * all copies or substantial portions of the Software.
jamie@1:  *
jamie@141:  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
jamie@141:  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
jamie@141:  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
jamie@141:  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
jamie@141:  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
jamie@141:  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
jamie@141:  * IN THE SOFTWARE.
jamie@1:  *
jamie@1:  */
jamie@1: 
jamie@107: /* scalar.c: defines functions that extract a feature as a single value from an input vector */
jamie@1: 
jamie@5: #include <stdlib.h>
jamie@43: #include <string.h>
jamie@78: #include <stdio.h>
jamie@113: #include <math.h>
jamie@192: #include <limits.h>
jamie@232: #include <stdint.h>
jamie@113: 
jamie@195: #ifndef DBL_MAX
jamie@195: #include <float.h> /* on Linux DBL_MAX is in float.h */
jamie@195: #endif
jamie@195: 
jamie@161: #include "dywapitchtrack/dywapitchtrack.h"
jamie@161: 
jamie@259: #include "xtract/libxtract.h"
jamie@259: #include "xtract/xtract_helper.h"
jamie@113: #include "xtract_macros_private.h"
jamie@161: #include "xtract_globals_private.h"
jamie@1: 
jamie@146: int xtract_mean(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@25: 
jamie@1:     int n = N;
jamie@1: 
jamie@146:     *result = 0.0;
jamie@113: 
jamie@1:     while(n--)
jamie@140:         *result += data[n];
jamie@140: 
jamie@140:     *result /= N;
jamie@140: 
jamie@140:     return XTRACT_SUCCESS;
jamie@140: }
jamie@140: 
jamie@146: int xtract_variance(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@140: 
jamie@140:     int n = N;
jamie@140: 
jamie@146:     *result = 0.0;
jamie@140: 
jamie@140:     while(n--)
jamie@146:         *result += pow(data[n] - *(double *)argv, 2);
jamie@25: 
jamie@43:     *result = *result / (N - 1);
jamie@44: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@146: int xtract_standard_deviation(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@25: 
jamie@146:     *result = sqrt(*(double *)argv);
jamie@25: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@146: int xtract_average_deviation(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@25: 
jamie@1:     int n = N;
jamie@44: 
jamie@146:     *result = 0.0;
jamie@113: 
jamie@1:     while(n--)
jamie@146:         *result += fabs(data[n] - *(double *)argv);
jamie@25: 
jamie@1:     *result /= N;
jamie@1: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@146: int xtract_skewness(const double *data, const int N, const void *argv,  double *result)
jamie@140: {
jamie@25: 
jamie@1:     int n = N;
jamie@1: 
jamie@146:     double temp = 0.0;
jamie@25: 
jamie@146:     *result = 0.0;
jamie@113: 
jamie@140:     while(n--)
jamie@140:     {
jamie@146:         temp = (data[n] - ((double *)argv)[0]) / ((double *)argv)[1];
jamie@146:         *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@146: int xtract_kurtosis(const double *data, const int N, const void *argv,  double *result)
jamie@140: {
jamie@25: 
jamie@1:     int n = N;
jamie@1: 
jamie@146:     double temp = 0.0;
jamie@113: 
jamie@146:     *result = 0.0;
jamie@25: 
jamie@140:     while(n--)
jamie@140:     {
jamie@146:         temp = (data[n] - ((double *)argv)[0]) / ((double *)argv)[1];
jamie@146:         *result += pow(temp, 4);
jamie@42:     }
jamie@25: 
jamie@42:     *result /= N;
jamie@146:     *result -= 3.0;
jamie@44: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@146: int xtract_spectral_centroid(const double *data, const int N, const void *argv,  double *result)
jamie@140: {
jamie@25: 
jamie@37:     int n = (N >> 1);
jamie@11: 
jamie@146:     const double *freqs, *amps;
jamie@146:     double FA = 0.0, A = 0.0;
jamie@11: 
jamie@52:     amps = data;
jamie@52:     freqs = data + n;
jamie@25: 
jamie@140:     while(n--)
jamie@140:     {
jamie@140:         FA += freqs[n] * amps[n];
jamie@140:         A += amps[n];
jamie@25:     }
jamie@25: 
jamie@146:     if(A == 0.0)
jamie@146:         *result = 0.0;
jamie@113:     else
jamie@113:         *result = FA / A;
jamie@11: 
jamie@56:     return XTRACT_SUCCESS;
jamie@11: }
jamie@11: 
jamie@146: int xtract_spectral_mean(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@52: 
jamie@52:     return xtract_spectral_centroid(data, N, argv, result);
jamie@52: 
jamie@52: }
jamie@52: 
jamie@146: int xtract_spectral_variance(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@52: 
jamie@53:     int m;
jamie@146:     double A = 0.0;
jamie@146:     const double *freqs, *amps;
jamie@52: 
jamie@53:     m = N >> 1;
jamie@52: 
jamie@53:     amps = data;
jamie@53:     freqs = data + m;
jamie@52: 
jamie@146:     *result = 0.0;
jamie@113: 
jamie@140:     while(m--)
jamie@140:     {
jamie@123:         A += amps[m];
jamie@146:         *result += pow(freqs[m] - ((double *)argv)[0], 2) * amps[m];
jamie@53:     }
jamie@53: 
jamie@123:     *result = *result / A;
jamie@52: 
jamie@56:     return XTRACT_SUCCESS;
jamie@52: }
jamie@52: 
jamie@146: int xtract_spectral_standard_deviation(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@52: 
jamie@146:     *result = sqrt(*(double *)argv);
jamie@52: 
jamie@56:     return XTRACT_SUCCESS;
jamie@52: }
jamie@52: 
jamie@146: /*int xtract_spectral_average_deviation(const double *data, const int N, const void *argv, double *result){
jamie@52: 
jamie@53:     int m;
jamie@146:     double A = 0.0;
jamie@146:     const double *freqs, *amps;
jamie@52: 
jamie@53:     m = N >> 1;
jamie@52: 
jamie@53:     amps = data;
jamie@53:     freqs = data + m;
jamie@52: 
jamie@146:     *result = 0.0;
jamie@113: 
jamie@53:     while(m--){
jamie@123:         A += amps[m];
jamie@146:         *result += fabs((amps[m] * freqs[m]) - *(double *)argv);
jamie@53:     }
jamie@53: 
jamie@53:     *result /= A;
jamie@52: 
jamie@56:     return XTRACT_SUCCESS;
jamie@123: }*/
jamie@52: 
jamie@146: int xtract_spectral_skewness(const double *data, const int N, const void *argv,  double *result)
jamie@140: {
jamie@52: 
jamie@53:     int m;
jamie@146:     const double *freqs, *amps;
jamie@52: 
jamie@53:     m = N >> 1;
jamie@53: 
jamie@53:     amps = data;
jamie@53:     freqs = data + m;
jamie@52: 
jamie@146:     *result = 0.0;
jamie@113: 
jamie@123:     while(m--)
jamie@146:         *result += pow(freqs[m] - ((double *)argv)[0], 3) * amps[m];
jamie@52: 
jamie@146:     *result /= pow(((double *)argv)[1], 3);
jamie@52: 
jamie@56:     return XTRACT_SUCCESS;
jamie@52: }
jamie@52: 
jamie@146: int xtract_spectral_kurtosis(const double *data, const int N, const void *argv,  double *result)
jamie@140: {
jamie@52: 
jamie@53:     int m;
jamie@146:     const double *freqs, *amps;
jamie@52: 
jamie@53:     m = N >> 1;
jamie@53: 
jamie@53:     amps = data;
jamie@53:     freqs = data + m;
jamie@52: 
jamie@146:     *result = 0.0;
jamie@113: 
jamie@123:     while(m--)
jamie@146:         *result += pow(freqs[m] - ((double *)argv)[0], 4) * amps[m];
jamie@52: 
jamie@146:     *result /= pow(((double *)argv)[1], 4);
jamie@146:     *result -= 3.0;
jamie@52: 
jamie@56:     return XTRACT_SUCCESS;
jamie@52: }
jamie@52: 
jamie@146: int xtract_irregularity_k(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@25: 
jamie@1:     int n,
jamie@140:         M = N - 1;
jamie@140: 
jamie@146:     *result = 0.0;
jamie@140: 
jamie@1:     for(n = 1; n < M; n++)
jamie@146:         *result += fabs(data[n] - (data[n-1] + data[n] + data[n+1]) / 3.0);
jamie@1: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@146: int xtract_irregularity_j(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@25: 
sean@196:     int n = N - 1;
jamie@1: 
jamie@146:     double num = 0.0, den = 0.0;
jamie@1: 
jamie@140:     while(n--)
jamie@140:     {
jamie@146:         num += pow(data[n] - data[n+1], 2);
jamie@146:         den += pow(data[n], 2);
jamie@1:     }
jamie@25: 
jamie@146:     *result = (double)(num / den);
jamie@1: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@146: int xtract_tristimulus_1(const double *data, const int N, const void *argv, double *result)
jamie@140: {
csaba@280:     int n = N >> 1, i;
csaba@280:     double den = 0.0, p1 = 0.0, fund = 0.0, temp = 0.0, h = 0.0;
sean@227:     const double *freqs;
jamie@1: 
sean@227:     fund = *(double *)argv;
sean@227:     freqs = data + n;
jamie@1: 
sean@227:     for(i = 0; i < n; i++)
jamie@140:     {
sean@227:         if((temp = data[i]))
jamie@140:         {
jamie@140:             den += temp;
sean@245:             h = floor(freqs[i] / fund + 0.5);
csaba@280:             if(h > 0 && h < 2 && (int)h == 1)
sean@227:                 p1 += temp;
jamie@140:         }
jamie@42:     }
jamie@42: 
jamie@146:     if(den == 0.0 || p1 == 0.0)
jamie@140:     {
jamie@146:         *result = 0.0;
jamie@113:         return XTRACT_NO_RESULT;
jamie@113:     }
jamie@140:     else
jamie@140:     {
jamie@113:         *result = p1 / den;
jamie@113:         return XTRACT_SUCCESS;
jamie@113:     }
jamie@1: }
jamie@1: 
jamie@146: int xtract_tristimulus_2(const double *data, const int N, const void *argv, double *result)
jamie@140: {
csaba@280:     int n = N >> 1, i;
csaba@280:     double den, p2, p3, p4, ps, fund, temp, h;
csaba@280:     den = p2 = p3 = p4 = ps = fund = temp = h = 0.0;
sean@228:     const double *freqs;
jamie@1: 
sean@228:     fund = *(double *)argv;
sean@228:     freqs = data + n;
jamie@1: 
sean@228:     for(i = 0; i < n; i++)
jamie@140:     {
sean@228:         if((temp = data[i]))
jamie@140:         {
jamie@140:             den += temp;
sean@245:             h = floor(freqs[i] / fund + 0.5);
csaba@280:             if (h > 1 && h < 5)
sean@228:             {
csaba@280:                 switch ((int)h)
csaba@280:                 {
csaba@280:                     case 2:
csaba@280:                         p2 += temp;
sean@228:                     break;
sean@228: 
csaba@280:                     case 3:
csaba@280:                         p3 += temp;
sean@228:                     break;
sean@228: 
csaba@280:                     case 4:
csaba@280:                         p4 += temp;
sean@228: 
csaba@280:                     default:
csaba@280:                         break;
csaba@280:                 }
sean@228:             }
jamie@140:         }
jamie@42:     }
jamie@42: 
jamie@113:     ps = p2 + p3 + p4;
jamie@25: 
jamie@146:     if(den == 0.0 || ps == 0.0)
jamie@140:     {
jamie@146:         *result = 0.0;
jamie@113:         return XTRACT_NO_RESULT;
jamie@113:     }
jamie@140:     else
jamie@140:     {
jamie@113:         *result = ps / den;
jamie@113:         return XTRACT_SUCCESS;
jamie@113:     }
jamie@113: 
jamie@1: }
jamie@1: 
jamie@146: int xtract_tristimulus_3(const double *data, const int N, const void *argv, double *result)
jamie@140: {
csaba@280:     int n = N >> 1, i;
csaba@280:     double den = 0.0, num = 0.0, fund = 0.0, temp = 0.0, h = 0.0;
sean@229:     const double *freqs;
jamie@25: 
sean@229:     fund = *(double *)argv;
sean@229:     freqs = data + n;
jamie@1: 
sean@229:     for(i = 0; i < n; i++)
jamie@140:     {
sean@229:         if((temp = data[i]))
jamie@140:         {
jamie@140:             den += temp;
sean@245:             h = floor(freqs[i] / fund + 0.5);
sean@229:             if(h >= 5)
jamie@140:                 num += temp;
jamie@140:         }
jamie@42:     }
jamie@25: 
jamie@146:     if(den == 0.0 || num == 0.0)
jamie@140:     {
jamie@146:         *result = 0.0;
jamie@113:         return XTRACT_NO_RESULT;
jamie@113:     }
jamie@140:     else
jamie@140:     {
jamie@113:         *result = num / den;
jamie@113:         return XTRACT_SUCCESS;
jamie@113:     }
jamie@1: }
jamie@1: 
jamie@146: int xtract_smoothness(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@25: 
jamie@184:     int n; 
jamie@184:     int M = N - 1;
jamie@184:     double prev = 0.0;
jamie@184:     double current = 0.0;
jamie@184:     double next = 0.0;
jamie@184:     double temp = 0.0;
jamie@1: 
jamie@184:     
jamie@59: 
jamie@140:     for(n = 1; n < M; n++)
jamie@140:     {
jamie@184:         if(n == 1)
jamie@184:         {
jamie@184:             prev = data[n-1] <= 0 ? XTRACT_LOG_LIMIT : data[n-1];
jamie@184:             current = data[n] <= 0 ? XTRACT_LOG_LIMIT : data[n];
jamie@184:         }
jamie@184:         else
jamie@184:         {
jamie@184:             prev = current;
jamie@184:             current = next;
jamie@184:         }
jamie@184:         
jamie@184:         next = data[n+1] <= 0 ? XTRACT_LOG_LIMIT : data[n+1];
jamie@184:         
jamie@184:         temp += fabs(20.0 * log(current) - (20.0 * log(prev) +
jamie@184:                          20.0 * log(current) + 20.0 * log(next)) / 3.0);
jamie@25:     }
jamie@184:     
jamie@184:     *result = temp;
jamie@44: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@146: int xtract_spread(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@1: 
jamie@140:     return xtract_spectral_variance(data, N, argv, result);
jamie@1: }
jamie@1: 
jamie@146: int xtract_zcr(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@1: 
jamie@1:     int n = N;
jamie@25: 
jamie@1:     for(n = 1; n < N; n++)
jamie@140:         if(data[n] * data[n-1] < 0) (*result)++;
jamie@25: 
jamie@146:     *result /= (double)N;
jamie@25: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@146: int xtract_rolloff(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@1: 
jamie@1:     int n = N;
jamie@146:     double pivot, temp, percentile;
jamie@42: 
jamie@146:     pivot = temp = 0.0;
jamie@146:     percentile = ((double *)argv)[1];
jamie@25: 
jamie@140:     while(n--) pivot += data[n];
jamie@25: 
jamie@146:     pivot *= percentile / 100.0;
jamie@25: 
jamie@42:     for(n = 0; temp < pivot; n++)
jamie@140:         temp += data[n];
jamie@1: 
jamie@146:     *result = n * ((double *)argv)[0];
jamie@146:     /* *result = (n / (double)N) * (((double *)argv)[1] * .5); */
jamie@25: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@146: int xtract_loudness(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@25: 
jamie@47:     int n = N, rv;
jamie@25: 
jamie@146:     *result = 0.0;
jamie@113: 
jamie@140:     if(n > XTRACT_BARK_BANDS)
jamie@140:     {
jamie@140:         n = XTRACT_BARK_BANDS;
jamie@140:         rv = XTRACT_BAD_VECTOR_SIZE;
jamie@93:     }
jamie@47:     else
jamie@140:         rv = XTRACT_SUCCESS;
jamie@1: 
jamie@1:     while(n--)
jamie@146:         *result += pow(data[n], 0.23);
jamie@38: 
jamie@47:     return rv;
jamie@1: }
jamie@1: 
jamie@146: int xtract_flatness(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@1: 
jamie@113:     int n, count, denormal_found;
jamie@1: 
jamie@140:     double num, den, temp;
jamie@25: 
jamie@146:     num = 1.0;
jamie@146:     den = temp = 0.0;
jamie@43: 
jamie@113:     denormal_found = 0;
jamie@113:     count = 0;
jamie@113: 
jamie@140:     for(n = 0; n < N; n++)
jamie@140:     {
jamie@146:         if((temp = data[n]) != 0.0)
jamie@140:         {
jamie@140:             if (xtract_is_denormal(num))
jamie@140:             {
jamie@113:                 denormal_found = 1;
jamie@113:                 break;
jamie@113:             }
jamie@113:             num *= temp;
jamie@113:             den += temp;
jamie@113:             count++;
jamie@113:         }
jamie@1:     }
jamie@44: 
jamie@140:     if(!count)
jamie@140:     {
jamie@146:         *result = 0.0;
jamie@113:         return XTRACT_NO_RESULT;
jamie@113:     }
jamie@25: 
jamie@146:     num = pow(num, 1.0 / (double)N);
jamie@146:     den /= (double)N;
jamie@44: 
jamie@44: 
jamie@146:     *result = (double) (num / den);
jamie@113: 
jamie@113:     if(denormal_found)
jamie@113:         return XTRACT_DENORMAL_FOUND;
jamie@113:     else
jamie@113:         return XTRACT_SUCCESS;
jamie@140: 
jamie@113: }
jamie@113: 
jamie@146: int xtract_flatness_db(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@113: 
jamie@146:     double flatness;
jamie@113: 
jamie@146:     flatness = *(double *)argv;
jamie@113: 
jamie@140:     if (flatness <= 0)
jamie@115:         flatness = XTRACT_LOG_LIMIT;
jamie@113: 
jamie@182:     *result = 10 * log10(flatness);
jamie@25: 
jamie@56:     return XTRACT_SUCCESS;
jamie@44: 
jamie@1: }
jamie@1: 
jamie@146: int xtract_tonality(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@25: 
jamie@146:     double sfmdb;
jamie@25: 
jamie@146:     sfmdb = *(double *)argv;
jamie@1: 
jamie@146:     *result = XTRACT_MIN(sfmdb / -60.0, 1);
jamie@25: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@146: int xtract_crest(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@25: 
jamie@146:     double max, mean;
jamie@45: 
jamie@146:     max = mean = 0.0;
jamie@45: 
jamie@146:     max = *(double *)argv;
jamie@146:     mean = *((double *)argv+1);
jamie@45: 
jamie@45:     *result = max / mean;
jamie@45: 
jamie@56:     return XTRACT_SUCCESS;
jamie@25: 
jamie@1: }
jamie@1: 
jamie@146: int xtract_noisiness(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@25: 
jamie@146:     double h, i, p; /*harmonics, inharmonics, partials */
jamie@45: 
jamie@146:     i = p = h = 0.0;
jamie@45: 
jamie@146:     h = *(double *)argv;
jamie@146:     p = *((double *)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@146: int xtract_rms_amplitude(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@1: 
jamie@1:     int n = N;
jamie@1: 
jamie@146:     *result = 0.0;
jamie@113: 
jamie@56:     while(n--) *result += XTRACT_SQ(data[n]);
jamie@1: 
jamie@146:     *result = sqrt(*result / (double)N);
jamie@25: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@146: int xtract_spectral_inharmonicity(const double *data, const int N, const void *argv, double *result)
jamie@140: {
sean@223:     int n = N >> 1, h = 0;
jamie@146:     double num = 0.0, den = 0.0, fund;
jamie@146:     const double *freqs, *amps;
jamie@1: 
jamie@146:     fund = *(double *)argv;
jamie@52:     amps = data;
jamie@52:     freqs = data + n;
jamie@25: 
jamie@140:     while(n--)
jamie@140:     {
jamie@140:         if(amps[n])
jamie@140:         {
sean@245:             h = floor(freqs[n] / fund + 0.5);
sean@223:             num += fabs(freqs[n] - h * fund) * XTRACT_SQ(amps[n]);
jamie@140:             den += XTRACT_SQ(amps[n]);
jamie@140:         }
jamie@1:     }
jamie@1: 
jamie@140:     *result = (2 * num) / (fund * den);
jamie@25: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@1: 
jamie@146: int xtract_power(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@1: 
jamie@56:     return XTRACT_FEATURE_NOT_IMPLEMENTED;
jamie@25: 
jamie@1: }
jamie@1: 
jamie@146: int xtract_odd_even_ratio(const double *data, const int N, const void *argv, double *result)
jamie@140: {
sean@225:     int n = N >> 1, h = 0;
sean@225:     double odd = 0.0, even = 0.0, fund, temp;
sean@225:     const double *freqs;
jamie@1: 
sean@225:     fund = *(double *)argv;
sean@225:     freqs = data + n;
jamie@1: 
sean@225:     while(n--)
jamie@140:     {
jamie@140:         if((temp = data[n]))
jamie@140:         {
sean@245:             h = floor(freqs[n] / fund + 0.5);
sean@225:             if(XTRACT_IS_ODD(h))
jamie@140:             {
jamie@140:                 odd += temp;
jamie@140:             }
jamie@140:             else
jamie@140:             {
jamie@140:                 even += temp;
jamie@140:             }
jamie@140:         }
jamie@1:     }
jamie@1: 
jamie@146:     if(odd == 0.0 || even == 0.0)
jamie@140:     {
jamie@146:         *result = 0.0;
jamie@113:         return XTRACT_NO_RESULT;
jamie@113:     }
jamie@140:     else
jamie@140:     {
jamie@113:         *result = odd / even;
jamie@113:         return XTRACT_SUCCESS;
jamie@113:     }
jamie@1: }
jamie@1: 
jamie@146: int xtract_sharpness(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@1: 
jamie@48:     int n = N, rv;
jamie@146:     double sl, g; /* sl = specific loudness */
jamie@140:     double temp;
jamie@48: 
jamie@146:     sl = g = 0.0;
jamie@146:     temp = 0.0;
jamie@48: 
jamie@140:     if(n > XTRACT_BARK_BANDS)
jamie@140:         rv = XTRACT_BAD_VECTOR_SIZE;
jamie@48:     else
jamie@140:         rv = XTRACT_SUCCESS;
jamie@48: 
jamie@48: 
jamie@140:     while(n--)
jamie@140:     {
jamie@146:         sl = pow(data[n], 0.23);
jamie@146:         g = (n < 15 ? 1.0 : 0.066 * exp(0.171 * n));
jamie@140:         temp += n * g * sl;
jamie@48:     }
jamie@48: 
jamie@146:     temp = 0.11 * temp / (double)N;
jamie@146:     *result = (double)temp;
jamie@48: 
jamie@48:     return rv;
jamie@25: 
jamie@1: }
jamie@1: 
jamie@146: int xtract_spectral_slope(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@1: 
jamie@146:     const double *freqs, *amps;
jamie@146:     double f, a,
jamie@140:           F, A, FA, FXTRACT_SQ; /* sums of freqs, amps, freq * amps, freq squared */
jamie@140:     int n, M;
jamie@140: 
jamie@146:     F = A = FA = FXTRACT_SQ = 0.0;
jamie@48:     n = M = N >> 1;
jamie@48: 
jamie@52:     amps = data;
jamie@52:     freqs = data + n;
jamie@48: 
jamie@140:     while(n--)
jamie@140:     {
jamie@140:         f = freqs[n];
jamie@140:         a = amps[n];
jamie@140:         F += f;
jamie@140:         A += a;
jamie@140:         FA += f * a;
jamie@140:         FXTRACT_SQ += f * f;
jamie@48:     }
jamie@48: 
jamie@146:     *result = (1.0 / A) * (M * FA - F * A) / (M * FXTRACT_SQ - F * F);
jamie@48: 
jamie@56:     return XTRACT_SUCCESS;
jamie@25: 
jamie@1: }
jamie@1: 
jamie@146: int xtract_lowest_value(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@25: 
jamie@45:     int n = N;
jamie@45: 
jamie@192:     *result = DBL_MAX;
jamie@45: 
jamie@140:     while(n--)
jamie@140:     {
jamie@192:         if(data[n] > *(double *)argv)
jamie@140:             *result = XTRACT_MIN(*result, data[n]);
jamie@45:     }
jamie@45: 
jamie@192:     if (*result == DBL_MAX)
jamie@192:         return XTRACT_NO_RESULT;
jamie@192:         
jamie@56:     return XTRACT_SUCCESS;
jamie@45: }
jamie@45: 
jamie@146: int xtract_highest_value(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@45: 
jamie@1:     int n = N;
jamie@1: 
jamie@46:     *result = data[--n];
jamie@44: 
jamie@140:     while(n--)
jamie@140:         *result = XTRACT_MAX(*result, data[n]);
jamie@44: 
jamie@56:     return XTRACT_SUCCESS;
jamie@1: }
jamie@1: 
jamie@45: 
jamie@146: int xtract_sum(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@45: 
jamie@45:     int n = N;
jamie@45: 
jamie@146:     *result = 0.0;
jamie@113: 
jamie@45:     while(n--)
jamie@140:         *result += *data++;
jamie@45: 
jamie@56:     return XTRACT_SUCCESS;
jamie@45: 
jamie@45: }
jamie@45: 
jamie@146: int xtract_nonzero_count(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@59: 
jamie@59:     int n = N;
jamie@140: 
jamie@146:     *result = 0.0;
jamie@59: 
jamie@59:     while(n--)
jamie@140:         *result += (*data++ ? 1 : 0);
jamie@59: 
jamie@59:     return XTRACT_SUCCESS;
jamie@59: 
jamie@59: }
jamie@59: 
sean@198: int xtract_hps(const double *data, const int N, const void *argv, double *result)
jamie@140: {
sean@198:     int n, M, i, peak_index, position1_lwr;
sean@198:     double tempProduct, peak, largest1_lwr, ratio1;
jamie@1: 
sean@198:     n = N / 2;
jamie@1: 
sean@198:     M = ceil(n / 3.0);
jamie@25: 
sean@198:     if (M <= 1)
jamie@140:     {
sean@198:         /* Input data is too short. */
sean@198:         *result = 0;
sean@198:         return XTRACT_NO_RESULT;
jamie@1:     }
jamie@25: 
sean@245:     peak_index = 0;
sean@245: 
sean@198:     tempProduct = peak = 0;
sean@198:     for (i = 0; i < M; ++i)
sean@198:     {
sean@198:         tempProduct = data [i] * data [i * 2] * data [i * 3];
jamie@25: 
sean@198:         if (tempProduct > peak)
jamie@140:         {
sean@198:             peak = tempProduct;
sean@198:             peak_index = i;
jamie@140:         }
jamie@1:     }
jamie@1: 
jamie@1:     largest1_lwr = position1_lwr = 0;
jamie@1: 
sean@198:     for(i = 0; i < N; ++i)
jamie@140:     {
sean@198:         if(data[i] > largest1_lwr && i != peak_index)
jamie@140:         {
sean@198:             largest1_lwr = data[i];
sean@198:             position1_lwr = i;
jamie@140:         }
jamie@1:     }
jamie@1: 
jamie@1:     ratio1 = data[position1_lwr] / data[peak_index];
jamie@1: 
jamie@140:     if(position1_lwr > peak_index * 0.4 && position1_lwr <
jamie@140:             peak_index * 0.6 && ratio1 > 0.1)
jamie@140:         peak_index = position1_lwr;
jamie@1: 
sean@198:     *result = data [n + peak_index];
sean@196: 
sean@196:     return XTRACT_SUCCESS;
jamie@1: }
jamie@5: 
jamie@146: int xtract_f0(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@5: 
jamie@78:     int M, tau, n;
jamie@146:     double sr;
jamie@43:     size_t bytes;
jamie@146:     double f0, err_tau_1, err_tau_x, array_max,
jamie@140:           threshold_peak, threshold_centre,
jamie@140:           *input;
jamie@22: 
jamie@146:     sr = *(double *)argv;
jamie@78:     if(sr == 0)
jamie@146:         sr = 44100.0;
jamie@43: 
andrea@211:     input = (double*)malloc(bytes = N * sizeof(double));
andrea@211:     input = (double*)memcpy(input, data, bytes);
jamie@146:     /*  threshold_peak = *((double *)argv+1);
jamie@146:     threshold_centre = *((double *)argv+2);
jamie@146:     printf("peak: %.2\tcentre: %.2\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@140:     for(n = 0; n < N; n++)
jamie@140:     {
jamie@140:         if (input[n] > array_max)
jamie@140:             array_max = input[n];
jamie@12:     }
jamie@25: 
jamie@25:     threshold_peak *= array_max;
jamie@25: 
jamie@25:     /* peak clip */
jamie@140:     for(n = 0; n < N; n++)
jamie@140:     {
jamie@140:         if(input[n] > threshold_peak)
jamie@140:             input[n] = threshold_peak;
jamie@140:         else if(input[n] < -threshold_peak)
jamie@140:             input[n] = -threshold_peak;
jamie@25:     }
jamie@25: 
jamie@25:     threshold_centre *= array_max;
jamie@25: 
jamie@25:     /* Centre clip */
jamie@140:     for(n = 0; n < N; n++)
jamie@140:     {
jamie@140:         if (input[n] < threshold_centre)
jamie@140:             input[n] = 0;
jamie@140:         else
jamie@140:             input[n] -= threshold_centre;
jamie@25:     }
jamie@25: 
jamie@25:     /* Estimate fundamental freq */
jamie@25:     for (n = 1; n < M; n++)
jamie@146:         err_tau_1 = err_tau_1 + fabs(input[n] - input[n+1]);
jamie@140:     /* 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@140:     for (tau = 2; tau < M; tau++)
jamie@140:     {
jamie@140:         err_tau_x = 0;
jamie@140:         for (n = 1; n < M; n++)
jamie@140:         {
jamie@146:             err_tau_x = err_tau_x + fabs(input[n] - input[n+tau]);
jamie@140:         }
jamie@140:         if (err_tau_x < err_tau_1)
jamie@140:         {
jamie@140:             f0 = sr / (tau + (err_tau_x / err_tau_1));
jamie@140:             *result = f0;
jamie@140:             free(input);
jamie@140:             return XTRACT_SUCCESS;
jamie@140:         }
jamie@25:     }
jamie@43:     *result = -0;
jamie@43:     free(input);
jamie@56:     return XTRACT_NO_RESULT;
jamie@5: }
jamie@43: 
jamie@146: int xtract_failsafe_f0(const double *data, const int N, const void *argv, double *result)
jamie@140: {
jamie@44: 
jamie@146:     double *spectrum = NULL, argf[2], *peaks = NULL, return_code, sr;
jamie@44: 
jamie@43:     return_code = xtract_f0(data, N, argv, result);
jamie@44: 
jamie@140:     if(return_code == XTRACT_NO_RESULT)
jamie@140:     {
jamie@146:         sr = *(double *)argv;
jamie@140:         if(sr == 0)
jamie@146:             sr = 44100.0;
jamie@146:         spectrum = (double *)malloc(N * sizeof(double));
jamie@146:         peaks = (double *)malloc(N * sizeof(double));
jamie@140:         argf[0] = sr;
jamie@140:         argf[1] = XTRACT_MAGNITUDE_SPECTRUM;
jamie@140:         xtract_spectrum(data, N, argf, spectrum);
jamie@146:         argf[1] = 10.0;
jamie@140:         xtract_peak_spectrum(spectrum, N >> 1, argf, peaks);
jamie@146:         argf[0] = 0.0;
jamie@140:         xtract_lowest_value(peaks+(N >> 1), N >> 1, argf, result);
jamie@44: 
jamie@140:         free(spectrum);
jamie@140:         free(peaks);
jamie@43:     }
jamie@43: 
jamie@56:     return XTRACT_SUCCESS;
jamie@43: 
jamie@43: }
jamie@44: 
jamie@161: int xtract_wavelet_f0(const double *data, const int N, const void *argv, double *result)
jamie@161: {
jamie@169:     /* double sr = *(double *)argv; */
jamie@161: 
jamie@161:     *result = dywapitch_computepitch(&wavelet_f0_state, data, 0, N);
jamie@161: 
jamie@161:     if (*result == 0.0)
jamie@161:     {
jamie@161:         return XTRACT_NO_RESULT;
jamie@161:     }
jamie@161: 
jamie@161:     return XTRACT_SUCCESS;
jamie@161: }
jamie@161: 
jamie@205: int xtract_midicent(const double *data, const int N, const void *argv, double *result)
jamie@205: {
jamie@205:     double f0 = *(double *)argv;
jamie@205:     double note = 0.0;
jamie@208:       
jamie@205:     note = 69 + log(f0 / 440.f) * 17.31234;
jamie@205:     note *= 100;
andrea@211:     note = floor( 0.5f + note ); // replace -> round(note);
jamie@161: 
jamie@209:     *result = note;
jamie@205:     
jamie@206:     if (note > 12700 || note < 0)
jamie@206:     {
jamie@206:         return XTRACT_ARGUMENT_ERROR;
jamie@206:     }
jamie@206:     
jamie@205:     return XTRACT_SUCCESS;
jamie@205: }
sean@198: 
jamie@219: int xtract_peak(const double *data, const int N, const void *argv, double *result)
jamie@215: {
jamie@215:     double threshold = *(double *)argv;
jamie@215:     double current = data[N - 1];
jamie@215:     double average = 0.0;
jamie@215:     double maximum = -DBL_MAX;
jamie@215:     
jamie@215:     for (uint32_t n = 0; n < N; ++n)
jamie@215:     {
jamie@215:         average += data[n];
jamie@215:         if (data[n] > maximum)
jamie@215:         {
jamie@215:             maximum = data[n];
jamie@215:         }
jamie@215:     }
jamie@215:     
jamie@215:     average /= (double)N;
jamie@241:         
jamie@215:     if (current != maximum)
jamie@215:     {
jamie@215:         return XTRACT_NO_RESULT;
jamie@215:     }
jamie@215:     
jamie@215:     if (current < average + threshold)
jamie@215:     {
jamie@215:         return XTRACT_NO_RESULT;
jamie@215:     }
jamie@215:     
jamie@241:     *result = current;
jamie@241:     
jamie@215:     return XTRACT_SUCCESS;
jamie@215:     
jamie@215: }
jamie@215: 
jamie@215: