## Signal Processing Library

Program name: mcepstra

Language: C

In file: parameter_lib.c

Objective: Calculate the Mel Cepstral Coefficients from the Cepstral vectors

Unfortunately this routine seems to have a slight problem.
** I based my work on what is explained in:
** http://limcserver.dee.ufcg.edu.br/~rrbrandt/dsp/lib_src/mcepstral.php
** I think my work is correct, but I am getting bad results.
** For my PhD work, I ended up abandoning this routine, and using a routine
** from Voicebox with MatLab. Maybe someday I can go back and find an error
** and fix this routine. If someone else finds an error, please point it out. Thanks.

Usage: void mel_cepstral(double vector_i[], double vector_o[], int order, double gain);

Parameters:

• vector_i - The Cepstral Coefficients input vector
• vector_o - The output array to put the calculated data
• order - The order
• gain - Gain to use in the routine

Mathematical Description:

`void mel_cepstral(double vector_i[], double vector_o[], int order, double gain){        int i, k, m, lpc_index, index;        double alpha;        double previous_am[MAX_ALLOC_SIZE];        double am[MAX_ALLOC_SIZE];        double K, sum;        /*        ** Information for this routine taken from:        ** "Recursive Calculation of Mel-Cepstrum from LP Coefficients"        ** By Keiichi Tokuda, Takao Kobayashi and Satoshi Imai        ** April 1994        */        alpha = 0.455;        K = gain;        if((order+1) >= MAX_ALLOC_SIZE) {                fprintf(stderr, "Please increase MAX_ALLOC_SIZE in window_lib.h to at least %d\n", order+2);                exit(1);        }        /* In equations 17, 18 and 19 in the article, we see the following variables:        ** a() = the array of LP coefficients denoted here by "vector_i"        ** a_tilde() = the current line of the matrix denoted here by "am"        ** a_tilde(i-1)() = the previous line of the matrix denoted here by "previous_am"        ** alpha = a constant denoted here by "alpha" and set to 0.455        ** K =        ** c_tilde() = output coefficient vector denoted here by "vector_o"        ** M = number of LP Coefficients denoted here by "p";        */        for(i = 0; i <= order; i++) {                /* start with a zero vector */                vector_o[i] = 0.0;                previous_am[i] = 0.0;        }        for(i = -order; i <= 0; i++) {                DEBUG_WRITE2(16, " loop i=%d\n", i);                /* added p to i into index so that index is always positive */                index = i + order;                if(index < 0) {                        fprintf(stderr, "Error with index (too small)\n");                        exit(1);                }                if(index > order) {                        fprintf(stderr, "Error with index (too large)\n");                        exit(1);                }                if(i < 0) lpc_index = -i;                if(i == 0) lpc_index = 0;                if(lpc_index < 0) {                        fprintf(stderr, "Error with lpc_index (too small)\n");                        exit(1);                }                if(lpc_index > order) {                        fprintf(stderr, "Error with lpc_index (too large)\n");                        exit(1);                }                if(index == 0) {                        am = vector_i[lpc_index];                        for(m = 1; m <= order; m++) {                                am = 0;                        }                } else if(index == 1) {                        am = vector_i[lpc_index] + alpha * previous_am;                        am = (1 - (alpha * alpha)) * previous_am;                        for(m = 2; m <= order; m++) {                                am[m] = alpha * (-am[m-1]);                        }                } else {                        am = vector_i[lpc_index] + alpha * previous_am;                        am = (1 - (alpha * alpha)) * previous_am + alpha * previous_am;                        for(m = 2; m <= order; m++) {                                am[m] = previous_am[m-1] + alpha * (previous_am[m] - am[m-1]);                        }                }                for(k = 0; k <= order; k++ ) {                        previous_am[k] = am[k];                }        }        vector_o = log(K/ am);        for(m = 1; m <= order; m++) {                sum = 0.0;                for(k = 1; k <= (m - 1); k++) {                        sum += ((double) k / (double) m) * vector_o[k] * am[m-k];                }                vector_o[m] = -am[m] - sum;        }} /* End of the mel cepstral function */`