/*
 *  SDSPLDC - Levinson-Durbin solver for complex correlation functions.
 *  DSP Blockset S-Function to solve a Hermitian 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.11 $ $Date: 2000/06/14 14:27:59 $
 */

#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   Nin = ssGetNumInputs(S);
    const int_T   N   = Nin/2;
    const real_T *ur  = u;
    const real_T *ui  = u + N;
    real_T       *yr  = y;
    real_T       *ykr = yr  + N;
    real_T       *yi  = ykr + N-1;
    real_T       *yki = yi  + N;
    real_T        E   = ur[0]; /* Assume lag-0 is purely real */
    int_T         i;

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

        for (j=1; j<i; j++) {
            /* k = y * r[reverse order] */
            kr += yr[j] * ur[i-j] - yi[j] * ui[i-j];
            ki += yi[j] * ur[i-j] + yr[j] * ui[i-j];
        }
        kr /= -E;    
        ki /= -E;
        E *= (1 - kr*kr - ki*ki);

        /* Update polynomial: */
        for (j=1; j<=(i-1)/2; j++) {
            real_T tr = yr[j];
            real_T ti = yi[j];

            /* ynew = yold + k * conj(yold[reverse order]) */
            yr[j]   += kr * yr[i-j] + ki * yi[i-j];
            yi[j]   += ki * yr[i-j] - kr * yi[i-j];

            yr[i-j] += kr * tr + ki * ti;
            yi[i-j] += ki * tr - kr * ti;
        }
        if (i%2 == 0) {
            real_T tr = yr[i/2];
            real_T ti = yi[i/2];
            yr[i/2] += kr * tr + ki * ti;
            yi[i/2] += ki * tr - kr * ti;
        } 

        yr[i] = ykr[i-1] = kr;
        yi[i] = yki[i-1] = ki;
    }
    yr[0] = 1.0;
    yi[0] = 0.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) {
        ssSetErrorStatus(S, "Invalid output port width");
        return(1);
    }
    return((outputPortWidth+2)/2);
}

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


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