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@30: jamie@140: #include "fftsg.h" jamie@140: jamie@98: #include "xtract/libxtract.h" jamie@98: #include "xtract_macros_private.h" jamie@98: jamie@85: #ifndef roundf jamie@140: float roundf(float f) jamie@140: { jamie@140: if (f - (int)f >= 0.5) jamie@140: return (float)((int)f + 1); jamie@140: else jamie@140: return (float)((int)f); jamie@140: } jamie@85: #endif jamie@85: jamie@113: #ifndef powf jamie@120: #define powf pow jamie@113: #endif jamie@113: jamie@113: #ifndef expf jamie@120: #define expf exp jamie@113: #endif jamie@113: jamie@113: #ifndef sqrtf jamie@120: #define sqrtf sqrt jamie@113: #endif jamie@113: jamie@113: #ifndef fabsf jamie@120: #define fabsf fabs jamie@113: #endif jamie@113: jamie@98: #include "xtract_globals_private.h" jamie@98: #include "xtract_macros_private.h" jamie@1: jamie@140: int xtract_spectrum(const float *data, const int N, const void *argv, float *result) jamie@140: { jamie@1: jamie@140: int vector = 0; jamie@140: int withDC = 0; jamie@140: int normalise = 0; jamie@140: float q = 0.f; jamie@140: float temp = 0.f; jamie@140: float max = 0.f; jamie@140: float NxN = XTRACT_SQ(N); jamie@140: float *marker = NULL; jamie@140: size_t bytes = N * sizeof(double); jamie@140: double *rfft = NULL; jamie@140: unsigned int n = 0; jamie@140: unsigned int m = 0; jamie@140: unsigned int nx2 = 0; jamie@140: unsigned int M = N >> 1; jamie@98: jamie@140: rfft = (double *)malloc(bytes); jamie@140: //memcpy(rfft, data, bytes); jamie@1: jamie@140: for(n = 0; n < N; ++n){ jamie@140: rfft[n] = (double)data[n]; jamie@140: } jamie@1: jamie@56: q = *(float *)argv; jamie@54: vector = (int)*((float *)argv+1); jamie@70: withDC = (int)*((float *)argv+2); jamie@105: normalise = (int)*((float *)argv+3); jamie@105: jamie@56: XTRACT_CHECK_q; jamie@46: jamie@140: if(!ooura_data_spectrum.initialised) 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@140: /* ooura is in-place jamie@140: * the output format seems to be jamie@140: * a[0] - DC, a[1] - nyquist, a[2...N-1] - remaining bins jamie@140: */ jamie@140: rdft(N, 1, rfft, ooura_data_spectrum.ooura_ip, jamie@140: ooura_data_spectrum.ooura_w); 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@140: if(!withDC && n == 0) jamie@140: { jamie@140: continue; jamie@140: } jamie@140: nx2 = n * 2; jamie@140: temp = XTRACT_SQ(rfft[nx2]) + XTRACT_SQ(rfft[nx2+1]); jamie@140: if (temp > XTRACT_LOG_LIMIT) jamie@140: { jamie@140: temp = logf(sqrtf(temp) / (float)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@140: if(!withDC && n == 0) jamie@140: { jamie@140: ++n; jamie@120: } jamie@140: result[m] = (XTRACT_SQ(rfft[n]) + XTRACT_SQ(rfft[N - n])) / 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@140: if(!withDC && n == 0) jamie@140: { jamie@140: ++n; jamie@120: } jamie@140: if ((temp = XTRACT_SQ(rfft[n]) + XTRACT_SQ(rfft[N - n])) > jamie@140: XTRACT_LOG_LIMIT) jamie@140: temp = logf(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@140: marker = &result[M-1]; jamie@120: } jamie@140: if(n==1 && withDC) /* discard Nyquist */ jamie@140: { jamie@140: ++n; jamie@140: } jamie@120: jamie@140: *marker = (float)(sqrt(XTRACT_SQ(rfft[n*2]) + jamie@140: XTRACT_SQ(rfft[n*2+1])) / (double)N); jamie@111: jamie@140: XTRACT_SET_FREQUENCY; jamie@140: XTRACT_GET_MAX; jamie@140: 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@140: free(rfft); jamie@120: jamie@56: return XTRACT_SUCCESS; jamie@1: } jamie@1: jamie@140: int xtract_autocorrelation_fft(const float *data, const int N, const void *argv, float *result) jamie@140: { jamie@120: jamie@140: double *rfft = NULL; jamie@140: int n = 0; jamie@140: int M = 0; 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@140: memcpy(rfft, data, N * sizeof(float)); jamie@75: jamie@140: rdft(M, 1, rfft, ooura_data_autocorrelation_fft.ooura_ip, jamie@140: ooura_data_autocorrelation_fft.ooura_w); jamie@1: jamie@140: 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@140: rfft[n*2+1] = 0.f; 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@75: /* Normalisation factor */ jamie@75: M = M * N; jamie@75: jamie@75: for(n = 0; n < N; n++) jamie@140: result[n] = rfft[n] / (float)M; jamie@75: jamie@140: free(rfft); jamie@38: jamie@56: return XTRACT_SUCCESS; jamie@1: } jamie@1: jamie@140: int xtract_mfcc(const float *data, const int N, const void *argv, float *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: { danstowell@68: result[filter] = 0.f; jamie@140: for(n = 0; n < N; n++) jamie@140: { jamie@71: result[filter] += data[n] * f->filters[filter][n]; jamie@30: } jamie@113: result[filter] = logf(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@140: int xtract_dct(const float *data, const int N, const void *argv, float *result) jamie@140: { jamie@120: jamie@140: int n; jamie@140: int m; jamie@140: float *temp = calloc(N, sizeof(float)); jamie@120: jamie@140: for (n = 0; n < N; ++n) jamie@140: { jamie@140: for(m = 1; m <= N; ++m) { jamie@140: temp[n] += data[m - 1] * cos(M_PI * (n / (float)N) * (m - 0.5)); jamie@140: } jamie@140: } jamie@120: jamie@140: return XTRACT_SUCCESS; jamie@30: } jamie@30: jamie@140: int xtract_autocorrelation(const float *data, const int N, const void *argv, float *result) jamie@140: { jamie@30: jamie@30: /* Naive time domain implementation */ jamie@120: jamie@30: int n = N, i; jamie@120: jamie@30: float 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@140: int xtract_amdf(const float *data, const int N, const void *argv, float *result) jamie@140: { jamie@1: jamie@1: int n = N, i; jamie@120: jamie@6: float md, temp; jamie@1: jamie@140: while(n--) jamie@140: { jamie@120: md = 0.f; 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@113: result[n] = md / (float)N; jamie@1: } jamie@38: jamie@56: return XTRACT_SUCCESS; jamie@1: } jamie@1: jamie@140: int xtract_asdf(const float *data, const int N, const void *argv, float *result) jamie@140: { jamie@120: jamie@1: int n = N, i; jamie@120: jamie@1: float sd; jamie@1: jamie@140: while(n--) jamie@140: { jamie@120: sd = 0.f; 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@113: result[n] = sd / (float)N; jamie@1: } jamie@38: jamie@56: return XTRACT_SUCCESS; jamie@1: } jamie@1: jamie@140: int xtract_bark_coefficients(const float *data, const int N, const void *argv, float *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@110: result[band] = 0.f; 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@140: int xtract_peak_spectrum(const float *data, const int N, const void *argv, float *result) jamie@140: { jamie@1: jamie@56: float 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@56: threshold = max = y = y2 = y3 = p = q = 0.f; jamie@120: jamie@140: if(argv != NULL) jamie@140: { jamie@56: q = ((float *)argv)[0]; jamie@55: threshold = ((float *)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@98: input = (float *)calloc(N, sizeof(float)); jamie@98: jamie@98: bytes = N * sizeof(float); jamie@43: jamie@43: if(input != NULL) jamie@120: input = 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@140: result[N + n] = q * (n + (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@43: free(input); jamie@56: return (rv ? rv : XTRACT_SUCCESS); jamie@1: } jamie@120: jamie@140: int xtract_harmonic_spectrum(const float *data, const int N, const void *argv, float *result) jamie@140: { jamie@120: jamie@140: int n = (N >> 1), M = n; jamie@38: jamie@43: const float *freqs, *amps; jamie@55: float f0, threshold, ratio, nearest, distance; jamie@38: jamie@52: amps = data; jamie@52: freqs = data + n; jamie@38: f0 = *((float *)argv); jamie@55: threshold = *((float *)argv+1); jamie@38: jamie@38: ratio = nearest = distance = 0.f; jamie@38: jamie@140: while(n--) jamie@140: { jamie@140: if(freqs[n]) jamie@140: { jamie@120: ratio = freqs[n] / f0; jamie@120: nearest = roundf(ratio); jamie@120: distance = fabs(nearest - ratio); jamie@120: if(distance > threshold) jamie@120: result[n] = result[M + n] = 0.f; 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@120: result[n] = result[M + n] = 0.f; jamie@38: } jamie@56: return XTRACT_SUCCESS; jamie@38: } jamie@120: jamie@140: int xtract_lpc(const float *data, const int N, const void *argv, float *result) jamie@140: { jamie@104: jamie@104: int i, j, k, M, L; jamie@140: float r = 0.f, jamie@104: error = 0.f; jamie@104: jamie@104: float *ref = NULL, jamie@140: *lpc = NULL ; jamie@104: jamie@104: error = data[0]; jamie@104: k = N; /* The length of *data */ 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@113: memset(result, 0, M * sizeof(float)); jamie@104: return XTRACT_NO_RESULT; jamie@104: } jamie@113: jamie@104: memset(result, 0, M * sizeof(float)); 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@104: float 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@140: int xtract_lpcc(const float *data, const int N, const void *argv, float *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@104: float 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@120: //cep_length = (int)((float *)argv)[0]; jamie@104: jamie@104: memset(result, 0, cep_length * sizeof(float)); jamie@104: jamie@140: for (n = 1; n <= order && n <= cep_length; n++) jamie@140: { jamie@104: sum = 0.f; 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@104: sum = 0.f; 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@104: //int xtract_lpcc_s(const float *data, const int N, const void *argv, float *result){ jamie@104: // return XTRACT_SUCCESS; jamie@104: //} jamie@104: jamie@140: int xtract_subbands(const float *data, const int N, const void *argv, float *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@115: result[n] = 0.f; 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: