%% $RCSfile: sdspstdvar2.tlc,v $
%% $Revision: 1.9 $ 
%% $Date: 2000/06/15 20:37:50 $
%%
%% Copyright 1995-2000 The MathWorks, Inc. 
%%
%% Abstract: Find the stdvar value of a vector

%implements sdspstdvar2 "C"


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

  %assign INPORT_DATA  = 0
  %assign INPORT_RESET = 1
  %assign OUTPORT_DATA = 0
  %%
  %assign fcn     = SFcnParamSettings.Fcn
  %assign FcnVar  = 0
  %assign FcnStd  = 1
  %%
  %assign running = CAST("Boolean",(SFcnParamSettings.isRunning == "Yes"))
  %assign cplx    = (LibBlockInputSignalIsComplex(INPORT_DATA) != 0)
  %assign dType   = cplx ? "creal_T" : "real_T"
  %%
  %assign inWidth   = LibBlockInputSignalWidth(INPORT_DATA)
  %assign numDims = LibBlockInputSignalNumDimensions(INPORT_DATA)
  %assign dims    = LibBlockInputSignalDimensions(INPORT_DATA)  
  %%
  %assign isFullMatrix = LibBlockInputSignalIsFullMatrix(INPORT_DATA)
  %assign framebased   = LibBlockInputSignalIsFrameData(INPORT_DATA)
  %assign nChans       = (framebased || isFullMatrix) ? dims[1] : 1
  %assign nSamps       = (framebased || isFullMatrix) ? dims[0] : inWidth
  %assign isFramedRow  = (framebased && (nSamps == 1) && inWidth > 1)

  %assign block = block + ...
                  INPORT_DATA + INPORT_RESET + OUTPORT_DATA + ...
                  cplx        + dType        + inWidth      + ...
                  framebased  + nChans       + nSamps       + ...
                  running     + fcn          + FcnVar       + ...
                  FcnStd      + isFramedRow

%endfunction %% BlockInstanceSetup


%% Function: InitializeConditions =============================================
%%
%function InitializeConditions(block, system) Output
  %%
  %if running
    /* DSP Blockset Stdvar (%<ParamSettings.FunctionName>) - %<Name> */
    %%
    %<LibBlockDWork(Iteration,"", "", 0)> = 0;
    %%
  %endif
%endfunction %% InitializeConditions


%% Function: Outputs ==========================================================
%%
%function Outputs(block, system) Output
  /* DSP Blockset Stdvar (%<ParamSettings.FunctionName>) - %<Name> */
  %%
  %if !running
      %%
      %if inWidth==1
        %%
        %<LibBlockOutputSignal(OUTPORT_DATA, "","", 0)> = 0.0;
                
      %elseif (isFramedRow)
          {
	      int_T j;
	      for (j=0;j < %<nChans>;j++) {
	          %<LibBlockOutputSignal(OUTPORT_DATA,"j","",0)> = 0.0;
	      }		
	  }
      %else	
          %%
          %if !cplx
              %<Stdvar_r(block)>
          %else
              %<Stdvar_c(block)>
          %endif
      %endif %% if inWidth==1  
  %else
    %%
    %% Block could be multirate if and only if it has more than
    %% one input port (eg, a reset input).  If it is multirate,
    %% generate "sample hit" conditional so code executes at the
    %% block's base rate, and not the overall model base rate.
    %%
    %assign differentPortRates = !LibIsSFcnSingleRate(block)
    %%
    %if differentPortRates
      if (%<LibIsSFcnSampleHit("InputPortIdx0")>) { %<LibTaskSFcnComment("InputPortIdx0")>
    %endif
      %%
      {
          typedef struct {
          %<dType> sum; 
          real_T sqsum;
          } sv;

          sv       *sv_data  = (sv *)%<LibBlockDWorkAddr(StdVarData, "", "", 0)>;
          uint32_T *iter_cnt = %<LibBlockDWorkAddr(Iteration, "", "", 0)>;

          %<Running_reset(block)>
          %%
          %if !framebased
              %if !cplx
                  %<Stdvar_run_sr(block)>
              %else
                  %<Stdvar_run_sc(block)>
              %endif
          %else
              %if !cplx
                  %<Stdvar_run_fr(block)>
              %else
                  %<Stdvar_run_fc(block)>
              %endif
          %endif
      }
    %if differentPortRates
      }
    %endif
  %endif

%endfunction


%% Function: Stdvar_r =========================================================
%%
%% Non-running real
%%
%function Stdvar_r(block) Output
%%
  {	 
    real_T *u = %<LibBlockInputSignalAddr(INPORT_DATA,"","",0)>;
    %if nChans > 1
        int j;
        for (j=0; j < %<nChans>; j++) {
    %endif
    %%
    real_T sx  = 0.0;
    real_T sx2 = 0.0;
    int_T  i   = %<nSamps>;

    while(i-- > 0) {
        const real_T val = *u++;
	    sx  += val;
	    sx2 += val * val;
    }
	%if fcn == FcnVar
	    /* Variance */
        %if nChans > 1
            %<LibBlockOutputSignal(OUTPORT_DATA,"j","",0)> = (sx2 - sx*sx / %<nSamps>) / %<nSamps-1>;
            } %% nChans for loop
        %else
	        %<LibBlockOutputSignal(OUTPORT_DATA, "","",0)> = (sx2 - sx*sx / %<nSamps>) / %<nSamps-1>;
        %endif
	%else 
	    /* Standard Deviation */                   
        %if nChans > 1
            %<LibBlockOutputSignal(OUTPORT_DATA,"j","",0)> = sqrt((sx2 - sx*sx / %<nSamps>) / %<nSamps-1>);
            } %% nChans for loop
        %else
	        %<LibBlockOutputSignal(OUTPORT_DATA, "","",0)> = sqrt((sx2 - sx*sx / %<nSamps>) / %<nSamps-1>);
        %endif
    %endif
 }

%endfunction  %% Stdvar_r


%% Function: Stdvar_c =========================================================
%%
%% Non-running complex
%%
%function Stdvar_c(block) Output
%%
  {	 
    creal_T *u = %<LibBlockInputSignalAddr(INPORT_DATA,"","",0)>;
    %if nChans > 1
        int j;
        for (j=0; j < %<nChans>; j++) {
    %endif
    creal_T sx  = {0.0, 0.0};
    real_T  sx2 = 0.0;
    int_T   i   = %<nSamps>;

    while(i-- > 0) {
        sx.re += u->re;
        sx.im += u->im;
        sx2 += CMAGSQ(*u);
        u++;
    }

	%if fcn == FcnVar
	  /* Variance */
        %if nChans > 1
            %<LibBlockOutputSignal(OUTPORT_DATA,"j","",0)> = (sx2 - CMAGSQ(sx) / %<nSamps>) / %<nSamps-1>;
        } %% nChans for loop                    
        %else
	        %<LibBlockOutputSignal(OUTPORT_DATA, "","",0)> = (sx2 - CMAGSQ(sx) / %<nSamps>) / %<nSamps-1>;
        %endif
	%else 
	  /* Standard Deviation */
        %if nChans > 1
            %<LibBlockOutputSignal(OUTPORT_DATA,"j","",0)> = sqrt((sx2 - CMAGSQ(sx) / %<nSamps>) / %<nSamps-1>);
        }  %% nChans for loop                     
        %else
	        %<LibBlockOutputSignal(OUTPORT_DATA, "","",0)> = sqrt((sx2 - CMAGSQ(sx) / %<nSamps>) / %<nSamps-1>);
        %endif
	%endif
    }

%endfunction  %% Stdvar_c


%% Function: Stdvar_run_sr ====================================================
%%
%% Running sample-based real
%%
%function Stdvar_run_sr(block) Output
  %% Real Running
  %%
  if(*iter_cnt == 0) {
      %if inWidth > 1
	  /* Single time step - special case for output value: */
	  real_T *u = %<LibBlockInputSignalAddr(INPORT_DATA,"","",0)>;
	  real_T *y = %<LibBlockOutputSignalAddr(OUTPORT_DATA,"", "", 0)>;
	  int_T i = %<inWidth>;
	  
	  while(i-- > 0) {
		real_T val = *u++;
		sv_data->sum       = val;
		(sv_data++)->sqsum = val * val;
		*y++ = 0.0;
	  }
      %else
          %assign u = LibBlockInputSignal(INPORT_DATA,"","",0)
          %assign y = LibBlockOutputSignal(OUTPORT_DATA,"", "", 0)
	  sv_data->sum   = %<u>;
	  sv_data->sqsum = %<u> * %<u>;
	  %<y> = 0.0;
      %endif
      %%
      (*iter_cnt)++;

  } else {

      (*iter_cnt)++;

      /* Protect against divide by zeros after increment. */
      if (*iter_cnt == 0) {
	  real_T *u = %<LibBlockInputSignalAddr(INPORT_DATA,"","",0)>;
          int_T i;
          *iter_cnt = 1;

          /* Clear states when counter rolls over */
	  for(i=0; i<%<inWidth>; i++) {
		real_T val         = u[i];
		sv_data->sum       = val;
		(sv_data++)->sqsum = val * val;
	  }
      }

      {
	%assign yout = (inWidth > 1) ? "*y++" : "*y"
	%assign uptr = (inWidth > 1) ? "*u++" : "*u"
        %%
	real_T *u     = %<LibBlockInputSignalAddr(INPORT_DATA,"","", 0)>;
	real_T *y     = %<LibBlockOutputSignalAddr(OUTPORT_DATA,"", "", 0)>;
	uint32_T iter = *iter_cnt;	    
        %%
	%if inWidth > 1
	    int_T   i;
	  
            for(i=0; i<%<inWidth>; i++) {
	%endif
	        %%
	        /* Update the state: */
	        real_T val = %<uptr>;
	        sv_data->sum += val;
	        sv_data->sqsum += val * val;
                %%
	        %if fcn == FcnVar
		    /* Variance */
		    %<yout> = (sv_data->sqsum - sv_data->sum * sv_data->sum / iter) / (iter-1);
	        %else 
		    /* Standard Deviation */
		    %<yout> = sqrt((sv_data->sqsum - sv_data->sum * sv_data->sum / iter) / (iter-1));
	        %endif
	        %%
	        %if inWidth > 1
		    sv_data++;
	        %endif
	        %%
	%if inWidth > 1
	    }
	%endif
      }
  }  /* if(*iter_cnt == 0) */

%endfunction  %% Stdvar_run_sr


%% Function: Stdvar_run_sc ====================================================
%%
%% Running sample-based complex
%%
%function Stdvar_run_sc(block) Output
%% COMPLEX Running
%%
  if(*iter_cnt == 0) {
	%if inWidth > 1
	    /* Single time step - special case for output value: */
	    creal_T *u = %<LibBlockInputSignalAddr(INPORT_DATA,"","",0)>;
	    real_T  *y = %<LibBlockOutputSignalAddr(OUTPORT_DATA,"", "", 0)>;
	    int_T    i = %<inWidth>;
      
	    while(i-- > 0) {
		    creal_T val = *u++;
		    sv_data->sum       = val;
		    (sv_data++)->sqsum = CMAGSQ(val);
		    *y++ = 0.0;
	    }
	%else
	    %assign u = LibBlockInputSignal(INPORT_DATA,"","",0)
	    %assign y = LibBlockOutputSignal(OUTPORT_DATA,"", "", 0)
            %%
	    sv_data->sum   = %<u>;
	    sv_data->sqsum = CMAGSQ(%<u>);
	    %<y> = 0.0;
	%endif
        %%
	(*iter_cnt)++;

  } else {

      (*iter_cnt)++;

      /* Protect against divide by zeros after increment. */
      if (*iter_cnt == 0) {
	    creal_T *u = %<LibBlockInputSignalAddr(INPORT_DATA,"","",0)>;
            int_T i;
            *iter_cnt = 1;

            /* Clear the states when the counter rolls over */
	    for(i=0; i<%<inWidth>; i++) {
		    creal_T val        = u[i];
		    sv_data->sum       = val;
		    (sv_data++)->sqsum = CMAGSQ(val);
	    }
      }

      {
	%assign yout  = (inWidth !=1) ? "*y++" : "*y"
	%assign uptr  = (inWidth !=1) ? "*u++" : "*u"
        %%
	creal_T *u    = %<LibBlockInputSignalAddr(INPORT_DATA,"","", 0)>;
	real_T  *y    = %<LibBlockOutputSignalAddr(OUTPORT_DATA,"", "", 0)>;
	uint32_T iter = *iter_cnt;	    
	%if inWidth > 1
	    int_T   i;
	  
            for(i=0; i<%<inWidth>; i++) {
	%endif
	        %%
	        /* Update the state: */
	        creal_T val = %<uptr>;
	        sv_data->sum.re += val.re;
	        sv_data->sum.im += val.im;
	        sv_data->sqsum  += CMAGSQ(val);

	        %if fcn == FcnVar
		    /* Variance */
		    %<yout> = (sv_data->sqsum - CMAGSQ(sv_data->sum) / iter) / (iter-1);
	        %else 
		    /* Standard Deviation */
		    %<yout> = sqrt((sv_data->sqsum - CMAGSQ(sv_data->sum) / iter) / (iter-1));
	        %endif
	        %%
	        %if inWidth > 1
		    sv_data++;
	        %endif	      
	        %%
	%if inWidth > 1
	    }   
	%endif
      }
  }

