Mercurial > hg > libxtract
changeset 104:a32738e9d955
- Fixes to descriptors.c where no break statement was given for certain cases is switch conditionals
- Added LPC and LPCC extraction functions. LPC implements Durbin method as described in Rabiner and Juang and implemented in Dr. Dobbs 1994 edition by Jutta Degener
| author | Jamie Bullock <jamie@postlude.co.uk> |
|---|---|
| date | Mon, 24 Dec 2007 13:21:13 +0000 |
| parents | 1cbbe5b5e461 |
| children | f2af1c75e3ed |
| files | TODO examples/puredata/lpcc-test.pd examples/puredata/xtract~.c src/descriptors.c src/libxtract.c src/vector.c xtract/libxtract.h xtract/xtract_vector.h |
| diffstat | 8 files changed, 255 insertions(+), 7 deletions(-) [+] |
line wrap: on
line diff
--- a/TODO Fri Dec 21 11:05:20 2007 +0000 +++ b/TODO Mon Dec 24 13:21:13 2007 +0000 @@ -13,3 +13,5 @@ ...do other stuff and eventually... ...optimise! Use the fftw guru interface to create multipurpose global plans from xtract_fft_init() +Add LPC via Schur algorithm +If argv uses an integer, then it should point to a variable of type int. At the moment floats are used to store ints.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/examples/puredata/lpcc-test.pd Mon Dec 24 13:21:13 2007 +0000 @@ -0,0 +1,57 @@ +#N canvas 356 233 722 300 10; +#X obj 160 170 xtract~ autocorrelation; +#N canvas 0 0 450 300 graph1 0; +#X array lpcc 16 float 5; +#A 0 0 0 -0.827113 -0.827113 -0.827113 -0.485055 0.314974 0.657032 +0.810476 0.396513 0.0440753 -0.630112 -0.9456 -0.887411 -1.12995 0.596854 +; +#X coords 0 16 15 -16 200 140 1; +#X restore 458 97 graph; +#X obj 149 34 noise~; +#N canvas 9 103 450 300 lpc 0; +#X obj 323 25 block~ 32; +#X obj 140 52 inlet~; +#X obj 140 113 xtract~ lpc; +#N canvas 237 125 450 300 lpcc 0; +#X obj 26 130 xtract~ lpcc; +#N canvas 0 0 450 300 graph1 0; +#X array lpc 16 float 5; +#A 0 -1.03207 -1.03207 -1.03207 -21.8757 -21.8757 -21.8757 -21.8757 +-0.391539 -0.391539 -0.391539 -0.391539 0.448536 0.448536 0.448536 +0.448536 0.516603; +#X coords 0 16 15 -16 200 140 1; +#X restore 219 37 graph; +#X obj 23 12 block~ 16; +#X obj 26 51 tabreceive~ lpc; +#X msg 95 107 list 16; +#X obj 24 196 tabsend~ lpcc; +#X obj 95 85 loadbang; +#X connect 0 0 5 0; +#X connect 3 0 0 0; +#X connect 4 0 0 1; +#X connect 6 0 4 0; +#X restore 298 172 pd lpcc; +#X obj 137 220 tabsend~ lpc; +#X obj 139 152 a_blockswap~ 32; +#X connect 1 0 2 0; +#X connect 2 0 5 0; +#X connect 5 0 4 0; +#X restore 160 218 pd lpc; +#X obj 31 52 a_hann 16; +#X obj 31 25 loadbang; +#X obj 98 23 bng 15 250 50 0 empty empty empty 0 -6 0 10 -262144 -1 +-1; +#X obj 156 139 *~; +#X obj 27 -41 block~ 16 4; +#X obj 265 81 *~ 0.4; +#X obj 148 87 *~ 0.1; +#X obj 266 46 osc~ 1000; +#X connect 0 0 3 0; +#X connect 2 0 10 0; +#X connect 4 0 7 0; +#X connect 5 0 4 0; +#X connect 6 0 4 0; +#X connect 7 0 0 0; +#X connect 9 0 7 1; +#X connect 10 0 7 1; +#X connect 11 0 9 0;
--- a/examples/puredata/xtract~.c Fri Dec 21 11:05:20 2007 +0000 +++ b/examples/puredata/xtract~.c Mon Dec 24 13:21:13 2007 +0000 @@ -93,7 +93,7 @@ tmp_in = copybytes(in, N * sizeof(t_float)); tmp_out = getbytes(N * sizeof(t_float)); - if(x->feature == XTRACT_PEAK_SPECTRUM) + if(x->feature == XTRACT_PEAK_SPECTRUM || x->feature == XTRACT_LPC) N >>= 1; return_code = xtract[x->feature](tmp_in, N, x->argv, tmp_out); @@ -239,7 +239,9 @@ x->feature == XTRACT_BARK_COEFFICIENTS || x->feature == XTRACT_SPECTRUM || x->feature == XTRACT_PEAK_SPECTRUM || - x->feature == XTRACT_HARMONIC_SPECTRUM) + x->feature == XTRACT_HARMONIC_SPECTRUM || + x->feature == XTRACT_LPC || + x->feature == XTRACT_LPCC) x->feature_type = XTRACT_VECTOR; else if (x->feature == XTRACT_FLUX || x->feature == XTRACT_ATTACK_TIME || @@ -273,8 +275,8 @@ x->argv = getbytes(argc * sizeof(float)); - while(argc--) - ((t_float *)x->argv)[argc] = atom_getfloat(&argv[argc]); + while(argc--) + ((t_float *)x->argv)[argc] = atom_getfloat(&argv[argc]); /* }*/ }
--- a/src/descriptors.c Fri Dec 21 11:05:20 2007 +0000 +++ b/src/descriptors.c Mon Dec 24 13:21:13 2007 +0000 @@ -137,21 +137,26 @@ case XTRACT_LOWEST_VALUE: case XTRACT_TONALITY: case XTRACT_MFCC: + case XTRACT_LPC: + case XTRACT_LPCC: *argv_min = XTRACT_ANY; *argv_max = XTRACT_ANY; *argv_def = XTRACT_ANY; *argv_unit = XTRACT_ANY; + break; case XTRACT_SPECTRAL_INHARMONICITY: *argv_min = 0.f; *argv_max = XTRACT_SR_UPPER_LIMIT / 2; *argv_def = XTRACT_FUNDAMENTAL_DEFAULT; *argv_unit = XTRACT_HERTZ; + break; case XTRACT_F0: case XTRACT_FAILSAFE_F0: *argv_min = XTRACT_SR_LOWER_LIMIT; *argv_max = XTRACT_SR_UPPER_LIMIT; *argv_def = XTRACT_SR_DEFAULT; *argv_unit = XTRACT_HERTZ; + break; /* argc = 2 */; case XTRACT_ROLLOFF: *argv_min = XTRACT_FFT_BANDS_MIN; @@ -162,6 +167,7 @@ *(argv_max + 1) = 100.f; *(argv_def + 1) = 95.f; *(argv_unit + 1) = XTRACT_PERCENT; + break; case XTRACT_SPECTRUM: *argv_min = XTRACT_SR_LOWER_LIMIT / 2; *argv_max = XTRACT_SR_UPPER_LIMIT / 2; @@ -171,6 +177,7 @@ *(argv_max + 1) = 3 ; *(argv_def + 1) = 0; *(argv_unit + 1) = XTRACT_NONE; + break; case XTRACT_PEAK_SPECTRUM: *argv_min = XTRACT_SR_LOWER_LIMIT / 2; *argv_max = XTRACT_SR_UPPER_LIMIT / 2; @@ -180,6 +187,7 @@ *(argv_max + 1) = 100.f ; *(argv_def + 1) = 10.f ; *(argv_unit + 1) = XTRACT_PERCENT; + break; case XTRACT_HARMONIC_SPECTRUM: *argv_min = 0.f; *argv_max = XTRACT_SR_UPPER_LIMIT / 2; @@ -189,6 +197,7 @@ *(argv_max + 1) = 1.f ; *(argv_def + 1) = .1f ; *(argv_unit + 1) = XTRACT_NONE; + break; case XTRACT_NOISINESS: case XTRACT_SKEWNESS: case XTRACT_KURTOSIS: @@ -203,6 +212,7 @@ *(argv_max + 1) = XTRACT_NONE; *(argv_def + 1) = XTRACT_NONE; *(argv_unit + 1) = XTRACT_NONE; + break; case XTRACT_BARK_COEFFICIENTS: /* BARK_COEFFICIENTS is special because argc = BARK_BANDS */ default: @@ -210,6 +220,7 @@ *argv_max = XTRACT_NONE; *argv_def = XTRACT_NONE; *argv_unit = XTRACT_NONE; + break; } argv_donor = &d->argv.donor[0]; @@ -326,6 +337,12 @@ case XTRACT_MFCC: *data_format = XTRACT_SPECTRAL_MAGNITUDES; break; + case XTRACT_LPC: + *data_format = XTRACT_AUTOCORRELATION_COEFFS; + break; + case XTRACT_LPCC: + *data_format = XTRACT_LPC_COEFFS; + break; case XTRACT_SPECTRAL_INHARMONICITY: case XTRACT_HARMONIC_SPECTRUM: *data_format = XTRACT_SPECTRAL_PEAKS; @@ -414,6 +431,8 @@ case XTRACT_TONALITY: case XTRACT_LOUDNESS: case XTRACT_NONZERO_COUNT: + case XTRACT_LPC: + case XTRACT_LPCC: *data_unit = XTRACT_ANY; break; case XTRACT_SPECTRAL_MEAN: @@ -646,6 +665,25 @@ strcpy(desc, "Extract MFCC from a spectrum"); strcpy(p_desc, "Extract MFCC from an audio spectrum"); strcpy(author, "Rabiner"); + break; + case XTRACT_LPC: + strcpy(name, "lpc"); + strcpy(p_name, "Linear predictive coding coefficients"); + strcpy(desc, "Extract LPC from autocorrelation coefficients"); + strcpy(p_desc, + "Extract LPC from autocorrelation coefficients"); + strcpy(author, + "Rabiner and Juang as implemented by Jutta Degener"); + *year = 1994; + break; + case XTRACT_LPCC: + strcpy(name, "lpcc"); + strcpy(p_name, "Linear predictive coding cepstral coefficients"); + strcpy(desc, "Extract LPC cepstrum from LPC coefficients"); + strcpy(p_desc, + "Extract LPC cepstrum from LPC coefficients"); + strcpy(author, "Rabiner and Juang"); + *year = 1993; break; case XTRACT_BARK_COEFFICIENTS: strcpy(name, "bark_coefficients"); @@ -926,6 +964,10 @@ *argc = 1; *argv_type = XTRACT_MEL_FILTER; break; + case XTRACT_LPCC: + *argc = 1; + *argv_type = XTRACT_INT; + break; case XTRACT_BARK_COEFFICIENTS: *argc = XTRACT_BARK_BANDS; *argv_type = XTRACT_INT; @@ -961,6 +1003,7 @@ case XTRACT_ASDF: case XTRACT_NONZERO_COUNT: case XTRACT_ODD_EVEN_RATIO: + case XTRACT_LPC: default: *argc = 0; break; @@ -1019,6 +1062,8 @@ case XTRACT_SPECTRUM: case XTRACT_AUTOCORRELATION_FFT: case XTRACT_MFCC: + case XTRACT_LPC: + case XTRACT_LPCC: case XTRACT_DCT: case XTRACT_HARMONIC_SPECTRUM: *is_scalar = XTRACT_FALSE; @@ -1077,14 +1122,17 @@ *result_unit = XTRACT_HERTZ; *result_min = 0.f; *result_max = XTRACT_SR_UPPER_LIMIT / 2; + break; case XTRACT_ZCR: *result_unit = XTRACT_HERTZ; *result_min = 0.f; *result_max = XTRACT_ANY; + break; case XTRACT_ODD_EVEN_RATIO: *result_unit = XTRACT_NONE; *result_min = 0.f; *result_max = 1.f; + break; case XTRACT_LOUDNESS: case XTRACT_FLATNESS: case XTRACT_TONALITY: @@ -1093,10 +1141,13 @@ case XTRACT_POWER: case XTRACT_SHARPNESS: case XTRACT_SPECTRAL_SLOPE: + case XTRACT_LPC: + case XTRACT_LPCC: default: *result_unit = XTRACT_UNKNOWN; *result_min = XTRACT_UNKNOWN; *result_max = XTRACT_UNKNOWN; + break; } } else { @@ -1113,18 +1164,31 @@ case XTRACT_DCT: *result_format = XTRACT_ARBITRARY_SERIES; *result_unit = XTRACT_ANY; + break; case XTRACT_BARK_COEFFICIENTS: *result_format = XTRACT_BARK_COEFFS; *result_unit = XTRACT_UNKNOWN; /* FIX: check */ + break; case XTRACT_PEAK_SPECTRUM: case XTRACT_SPECTRUM: case XTRACT_HARMONIC_SPECTRUM: *result_format = XTRACT_SPECTRAL; *result_unit = XTRACT_ANY_AMPLITUDE_HERTZ; + break; case XTRACT_AUTOCORRELATION_FFT: + break; case XTRACT_MFCC: *result_format = XTRACT_MEL_COEFFS; *result_unit = XTRACT_UNKNOWN; /* FIX: check */ + break; + case XTRACT_LPC: + *result_format = XTRACT_LPC_COEFFS; + *result_unit = XTRACT_UNKNOWN; + break; + case XTRACT_LPCC: + *result_format = XTRACT_LPCC_COEFFS; + *result_unit = XTRACT_UNKNOWN; + break; default: break; }
--- a/src/libxtract.c Fri Dec 21 11:05:20 2007 +0000 +++ b/src/libxtract.c Mon Dec 24 13:21:13 2007 +0000 @@ -79,6 +79,8 @@ xtract_autocorrelation_fft, xtract_mfcc, xtract_dct, - xtract_harmonic_spectrum + xtract_harmonic_spectrum, + xtract_lpc, + xtract_lpcc };
--- a/src/vector.c Fri Dec 21 11:05:20 2007 +0000 +++ b/src/vector.c Mon Dec 24 13:21:13 2007 +0000 @@ -449,3 +449,90 @@ return XTRACT_SUCCESS; } +int xtract_lpc(const float *data, const int N, const void *argv, float *result){ + + int i, j, k, M, L; + float r = 0.f, + error = 0.f; + + float *ref = NULL, + *lpc = NULL ; + + error = data[0]; + k = N; /* The length of *data */ + L = N - 1; /* The number of LPC coefficients */ + M = L * 2; /* The length of *result */ + ref = result; + lpc = result+L; + + if(error == 0.0){ + for(i = 0; i < M; i++) + result[i] = 0.f; + return XTRACT_NO_RESULT; + } + + memset(result, 0, M * sizeof(float)); + + for (i = 0; i < L; i++) { + + /* Sum up this iteration's reflection coefficient. */ + r = -data[i + 1]; + for (j = 0; j < i; j++) + r -= lpc[j] * data[i - j]; + ref[i] = r /= error; + + /* Update LPC coefficients and total error. */ + lpc[i] = r; + for (j = 0; j < i / 2; j++) { + float tmp = lpc[j]; + lpc[j] = r * lpc[i - 1 - j]; + lpc[i - 1 - j] += r * tmp; + } + if (i % 2) lpc[j] += lpc[j] * r; + + error *= 1 - r * r; + } + + return XTRACT_SUCCESS; +} + +int xtract_lpcc(const float *data, const int N, const void *argv, float *result){ + + /* Given N lpc coefficients extract an LPC cepstrum of size argv[0] */ + /* Based on an an algorithm by rabiner and Juang */ + + int n, k; + float sum; + int order = N - 1; /* Eventually change this to Q = 3/2 p as suggested in Rabiner */ + int cep_length; + + if(argv == NULL) + cep_length = N - 1; + else + cep_length = (int)((float *)argv)[0]; + + memset(result, 0, cep_length * sizeof(float)); + + for (n = 1; n <= order && n <= cep_length; n++){ + sum = 0.f; + for (k = 1; k < n; k++) + sum += k * result[k-1] * data[n - k]; + result[n-1] = data[n] + sum / n; + } + + /* be wary of these interpolated values */ + for(n = order + 1; n <= cep_length; n++){ + sum = 0.f; + for (k = n - (order - 1); k < n; k++) + sum += k * result[k-1] * data[n - k]; + result[n-1] = sum / n; + } + + return XTRACT_SUCCESS; + +} +//int xtract_lpcc_s(const float *data, const int N, const void *argv, float *result){ +// return XTRACT_SUCCESS; +//} + +
--- a/xtract/libxtract.h Fri Dec 21 11:05:20 2007 +0000 +++ b/xtract/libxtract.h Mon Dec 24 13:21:13 2007 +0000 @@ -56,7 +56,7 @@ * @{ */ -#define XTRACT_FEATURES 54 +#define XTRACT_FEATURES 56 /** \brief Enumeration of features, elements are used as indixes to an array of pointers to feature extracton functions */ enum xtract_features_ { @@ -113,7 +113,9 @@ XTRACT_AUTOCORRELATION_FFT, XTRACT_MFCC, XTRACT_DCT, - XTRACT_HARMONIC_SPECTRUM + XTRACT_HARMONIC_SPECTRUM, + XTRACT_LPC, + XTRACT_LPCC }; /** \brief Enumeration of feature initialisation functions */ @@ -198,9 +200,12 @@ XTRACT_SPECTRAL_HARMONICS_MAGNITUDES, /* N spectral harmonic frequencies */ XTRACT_SPECTRAL_HARMONICS_FREQUENCIES, + XTRACT_AUTOCORRELATION_COEFFS, XTRACT_ARBITRARY_SERIES, XTRACT_AUDIO_SAMPLES, XTRACT_MEL_COEFFS, + XTRACT_LPC_COEFFS, + XTRACT_LPCC_COEFFS, XTRACT_BARK_COEFFS, XTRACT_NO_DATA } xtract_vector_t;
--- a/xtract/xtract_vector.h Fri Dec 21 11:05:20 2007 +0000 +++ b/xtract/xtract_vector.h Mon Dec 24 13:21:13 2007 +0000 @@ -126,6 +126,35 @@ */ int xtract_harmonic_spectrum(const float *data, const int N, const void *argv, float *result); +/** \brief Extract Linear Predictive Coding Coefficients + * + * Based on algorithm in Rabiner and Juang as implemented by Jutta Degener in Dr. Dobb's Journal December, 1994. + * + * Returns N-1 reflection (PARCOR) coefficients and N-1 LPC coefficients via *result + * + * \param *data: N autocorrelation values e.g the data pointed to by *result from xtract_autocorrelation() + * \param N: the number of autocorrelation values to be considered + * \param *argv: a pointer to NULL + * \param *result: a pointer to an array containing N-1 reflection coefficients and N-1 LPC coefficients. + * + * An array of size 2 * (N - 1) must be allocated, and *result must point to its first element. + */ +int xtract_lpc(const float *data, const int N, const void *argv, float *result); + +/** \brief Extract Linear Predictive Coding Cepstral Coefficients + * + * \param *data: a pointer to the first element in an array of LPC coeffiecients e.g. a pointer to the second half of the array pointed to by *result from xtract_lpc() + * \param N: the number of LPC coefficients to be considered + * \param *argv: a pointer to a float representing the order of the result vector. This must be a whole number. According to Rabiner and Juang the ratio between the number (p) of LPC coefficients and the order (Q) of the LPC cepstrum is given by Q ~ (3/2)p where Q > p. + * \param *result: a pointer to an array containing the resultant LPCC. + * + * An array of size Q, where Q is given by argv[0] must be allocated, and *result must point to its first element. + * + */ +int xtract_lpcc(const float *data, const int N, const void *argv, float *result); + + + /** @} */ #ifdef __cplusplus