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:

Comments:

User Comments

Code:

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[0] = vector_i[lpc_index];
                        for(m = 1; m <= order; m++) {
                                am[0] = 0;
                        }
                } else if(index == 1) {
                        am[0] = vector_i[lpc_index] + alpha * previous_am[0];
                        am[1] = (1 - (alpha * alpha)) * previous_am[0];
                        for(m = 2; m <= order; m++) {
                                am[m] = alpha * (-am[m-1]);
                        }
                } else {
                        am[0] = vector_i[lpc_index] + alpha * previous_am[0];
                        am[1] = (1 - (alpha * alpha)) * previous_am[0] + alpha * previous_am[1];
                        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[0] = log(K/ am[0]);
        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 */