real_peel(): maps LPC coefficients to waveguide radii

real_peel

/*----------------------------------------------------------------------------
    rt_lpc - real-time LPC analysis/synthesis and visualization

    Copyright (c) 2004 Ananya Misra, Ge Wang, Perry R. Cook.  All rights reserved.
      http://soundlab.cs.princeton.edu/

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
    U.S.A.
-----------------------------------------------------------------------------*/
//-----------------------------------------------------------------------------
// name: rt_lpc.cpp
// desc: real-time LPC analysis/synthesis/visualization, formerly iannith
//
// author: Ananya Misra ([email protected])
//         Ge Wang ([email protected])
//         Perry Cook ([email protected])
// library:
//         (STK) Synthesis ToolKit - Perry Cook and Gary Scavone
//         (FFT) CARL CMusic Distribution
// date: today
//-----------------------------------------------------------------------------

// .....some code .....

//------------------------------------------------------------------------------
// name: real_peel()
// desc: ...
//------------------------------------------------------------------------------
void real_peel( float * coef, int order, float * radii )
{
    int i,j;
    float rs, temp, tempcoeffs1[1024], tempcoeffs2[1024], reflections[1024];

    tempcoeffs1[0] = 1.0f;

    // initialize
    for( i = 0; i < order; i++ )
        tempcoeffs1[i+1] = -coef[i];

    for( i = 0; i < order; i++ )
    {
        for( j = 0; j <= order-i; j++ ) tempcoeffs2[order-i-j] = tempcoeffs1[j];
        reflections[i] = tempcoeffs2[0];
        rs = 1.0f - reflections[i]*reflections[i];
        if( rs <= 0.0f ) rs = 10000.0f; else rs = 1.0f / rs;
        for( j = 0; j < order-i; j++ )
            tempcoeffs1[j] = (tempcoeffs1[j] - reflections[i] * tempcoeffs2[j]) * rs;
    }

    radii[0] = 1.0f;
    for( i = 0; i <= order; i++ )
        radii[i+1] = radii[i] * sqrt((1-reflections[i])/(1.0f + reflections[i]));
    for( i = 0; i <= order; i++ )
        radii[i] = radii[i] / radii[order];
    for( i = 0; i <= order/2; i++ )
    {
        temp = radii[i];
        radii[i] = radii[order-i];
        radii[order-i] = temp;
    }
    radii[0] = .5f;
    //for( i = 0; i < order; i++ )
    //    printf("%f ",radii[i]);

    //printf("\n");
    //fflush(stdout);
}