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 jamie@43: #include jamie@1: jamie@43: int xtract_mean(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@42: *result += data[n]; jamie@25: jamie@1: *result /= N; jamie@38: jamie@38: return 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@38: return 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@1: *result = sqrt(*(float *)argv); jamie@25: jamie@38: return 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@38: return 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@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, 3); jamie@42: } jamie@1: jamie@42: *result /= N; jamie@44: jamie@38: return 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@38: return SUCCESS; jamie@1: } jamie@1: jamie@11: jamie@43: int xtract_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@25: freqs = data; jamie@38: amps = 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@38: return SUCCESS; jamie@11: } jamie@11: 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; jamie@1: 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@38: return 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@37: float 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@1: *result = num / den; jamie@1: jamie@38: return 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@38: return 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@38: return 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@38: return SUCCESS; jamie@1: } jamie@1: jamie@43: int xtract_smoothness(const float *data, const int N, const void *argv, float *result){ jamie@25: jamie@1: int n = N; jamie@1: jamie@43: float *input; jamie@43: jamie@43: input = (float *)malloc(N * sizeof(float)); jamie@43: input = memcpy(input, data, N * sizeof(float)); jamie@43: jamie@43: if (input[0] <= 0) input[0] = 1; jamie@43: if (input[1] <= 0) input[1] = 1; jamie@25: jamie@1: for(n = 2; n < N; n++){ jamie@43: if(input[n] <= 0) input[n] = 1; jamie@43: *result += abs(20 * log(input[n-1]) - (20 * log(input[n-2]) + jamie@43: 20 * log(input[n-1]) + 20 * log(input[n])) / 3); jamie@25: } jamie@43: jamie@43: free(input); jamie@44: jamie@38: return 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@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@38: return 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@38: return 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@42: float pivot, temp; jamie@42: jamie@42: pivot = temp = 0.f; jamie@25: jamie@1: while(n--) pivot += data[n]; jamie@25: jamie@42: pivot *= ((float *)argv)[0]; jamie@25: jamie@42: for(n = 0; temp < pivot; n++) jamie@42: temp += data[n]; jamie@1: jamie@42: *result = (n / (float)N) * (((float *)argv)[1] * .5); jamie@25: jamie@38: return 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@47: if(n > BARK_BANDS) jamie@47: rv = BAD_VECTOR_SIZE; jamie@47: else jamie@47: rv = SUCCESS; jamie@1: jamie@1: while(n--) jamie@25: *result += pow(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@44: num = pow(num, 1.f / N); jamie@1: den /= N; jamie@25: jamie@45: if(num < VERY_SMALL_NUMBER) jamie@45: num = VERY_SMALL_NUMBER; jamie@44: jamie@45: if(den < VERY_SMALL_NUMBER) jamie@45: den = VERY_SMALL_NUMBER; jamie@44: jamie@42: *result = num / den; jamie@25: jamie@38: return 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@45: sfmdb = (sfm > 0 ? ((10 * log10(sfm)) / -60) : 0); jamie@25: jamie@1: *result = MIN(sfmdb, 1); jamie@25: jamie@38: return 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@45: return 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@45: return 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@1: while(n--) *result += SQ(data[n]); jamie@1: jamie@1: *result = sqrt(*result / N); jamie@25: jamie@38: return SUCCESS; jamie@1: } jamie@1: jamie@43: int xtract_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@41: freqs = data; jamie@41: amps = data + n; jamie@25: jamie@1: while(n--){ jamie@41: num += abs(freqs[n] - n * fund) * SQ(amps[n]); jamie@41: den += SQ(amps[n]); jamie@1: } jamie@1: jamie@41: *result = (2 * num) / (fund * den); jamie@25: jamie@38: return 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@38: return 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@43: float num = 0.f, den = 0.f, temp, f0; jamie@1: jamie@43: f0 = *(float *)argv; jamie@44: jamie@43: for(n = 0; n < M; n++){ jamie@43: if((temp = data[n])){ jamie@43: if(((int)(rintf(temp / f0)) % 2) != 0){ jamie@43: num += data[M + n]; jamie@43: } jamie@43: else{ jamie@43: den += data[M + n]; jamie@43: } jamie@43: } jamie@1: } jamie@1: jamie@1: *result = num / den; jamie@25: jamie@38: return SUCCESS; jamie@1: } jamie@1: jamie@43: int xtract_sharpness(const float *data, const int N, const void *argv, float *result){ jamie@1: jamie@38: return FEATURE_NOT_IMPLEMENTED; jamie@25: jamie@1: } jamie@1: jamie@43: int xtract_slope(const float *data, const int N, const void *argv, float *result){ jamie@1: jamie@38: return FEATURE_NOT_IMPLEMENTED; 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@45: *result = MIN(*result, data[n]); jamie@45: } jamie@45: jamie@45: return 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@45: *result = MAX(*result, data[n]); jamie@44: jamie@38: return 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@45: return SUCCESS; jamie@45: jamie@45: } jamie@45: 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@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@38: return 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@25: int M, sr, tau, n; 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@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@25: return SUCCESS; jamie@25: } jamie@25: } jamie@43: *result = -0; jamie@43: free(input); jamie@25: return 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@43: float *magnitudes = NULL, argf[2], *peaks = NULL, return_code; jamie@44: jamie@43: return_code = xtract_f0(data, N, argv, result); jamie@44: jamie@43: if(return_code == NO_RESULT){ jamie@44: jamie@43: magnitudes = (float *)malloc(N * sizeof(float)); jamie@43: peaks = (float *)malloc(N * sizeof(float)); jamie@46: xtract_magnitude_spectrum(data, N, argv, magnitudes); jamie@43: argf[0] = 10.f; jamie@43: argf[1] = *(float *)argv; jamie@43: xtract_peaks(magnitudes, N, argf, peaks); jamie@43: argf[0] = 0.f; jamie@45: xtract_lowest_value(peaks, N >> 1, argf, result); jamie@44: jamie@43: free(magnitudes); jamie@43: free(peaks); jamie@43: } jamie@43: jamie@43: return SUCCESS; jamie@43: jamie@43: } jamie@44: