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@1: /* xtract_vector.c: defines functions that extract a feature as a single value from an input vector */ jamie@1: jamie@1: #include jamie@43: #include jamie@43: #include jamie@171: #include jamie@30: jamie@148: #include "fft.h" jamie@140: jamie@259: #include "xtract/libxtract.h" jamie@98: #include "xtract_macros_private.h" jamie@146: #include "xtract_globals_private.h" jamie@98: jamie@154: #ifndef M_PI jamie@154: #define M_PI 3.14159265358979323846264338327 jamie@154: #endif jamie@154: jamie@146: int xtract_spectrum(const double *data, const int N, const void *argv, double *result) jamie@140: { jamie@1: jamie@140: int vector = 0; jamie@140: int withDC = 0; jamie@140: int normalise = 0; jamie@146: double q = 0.0; jamie@146: double temp = 0.0; jamie@146: double max = 0.0; jamie@146: double NxN = XTRACT_SQ(N); jamie@146: double *marker = NULL; jamie@150: double real = 0.0; jamie@150: double imag = 0.0; jamie@140: unsigned int n = 0; jamie@140: unsigned int m = 0; jamie@140: unsigned int M = N >> 1; jamie@171: #ifdef USE_OOURA jamie@171: double *fft = NULL; jamie@171: #else jamie@150: DSPDoubleSplitComplex *fft = NULL; jamie@150: #endif jamie@98: jamie@146: q = *(double *)argv; jamie@146: vector = (int)*((double *)argv+1); jamie@146: withDC = (int)*((double *)argv+2); jamie@146: normalise = (int)*((double *)argv+3); jamie@105: jamie@56: XTRACT_CHECK_q; jamie@150: #ifdef USE_OOURA jamie@140: if(!ooura_data_spectrum.initialised) jamie@150: #else jamie@150: if(!vdsp_data_spectrum.initialised) jamie@150: #endif jamie@140: { jamie@140: fprintf(stderr, jamie@140: "libxtract: error: xtract_spectrum() failed, " jamie@140: "fft data unitialised.\n"); jamie@98: return XTRACT_NO_RESULT; jamie@98: } jamie@98: jamie@150: #ifdef USE_OOURA jamie@140: /* ooura is in-place jamie@147: * the output format is jamie@140: * a[0] - DC, a[1] - nyquist, a[2...N-1] - remaining bins jamie@140: */ andrea@211: fft = (double*)malloc(N * sizeof(double)); jamie@171: assert(fft != NULL); jamie@171: memcpy(fft, data, N * sizeof(double)); jamie@171: jamie@171: rdft(N, 1, fft, ooura_data_spectrum.ooura_ip, jamie@140: ooura_data_spectrum.ooura_w); jamie@150: #else jamie@150: fft = &vdsp_data_spectrum.fft; jamie@150: vDSP_ctozD((DSPDoubleComplex *)data, 2, fft, 1, N >> 1); jamie@150: vDSP_fft_zripD(vdsp_data_spectrum.setup, fft, 1, jamie@150: vdsp_data_spectrum.log2N, FFT_FORWARD); jamie@150: #endif jamie@54: jamie@140: switch(vector) jamie@140: { jamie@67: jamie@120: case XTRACT_LOG_MAGNITUDE_SPECTRUM: jamie@140: for(n = 0, m = 0; m < M; ++n, ++m) jamie@140: { jamie@150: marker = &result[m]; jamie@150: jamie@150: if(n==0 && !withDC) /* discard DC and keep Nyquist */ jamie@140: { jamie@150: ++n; jamie@150: #ifdef USE_OOURA jamie@150: marker = &result[M-1]; jamie@150: #endif jamie@140: } jamie@150: #ifdef USE_OOURA jamie@150: if(n==1 && withDC) /* discard Nyquist */ jamie@150: { jamie@150: ++n; jamie@150: } jamie@178: if(n == M) jamie@178: { jamie@178: XTRACT_SET_FREQUENCY; jamie@178: break; jamie@178: } jamie@150: jamie@171: real = fft[n*2]; jamie@171: imag = fft[n*2+1]; jamie@150: #else jamie@150: real = fft->realp[n]; jamie@150: imag = fft->realp[n]; jamie@150: #endif jamie@150: jamie@150: temp = XTRACT_SQ(real) + XTRACT_SQ(imag); jamie@140: if (temp > XTRACT_LOG_LIMIT) jamie@140: { jamie@146: temp = log(sqrt(temp) / (double)N); jamie@140: } jamie@140: else jamie@140: { jamie@140: temp = XTRACT_LOG_LIMIT_DB; jamie@140: } jamie@140: result[m] = jamie@140: /* Scaling */ jamie@140: (temp + XTRACT_DB_SCALE_OFFSET) / jamie@140: XTRACT_DB_SCALE_OFFSET; jamie@111: jamie@140: XTRACT_SET_FREQUENCY; jamie@140: XTRACT_GET_MAX; jamie@140: } jamie@140: break; jamie@67: jamie@140: case XTRACT_POWER_SPECTRUM: jamie@140: for(n = 0, m = 0; m < M; ++n, ++m) jamie@140: { jamie@150: jamie@150: marker = &result[m]; jamie@150: jamie@150: if(n==0 && !withDC) /* discard DC and keep Nyquist */ jamie@150: { jamie@150: ++n; jamie@150: #ifdef USE_OOURA jamie@150: marker = &result[M-1]; jamie@150: #endif jamie@150: } jamie@150: #ifdef USE_OOURA jamie@150: if(n==1 && withDC) /* discard Nyquist */ jamie@140: { jamie@140: ++n; jamie@120: } jamie@178: if(n == M) jamie@178: { jamie@178: XTRACT_SET_FREQUENCY; jamie@178: break; jamie@178: } jamie@150: jamie@171: real = fft[n*2]; jamie@171: imag = fft[n*2+1]; jamie@150: #else jamie@150: real = fft->realp[n]; jamie@150: imag = fft->realp[n]; jamie@150: #endif jamie@150: jamie@150: result[m] = (XTRACT_SQ(real) + XTRACT_SQ(imag)) / NxN; jamie@140: XTRACT_SET_FREQUENCY; jamie@140: XTRACT_GET_MAX; jamie@140: } jamie@140: break; jamie@120: jamie@140: case XTRACT_LOG_POWER_SPECTRUM: jamie@140: for(n = 0, m = 0; m < M; ++n, ++m) jamie@140: { jamie@150: marker = &result[m]; jamie@150: jamie@150: if(n==0 && !withDC) /* discard DC and keep Nyquist */ jamie@150: { jamie@150: ++n; jamie@150: #ifdef USE_OOURA jamie@150: marker = &result[M-1]; jamie@150: #endif jamie@150: } jamie@150: #ifdef USE_OOURA jamie@150: if(n==1 && withDC) /* discard Nyquist */ jamie@140: { jamie@140: ++n; jamie@120: } jamie@178: if(n == M) jamie@178: { jamie@178: XTRACT_SET_FREQUENCY; jamie@178: break; jamie@178: } jamie@150: jamie@171: real = fft[n*2]; jamie@171: imag = fft[n*2+1]; jamie@150: #else jamie@150: real = fft->realp[n]; jamie@150: imag = fft->realp[n]; jamie@150: #endif jamie@150: jamie@150: if ((temp = XTRACT_SQ(real) + XTRACT_SQ(imag)) > jamie@140: XTRACT_LOG_LIMIT) jamie@146: temp = log(temp / NxN); jamie@140: else jamie@140: temp = XTRACT_LOG_LIMIT_DB; jamie@67: jamie@140: result[m] = (temp + XTRACT_DB_SCALE_OFFSET) / jamie@140: XTRACT_DB_SCALE_OFFSET; jamie@140: XTRACT_SET_FREQUENCY; jamie@140: XTRACT_GET_MAX; jamie@140: } jamie@140: break; jamie@111: jamie@140: default: jamie@140: /* MAGNITUDE_SPECTRUM */ jamie@140: for(n = 0, m = 0; m < M; ++n, ++m) jamie@140: { jamie@140: marker = &result[m]; jamie@140: jamie@140: if(n==0 && !withDC) /* discard DC and keep Nyquist */ jamie@140: { jamie@140: ++n; jamie@150: #ifdef USE_OOURA jamie@140: marker = &result[M-1]; jamie@150: #endif jamie@120: } jamie@150: #ifdef USE_OOURA jamie@140: if(n==1 && withDC) /* discard Nyquist */ jamie@140: { jamie@140: ++n; jamie@140: } jamie@178: if(n == M) jamie@178: { jamie@178: XTRACT_SET_FREQUENCY; jamie@178: break; jamie@178: } jamie@178: jamie@171: real = fft[n*2]; jamie@171: imag = fft[n*2+1]; jamie@150: #else jamie@150: real = fft->realp[n]; jamie@150: imag = fft->realp[n]; jamie@150: #endif jamie@150: *marker = sqrt(XTRACT_SQ(real) + XTRACT_SQ(imag)) / (double)N; jamie@111: jamie@140: XTRACT_SET_FREQUENCY; jamie@140: XTRACT_GET_MAX; jamie@140: } jamie@140: break; jamie@70: } jamie@105: jamie@140: if(normalise) jamie@140: { jamie@105: for(n = 0; n < M; n++) jamie@105: result[n] /= max; jamie@105: } jamie@105: jamie@171: #ifdef USE_OOURA jamie@171: free(fft); jamie@171: #endif jamie@171: jamie@56: return XTRACT_SUCCESS; jamie@1: } jamie@1: jamie@146: int xtract_autocorrelation_fft(const double *data, const int N, const void *argv, double *result) jamie@140: { jamie@120: jamie@140: double *rfft = NULL; jamie@140: int n = 0; jamie@140: int M = 0; jamie@150: #ifndef USE_OOURA jamie@150: DSPDoubleSplitComplex *fft = NULL; jamie@154: double M_double = 0.0; jamie@150: #endif jamie@1: jamie@75: M = N << 1; jamie@43: jamie@75: /* Zero pad the input vector */ jamie@140: rfft = (double *)calloc(M, sizeof(double)); jamie@146: memcpy(rfft, data, N * sizeof(double)); jamie@75: jamie@150: #ifdef USE_OOURA jamie@140: rdft(M, 1, rfft, ooura_data_autocorrelation_fft.ooura_ip, jamie@140: ooura_data_autocorrelation_fft.ooura_w); jamie@1: jamie@150: for(n = 2; n < M; ++n) jamie@140: { jamie@140: rfft[n*2] = XTRACT_SQ(rfft[n*2]) + XTRACT_SQ(rfft[n*2+1]); jamie@146: rfft[n*2+1] = 0.0; jamie@75: } jamie@120: jamie@140: rfft[0] = XTRACT_SQ(rfft[0]); jamie@140: rfft[1] = XTRACT_SQ(rfft[1]); jamie@75: jamie@140: rdft(M, -1, rfft, ooura_data_autocorrelation_fft.ooura_ip, jamie@140: ooura_data_autocorrelation_fft.ooura_w); jamie@120: jamie@150: #else jamie@150: /* vDSP has its own autocorrelation function, but it doesn't fit the jamie@150: * LibXtract model, e.g. we can't guarantee it's going to use jamie@150: * an FFT for all values of N */ jamie@150: fft = &vdsp_data_autocorrelation_fft.fft; jamie@150: vDSP_ctozD((DSPDoubleComplex *)data, 2, fft, 1, N); jamie@150: vDSP_fft_zripD(vdsp_data_autocorrelation_fft.setup, fft, 1, jamie@150: vdsp_data_autocorrelation_fft.log2N, FFT_FORWARD); jamie@150: jamie@150: for(n = 0; n < N; ++n) jamie@150: { jamie@150: fft->realp[n] = XTRACT_SQ(fft->realp[n]) + XTRACT_SQ(fft->imagp[n]); jamie@150: fft->imagp[n] = 0.0; jamie@150: } jamie@150: jamie@150: vDSP_fft_zripD(vdsp_data_autocorrelation_fft.setup, fft, 1, jamie@150: vdsp_data_autocorrelation_fft.log2N, FFT_INVERSE); jamie@150: #endif jamie@150: jamie@75: /* Normalisation factor */ jamie@75: M = M * N; jamie@75: jamie@150: #ifdef USE_OOURA jamie@75: for(n = 0; n < N; n++) jamie@146: result[n] = rfft[n] / (double)M; jamie@140: free(rfft); jamie@150: #else jamie@150: M_double = (double)M; jamie@150: vDSP_ztocD(fft, 1, (DOUBLE_COMPLEX *)result, 2, N); jamie@150: vDSP_vsdivD(result, 1, &M_double, result, 1, N); jamie@150: #endif jamie@38: jamie@56: return XTRACT_SUCCESS; jamie@1: } jamie@1: jamie@146: int xtract_mfcc(const double *data, const int N, const void *argv, double *result) jamie@140: { jamie@30: jamie@30: xtract_mel_filter *f; jamie@30: int n, filter; jamie@30: jamie@30: f = (xtract_mel_filter *)argv; jamie@120: jamie@140: for(filter = 0; filter < f->n_filters; filter++) jamie@140: { jamie@146: result[filter] = 0.0; jamie@140: for(n = 0; n < N; n++) jamie@140: { jamie@71: result[filter] += data[n] * f->filters[filter][n]; jamie@30: } jamie@146: result[filter] = log(result[filter] < XTRACT_LOG_LIMIT ? XTRACT_LOG_LIMIT : result[filter]); jamie@30: } jamie@30: jamie@30: xtract_dct(result, f->n_filters, NULL, result); jamie@120: jamie@56: return XTRACT_SUCCESS; jamie@30: } jamie@30: jamie@146: int xtract_dct(const double *data, const int N, const void *argv, double *result) jamie@140: { jamie@120: jamie@148: int n; jamie@148: int m; andrea@211: double *temp = (double*)calloc(N, sizeof(double)); jamie@120: jamie@148: for (n = 0; n < N; ++n) jamie@148: { jamie@148: for(m = 1; m <= N; ++m) { jamie@148: temp[n] += data[m - 1] * cos(M_PI * (n / (double)N) * (m - 0.5)); jamie@140: } jamie@148: } jamie@120: jamie@148: memcpy(result, temp, N * sizeof(double)); jamie@148: free(temp); jamie@148: jamie@148: return XTRACT_SUCCESS; jamie@30: } jamie@30: jamie@146: int xtract_autocorrelation(const double *data, const int N, const void *argv, double *result) jamie@140: { jamie@30: jamie@30: /* Naive time domain implementation */ jamie@120: jamie@30: int n = N, i; jamie@120: jamie@146: double corr; jamie@30: jamie@140: while(n--) jamie@140: { jamie@120: corr = 0; jamie@140: for(i = 0; i < N - n; i++) jamie@140: { 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@146: int xtract_amdf(const double *data, const int N, const void *argv, double *result) jamie@140: { jamie@1: jamie@1: int n = N, i; jamie@120: jamie@146: double md, temp; jamie@1: jamie@140: while(n--) jamie@140: { jamie@146: md = 0.0; jamie@140: for(i = 0; i < N - n; i++) jamie@140: { jamie@6: temp = data[i] - data[i + n]; jamie@120: temp = (temp < 0 ? -temp : temp); jamie@120: md += temp; jamie@1: } jamie@146: result[n] = md / (double)N; jamie@1: } jamie@38: jamie@56: return XTRACT_SUCCESS; jamie@1: } jamie@1: jamie@146: int xtract_asdf(const double *data, const int N, const void *argv, double *result) jamie@140: { jamie@120: jamie@1: int n = N, i; jamie@120: jamie@146: double sd; jamie@1: jamie@140: while(n--) jamie@140: { jamie@146: sd = 0.0; jamie@140: for(i = 0; i < N - n; i++) jamie@140: { jamie@6: /*sd = 1;*/ jamie@56: sd += XTRACT_SQ(data[i] - data[i + n]); jamie@1: } jamie@146: result[n] = sd / (double)N; jamie@1: } jamie@38: jamie@56: return XTRACT_SUCCESS; jamie@1: } jamie@1: jamie@146: int xtract_bark_coefficients(const double *data, const int N, const void *argv, double *result) jamie@140: { jamie@1: jamie@1: int *limits, band, n; jamie@1: jamie@1: limits = (int *)argv; jamie@120: jamie@140: for(band = 0; band < XTRACT_BARK_BANDS - 1; band++) jamie@140: { jamie@146: result[band] = 0.0; 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@146: int xtract_peak_spectrum(const double *data, const int N, const void *argv, double *result) jamie@140: { jamie@1: jamie@146: double 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@146: threshold = max = y = y2 = y3 = p = q = 0.0; jamie@120: jamie@140: if(argv != NULL) jamie@140: { jamie@146: q = ((double *)argv)[0]; jamie@146: threshold = ((double *)argv)[1]; jamie@1: } jamie@49: else jamie@56: rv = XTRACT_BAD_ARGV; jamie@49: jamie@140: if(threshold < 0 || threshold > 100) jamie@140: { jamie@55: threshold = 0; jamie@56: rv = XTRACT_BAD_ARGV; jamie@1: } jamie@1: jamie@56: XTRACT_CHECK_q; jamie@49: jamie@146: input = (double *)calloc(N, sizeof(double)); jamie@98: jamie@146: bytes = N * sizeof(double); jamie@43: jamie@43: if(input != NULL) andrea@211: input = (double*)memcpy(input, data, bytes); jamie@43: else jamie@120: return XTRACT_MALLOC_FAILED; jamie@43: jamie@45: while(n--) jamie@56: max = XTRACT_MAX(max, input[n]); jamie@120: jamie@55: threshold *= .01 * max; jamie@1: jamie@1: result[0] = 0; jamie@59: result[N] = 0; jamie@1: jamie@140: for(n = 1; n < N; n++) jamie@140: { jamie@140: if(input[n] >= threshold) jamie@140: { jamie@140: if(input[n] > input[n - 1] && n + 1 < N && input[n] > input[n + 1]) jamie@140: { jamie@180: result[N + n] = q * (n + 1 + (p = .5 * ((y = input[n-1]) - jamie@140: (y3 = input[n+1])) / (input[n - 1] - 2 * jamie@140: (y2 = input[n]) + input[n + 1]))); jamie@52: result[n] = y2 - .25 * (y - y3) * p; jamie@1: } jamie@140: else jamie@140: { jamie@1: result[n] = 0; jamie@59: result[N + n] = 0; jamie@1: } jamie@1: } jamie@140: else jamie@140: { jamie@1: result[n] = 0; jamie@59: result[N + n] = 0; jamie@1: } jamie@140: } jamie@120: jamie@179: free(input); jamie@56: return (rv ? rv : XTRACT_SUCCESS); jamie@1: } jamie@120: jamie@146: int xtract_harmonic_spectrum(const double *data, const int N, const void *argv, double *result) jamie@140: { jamie@120: jamie@140: int n = (N >> 1), M = n; jamie@38: jamie@146: const double *freqs, *amps; jamie@146: double f0, threshold, ratio, nearest, distance; jamie@38: jamie@52: amps = data; jamie@52: freqs = data + n; jamie@146: f0 = *((double *)argv); jamie@146: threshold = *((double *)argv+1); jamie@38: jamie@146: ratio = nearest = distance = 0.0; jamie@38: jamie@140: while(n--) jamie@140: { jamie@140: if(freqs[n]) jamie@140: { jamie@120: ratio = freqs[n] / f0; andrea@211: nearest = floor( 0.5f + ratio); // replace -> nearest = round(ratio); andrea@211: distance = fabs(nearest - ratio); jamie@120: if(distance > threshold) jamie@146: result[n] = result[M + n] = 0.0; jamie@140: else jamie@140: { jamie@120: result[n] = amps[n]; jamie@120: result[M + n] = freqs[n]; jamie@120: } jamie@120: } jamie@120: else jamie@146: result[n] = result[M + n] = 0.0; jamie@38: } jamie@56: return XTRACT_SUCCESS; jamie@38: } jamie@120: jamie@146: int xtract_lpc(const double *data, const int N, const void *argv, double *result) jamie@140: { jamie@104: jamie@154: int i, j, M, L; jamie@146: double r = 0.0, jamie@146: error = 0.0; jamie@104: jamie@146: double *ref = NULL, jamie@140: *lpc = NULL ; jamie@104: jamie@104: error = data[0]; jamie@104: L = N - 1; /* The number of LPC coefficients */ jamie@104: M = L * 2; /* The length of *result */ jamie@104: ref = result; jamie@104: lpc = result+L; jamie@113: jamie@140: if(error == 0.0) jamie@140: { jamie@146: memset(result, 0, M * sizeof(double)); jamie@104: return XTRACT_NO_RESULT; jamie@104: } jamie@113: jamie@146: memset(result, 0, M * sizeof(double)); jamie@104: jamie@140: for (i = 0; i < L; i++) jamie@140: { jamie@104: jamie@104: /* Sum up this iteration's reflection coefficient. */ jamie@104: r = -data[i + 1]; jamie@140: for (j = 0; j < i; j++) jamie@104: r -= lpc[j] * data[i - j]; jamie@104: ref[i] = r /= error; jamie@104: jamie@104: /* Update LPC coefficients and total error. */ jamie@104: lpc[i] = r; jamie@140: for (j = 0; j < i / 2; j++) jamie@140: { jamie@146: double tmp = lpc[j]; jamie@104: lpc[j] = r * lpc[i - 1 - j]; jamie@104: lpc[i - 1 - j] += r * tmp; jamie@104: } jamie@104: if (i % 2) lpc[j] += lpc[j] * r; jamie@104: jamie@104: error *= 1 - r * r; jamie@104: } jamie@104: jamie@104: return XTRACT_SUCCESS; jamie@104: } jamie@104: jamie@146: int xtract_lpcc(const double *data, const int N, const void *argv, double *result) jamie@140: { jamie@104: jamie@104: /* Given N lpc coefficients extract an LPC cepstrum of size argv[0] */ jamie@104: /* Based on an an algorithm by rabiner and Juang */ jamie@104: jamie@104: int n, k; jamie@146: double sum; jamie@104: int order = N - 1; /* Eventually change this to Q = 3/2 p as suggested in Rabiner */ jamie@140: int cep_length; jamie@120: jamie@104: if(argv == NULL) jamie@115: cep_length = N - 1; /* FIX: if we're going to have default values, they should come from the descriptor */ jamie@104: else jamie@115: cep_length = *(int *)argv; jamie@146: //cep_length = (int)((double *)argv)[0]; jamie@104: jamie@146: memset(result, 0, cep_length * sizeof(double)); jamie@104: jamie@140: for (n = 1; n <= order && n <= cep_length; n++) jamie@140: { jamie@146: sum = 0.0; jamie@104: for (k = 1; k < n; k++) jamie@104: sum += k * result[k-1] * data[n - k]; jamie@104: result[n-1] = data[n] + sum / n; jamie@104: } jamie@104: jamie@104: /* be wary of these interpolated values */ jamie@140: for(n = order + 1; n <= cep_length; n++) jamie@140: { jamie@146: sum = 0.0; jamie@104: for (k = n - (order - 1); k < n; k++) jamie@104: sum += k * result[k-1] * data[n - k]; jamie@104: result[n-1] = sum / n; jamie@104: } jamie@104: jamie@104: return XTRACT_SUCCESS; jamie@104: jamie@104: } jamie@146: //int xtract_lpcc_s(const double *data, const int N, const void *argv, double *result){ jamie@104: // return XTRACT_SUCCESS; jamie@104: //} jamie@104: jamie@146: int xtract_subbands(const double *data, const int N, const void *argv, double *result) jamie@140: { jamie@104: jamie@114: int n, bw, xtract_func, nbands, scale, start, lower, *argi, rv; jamie@114: jamie@114: argi = (int *)argv; jamie@114: jamie@114: xtract_func = argi[0]; jamie@114: nbands = argi[1]; jamie@114: scale = argi[2]; jamie@114: start = argi[3]; jamie@114: jamie@114: if(scale == XTRACT_LINEAR_SUBBANDS) jamie@114: bw = floorf((N - start) / nbands); jamie@114: else jamie@114: bw = start; jamie@114: jamie@114: lower = start; jamie@115: rv = XTRACT_SUCCESS; jamie@114: jamie@140: for(n = 0; n < nbands; n++) jamie@140: { jamie@114: jamie@114: /* Bounds sanity check */ jamie@140: if(lower >= N || lower + bw >= N) jamie@140: { jamie@120: // printf("n: %d\n", n); jamie@146: result[n] = 0.0; jamie@114: continue; jamie@115: } jamie@114: jamie@114: rv = xtract[xtract_func](data+lower, bw, NULL, &result[n]); jamie@114: jamie@114: if(rv != XTRACT_SUCCESS) jamie@114: return rv; jamie@114: jamie@140: switch(scale) jamie@140: { jamie@140: case XTRACT_OCTAVE_SUBBANDS: jamie@140: lower += bw; jamie@140: bw = lower; jamie@140: break; jamie@140: case XTRACT_LINEAR_SUBBANDS: jamie@140: lower += bw; jamie@140: break; jamie@114: } jamie@114: jamie@114: } jamie@114: jamie@114: return rv; jamie@114: jamie@114: } jamie@114: jamie@114: jamie@114: