%% $RCSfile: sdspacf2.tlc,v $
%% $Revision: 1.4 $  $Date: 2000/06/15 20:37:38 $
%%
%% Copyright 1995-2000 The MathWorks, Inc.
%%
%% Abstract: Target file for the S-Function sdspacf.c
%%           Used for the AutoCorrelation Block

%implements "sdspacf2" "C"

%% Function: BlockInstanceSetup ===============================================
%%
%function BlockInstanceSetup(block, system) void

    %assign INPORT  = 0
    %assign OUTPORT = 0

    %assign cmplx   = LibBlockOutputSignalIsComplex(OUTPORT)
    %assign DType   = cmplx ? "creal_T" : "real_T"

    %assign framebased = LibBlockOutputSignalIsFrameData(OUTPORT)

    %assign numDims = LibBlockInputSignalNumDimensions(INPORT)
    %assign inDims  = LibBlockInputSignalDimensions(INPORT)
    %assign inCols  = (numDims == 2) ? inDims[1] : 1
    %assign inRows  = inDims[0]

    %assign numDims  = LibBlockOutputSignalNumDimensions(OUTPORT)
    %assign outDims  = LibBlockOutputSignalDimensions(OUTPORT)
    %assign outCols  = (numDims == 2) ? outDims[1] : 1
    %assign outRows  = outDims[0]

    %if(!framebased && inRows == 1 && inCols != 1)
        %% Treat single rows as a column vector when not framebased
        %assign inRows  = inCols
        %assign inCols  = 1
        %assign outRows = outCols
        %assign outCols = 1
    %endif

    %assign block = block  + INPORT + OUTPORT + cmplx + DType + ...
                    inCols + inRows + outCols + outRows 


%endfunction %% BlockInstanceSetup


%% Function: Outputs ==========================================================
%%
%function Outputs(block, system) Output
    /* DSP Blockset Autocorrelation (%<ParamSettings.FunctionName>) - %<Name> */
    %%
    %assign bias  = SFcnParamSettings.Bias
    %%
    %% Bias enumerations:
    %assign NO_BIAS  = 1
    %assign BIASED   = 2
    %assign UNBIASED = 3
    %assign COEFF    = 4    
   {
        const int_T inRows = %<inRows>;
        const int_T outRows = %<outRows>;        
        %<DType> *u = %<LibBlockInputSignalAddr(INPORT, "", "", 0)>;
        %<DType> *y = %<LibBlockOutputSignalAddr(OUTPORT,"","",0)>;
        int_T        i;
        %if bias == COEFF
            %% Only need norm variable for COEFF case:
            real_T  norm = 0.0; /* quiet bogus gcc warning */
        %endif
        %%
        %if outCols > 1
            const int_T inCols = %<inCols>;
            const int_T outCols = %<outCols>;
            int_T  ch;

          for(ch=0; ch<outCols; ch++) {
        %endif
            %%
            for (i=0; i<outRows; i++) {
                %%
                %<ACF_sum()> \
                %%
                %switch bias
                    %case NO_BIAS
                    %%                    
                    *y++ = sum;  /* No bias */
                    %break
                    %%
                %case BIASED
                    %%
                    %<ACF_biased()> \
                    %break 
                    %%
                %case UNBIASED
                    %%
                    %<ACF_unbiased()> \
                    %break
                    %%
                %case COEFF
                    %%
                    %<ACF_coeff()> \
                    %%
                %break
                %%
                %endswitch
            }
        %if outCols > 1
            u += inRows;   /* Top of input next channel */
          }
        %endif
     }

%endfunction   %% Outputs


%% Function: ACF_sum ===========================================
%%
%function ACF_sum() Output 
%%
    %<DType> *p0  = u;
    %<DType> *p1  = u+i;
    %%
    %if cmplx
        creal_T sum = {0.0,0.0};
    %else
        real_T  sum  = 0.0;
    %endif
    %%
    int_T  jcnt = inRows-i;

    while(jcnt-- > 0) {
        %if cmplx
            const creal_T *u0 = p0++;
            const creal_T *u1 = p1++;

            /* add: conj(u0) * u1 */
            sum.re += CMULT_XCONJ_RE(*u0, *u1);
            sum.im += CMULT_XCONJ_IM(*u0, *u1);
        %else
            sum += *p0++ * *p1++;
        %endif
    }
%endfunction   %% ACF_sum


%% Function: ACF_biased ===========================================
%%
%function ACF_biased() Output 
%%
    /* Biased ACF estimate: */
    %if cmplx
        y->re     = (1.0 / %<CAST("Real",inRows)>) * sum.re;
        (y++)->im = (1.0 / %<CAST("Real",inRows)>) * sum.im;
    %else
        *y++ = (1.0 / %<CAST("Real",inRows)>) * sum;
    %endif
%%
%endfunction   %% ACF_biased


%% Function: ACF_unbiased ===========================================
%%
%function ACF_unbiased() Output 
%%
    /* Unbiased ACF estimate: */
    %if cmplx
        y->re     = sum.re / (inRows-i);
        (y++)->im = sum.im / (inRows-i);
    %else
        *y++ = sum / (inRows - i);
    %endif
%%
%endfunction   %% ACF_unbiased


%% Function: ACF_coeff ===========================================
%%
%function ACF_coeff() Output 
%%
    /* Autocorrelation coefficient:
    * Set norm = 1/magnitude of zero-lag and use for all remaining lags
    */
    %if cmplx
        if (i == 0) {
        /* In this special case the imaginary part is zero
         * because u0 and u1 are the same point at i=0, 
         * CMULT_XCONJ_IM(*u0, *u1) = 0;
         */
            norm = (sum.re == 0.0) ? 1.0 : 1.0/sum.re;
        }
        y->re     = norm * sum.re;
        (y++)->im = norm * sum.im;
    %else
        if (i == 0) {
	    *y++ = 1.0;      /* Normalize the real element at zero lag to be 1 */
            norm = (sum == 0.0) ? 1.0 : 1.0/sum;
        } else {
	    *y++ = sum * norm;
        }
    %endif
%endfunction   %% ACF_coef



%% [EOF] sdspacf2.tlc