%endfunction  %% Stdvar_run_sc


%% Function: Stdvar_run_fr ====================================================
%%
%% Running frame-based real
%%
%function Stdvar_run_fr(block) Output
{
    real_T   *u = %<LibBlockInputSignalAddr(INPORT_DATA,"","",0)>;
    real_T   *y = %<LibBlockOutputSignalAddr(OUTPORT_DATA,"", "", 0)>;
    uint32_T  iter;
    int_T     i;

    /* Clear the states */
    if (*iter_cnt == 0) {
        for(i=0; i<%<nChans>; i++) {
            sv_data->sum       = 0.0;
            (sv_data++)->sqsum = 0.0;        
        }
        sv_data -= %<nChans>;
    }
        
    for(i=0; i<%<nChans>; i++) {
        real_T sx  = sv_data->sum;
        real_T sx2 = sv_data->sqsum;
        int_T  j;

        iter = *iter_cnt;

        for(j=0; j<%<nSamps>; j++) {
            const real_T val = *u++;
            sx  += val;
            sx2 += val * val;

            iter++;

            /* Protect against divide by zeros after increment. */
            if (iter == 0) {
                iter = 1;

                /* Clear the states */ 
                sx  = 0.0;
                sx2 = 0.0;
            }

            if(iter==1) {
                *y++ = 0.0;
            } else {
                const real_T yp = (sx2 - sx*sx / iter) / (iter - 1);
                %%
              	%if fcn == FcnVar
                    /* Variance */
                    *y++ = yp;
                %else
                    /* Standard Deviation */  
                    *y++ = sqrt(yp);
                %endif
            }
        }
        
        /* Update state: */
        sv_data->sum       = sx;
        (sv_data++)->sqsum = sx2;
    }
    *iter_cnt = iter;  /* Update iteration count DWork */
}

%endfunction  %% Stdvar_run_fr


%% Function: Stdvar_run_fc ====================================================
%%
%% Running frame-based complex
%%
%function Stdvar_run_fc(block) Output
{
    creal_T  *u = %<LibBlockInputSignalAddr(INPORT_DATA,"","",0)>;
    real_T   *y = %<LibBlockOutputSignalAddr(OUTPORT_DATA,"", "", 0)>;
    uint32_T  iter;
    int_T     i;        

    /* Clear the states */
    if (*iter_cnt == 0) {
        for(i=0; i<%<nChans>; i++) {
            sv_data->sum.re    = 0.0;
            sv_data->sum.im    = 0.0;
            (sv_data++)->sqsum = 0.0;
        }
        sv_data -= %<nChans>;
    }
    
    for(i=0; i<%<nChans>; i++) {
        creal_T sx  = sv_data->sum;
        real_T  sx2 = sv_data->sqsum;
        int_T   j;

        iter = *iter_cnt;
        
        for(j=0; j<%<nSamps>; j++) {
            const creal_T val = *u;
            sx.re += val.re;
            sx.im += val.im;
            sx2 += CMAGSQ(*u);
            u++;
            iter++;

            /* Protect against divide by zeros after increment. */
            if (iter == 0) {
                iter = 1;

                /* Clear the states */
                sx.re = 0.0;
                sx.im = 0.0;
                sx2   = 0.0;
            }

            if(iter==1) {
                *y++ = 0.0;
            } else {
                const real_T yp = (sx2 - CMAGSQ(sx) / iter) / (iter - 1);
                %%
                %if fcn == FcnVar
                /* Variance */
                    *y++ = yp;
                %else
                    /* Standard Deviation */  
                    *y++ = sqrt(yp);
                %endif
            }
        }
        
        /* Update state: */
        sv_data->sum       = sx;
        (sv_data++)->sqsum = sx2;
    }
    *iter_cnt = iter;  /* Update iteration count DWork */
}

%endfunction  %% Stdvar_run_fc



%% Function: Running_reset ====================================================
%%
%% Running reset code
%%
%function Running_reset(block) Output
%%
    %assign resetport = CAST("Boolean",(SFcnParamSettings.isRstPort == "Yes"))
    %%
    %if resetport
        if(%<LibBlockInputSignal(INPORT_RESET,"","",0)> != 0.0) { /* Reset iteration count. */
	    *iter_cnt = 0;
	}
    %endif
%%
%endfunction  %% Running_reset


%% [EOF] sdspstdvar2.tlc