/*
 *  SDSPLD - Levinson-Durbin solver for real correlation functions.
 *  DSP Blockset S-Function to solve a symmetric Toeplitz system of
 *  equations using the Levinson-Durbin recursion.  Input is a vector
 *  of autocorrelation coefficients, starting with lag 0 as the first
 *  element.  Recursion order is length(input)-1.
 *
 *  Copyright 1995-2000 The MathWorks, Inc.
 *  $Revision: 1.10 $ $Date: 2000/06/14 14:27:58 $
 */

#define S_FUNCTION_NAME sdspld

#include "dsp_sim.h"

static void mdlInitializeSizes(SimStruct *S)
{
    ssSetNumSFcnParams(S, 0);

#if defined(MATLAB_MEX_FILE)
    if (ssGetNumSFcnParams(S) != ssGetSFcnParamsCount(S)) return;
#endif

    ssSetNumInputs(        S, DYNAMICALLY_SIZED);
    ssSetNumOutputs(       S, DYNAMICALLY_SIZED);
    ssSetDirectFeedThrough(S, 1);
    ssSetNumSampleTimes(   S, 1);
    ssSetOptions(          S, SS_OPTION_EXCEPTION_FREE_CODE);
}


static void mdlInitializeSampleTimes(SimStruct *S)
{
    ssSetSampleTimeEvent(S, 0, INHERITED_SAMPLE_TIME);
    ssSetOffsetTimeEvent(S, 0, 0.0);
}


static void mdlInitializeConditions(real_T *x0, SimStruct *S)
{
}


static void mdlOutputs(real_T *y, const real_T *x, const real_T *u, 
                       SimStruct *S, int_T tid)
{
    const int_T N = ssGetNumInputs(S);
    real_T      E = u[0];
    int_T       i;

    for(i=1; i<N; i++) {
        int_T  j;
        real_T ki = u[i];

        /* Update reflection coefficient: */
        for (j=1; j<i; j++) {
            ki += y[j] * u[i-j];
        }
        ki /= -E;
        E *= (1 - ki*ki);

        /* Update polynomial: */
        for (j=1; j<=(i-1)/2; j++) {
            real_T t = y[j];
            y[j]   += ki * y[i-j];
            y[i-j] += ki * t;
        }
        if (i%2 == 0) {
            y[i/2] *= 1+ki;
        }
        y[i] = y[N+i-1] = ki;  /* Record reflection coefficient */
    }
    y[0] = 1.0;
}


static void mdlUpdate(real_T *x, const real_T *u, SimStruct *S, int_T tid)
{
}


static void mdlDerivatives(real_T *dx, const real_T *x, const real_T *u, 
                           SimStruct *S, int_T tid)
{
}


static void mdlTerminate(SimStruct *S)
{
}


#ifdef MATLAB_MEX_FILE
#define MDL_GET_INPUT_PORT_WIDTH
/* Function: mdlGetInputPortWidth =============================================
 *
 */
static int_T mdlGetInputPortWidth(SimStruct *S, int_T outputPortWidth)
{
    if (outputPortWidth % 2 == 0) {
        THROW_ERROR(S, "Invalid output port width");
    }
    return((outputPortWidth+1)/2);
}

#define MDL_GET_OUTPUT_PORT_WIDTH
/* Function: mdlGetOutputPortWidth ============================================
 *
 */
static int_T mdlGetOutputPortWidth(SimStruct *S, int_T inputPortWidth)
{
    return(2*inputPortWidth-1);
}
#endif


#include "dsp_trailer.c"
 
/* [EOF]: sdspld.c */