/*
 * SDSPSTAT DSP Blockset S-Function for some statistics functions:
 * mean, variance, and RMS.
 *
 * Function Type: 0=mean, 1=var, 2=std, 3=RMS
 *
 *  Updated: D. Orofino, 25-Jun-97
 *  Copyright 1995-2000 The MathWorks, Inc.
 *  $Revision: 1.7 $  $Date: 2000/05/05 20:16:07 $
 */

#define S_FUNCTION_NAME sdspstat

#include <math.h>

#include "dsp_sim.h"

#define FCN_TYPE ssGetArg(S,0)
#define NUM_ARGS 1

typedef enum {
    fcnMean,
    fcnVar,
    fcnStdDev,
    fcnRMS
} FcnType;

#ifdef MATLAB_MEX_FILE
#define MDL_CHECK_PARAMETERS
static void mdlCheckParameters(SimStruct *S)
{
    /*
     * One parameter, a scalar double (0 or 1, only):
     */
    if ( !mxIsDouble(FCN_TYPE)    ||
         (mxGetM(FCN_TYPE) != 1)  ||
         (mxGetN(FCN_TYPE) != 1)  ||
         (mxGetPr(FCN_TYPE)[0]<0) ||
         (mxGetPr(FCN_TYPE)[0]>3)  ) {
        ssSetErrorStatus(S, "Function type must be a scalar in the range [0,3].");
    }
}
#endif


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

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

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


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)
{
    FcnType   ftype = (FcnType)((int_T)mxGetPr(FCN_TYPE)[0]);
    UPtrsType uptr  = ssGetUPtrs(S);
    int_T     N     = ssGetNumInputs(S);
    int_T     i;

    switch(ftype) {
    case fcnMean:
        {
            /* Compute sum(x)/length(x) */
            real_T sum = **(uptr++);
            for (i=1; i<N; i++) {
                sum += **(uptr++);
            }
            *y = sum / N;
            break;
        }

    case fcnVar:
    case fcnStdDev:
        {
            /* Compute variance or standard deviation:
             *   var = (x - mean(x)).^2 / (N-1),
             * or equivalently,
             *   var = (sum(x.^2) - sum(x)*sum(x)/N) / (N-1)
             *
             * Either way,
             *   std = sqrt(var)
             */
            real_T sx  = 0.0;
            real_T sx2 = 0.0;

            if (N == 1) {
                *y = 0.0;
            } else {
                for (i=0; i<N; i++) {
                    real_T v = **(uptr++);
                    sx  += v;
                    sx2 += v*v;
                }
                *y = (sx2 - sx*sx/N) / (N-1);

                if (ftype == fcnStdDev) {
                    *y = sqrt(*y);
                }
            }
            break;
        }

    case fcnRMS:
        {
            /* Compute RMS:
             *   rms = sum(x.^2)/N
             */
            real_T sx2 = 0.0;

            for (i=0; i<N; i++) {
                real_T v = **(uptr++);
                sx2 += v*v;
            }
            *y = sqrt(sx2/N);
            break;
        }
    }
}


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)
{
}


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