%% $RCSfile: sdspdct3.tlc,v $
%% $Revision: 1.3 $ 
%% $Date: 2000/05/11 18:43:47 $
%%
%% Abstract: TLC Code generation for DSP Blockset DCT block.
%%           Generates in-line and function-call code as needed.
%%
%% Copyright (c) 1995-2000 by The MathWorks, Inc.
%%

%implements sdspdct3 "C"

%include "dsplib.tlc"

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

  %<DSPInitialize(block)>  
  %<DSPAddExternFcnDeclaration("dspfft")>

%endfunction %% BlockInstanceSetup


%% Function: Start ===========================================================
%%
%function Start(block, system) Output
  %%  
  %if !LibIsEqual(TID, "constant")
      %assign xargs = "y_nChans, u_nChans, y_nSamps, u, y"
      %%
      %<GenerateMdlCode(block, system, "Start", xargs)>
      %%
  %endif
  %%
%endfunction %%Start


%% Function: Outputs ==========================================================
%%
%function Outputs(block, system) Output
  %%
  %assign xargs =  "y_nChans, y_nSamps"
  %if SharedIOBuffers
    %assign xargs = xargs + ",u"
  %endif
  %%
  %<GenerateMdlCode(block, system, "Outputs", xargs)>
  %%
%endfunction %% Outputs


%% Function: Terminate ========================================================
%% Abstract:
%%
%function Terminate(block, system) Output
  %%
  %% Remove block instance from database:
  %<DSPTerminate(block)>
  %%
%endfunction %% Terminate


%% Function: CoreStartCode =============================================
%% Abstract:
%%   Generate the core code for the block's operation
%%
%function CoreStartCode(block,system) Output
  %%
    const real_T  den  = sqrt(2*u_nSamps)/2.0;
    const real_T  piN2 = 2.0 * atan(1.0) / u_nSamps;  
    int_T         i;
    
    for (i = 0; i < u_nSamps; i++) {
      dw_Weights[i].re =  cos(i*piN2) / den;
      dw_Weights[i].im = -sin(i*piN2) / den;
    }
    dw_Weights[0].re /= sqrt(2.0);

%endfunction  %% CoreStartCode


%% Function: CoreOutputsCode =============================================
%% Abstract:
%%   Generate the core code for the block's operation
%%
%function CoreOutputsCode(block,system) Output
%%
  if (u_nSamps == 1) {
    %<Outputs_Scalar_Body(block, system)>
  } else {
    %%
    %if !MergeComplexCode(block) && !block.CGInlineFcn 
      %% Complex and Real combined into one function
      if(u_cplx) {
	%<Outputs_Cplx_Body(block, system)>
      } else {
	%<Outputs_Real_Body(block, system)>  
      }
    %else
      %% Complex or Real input
      %if cplx[1] == 1		  
	%<Outputs_Cplx_Body(block, system)>
      %else
	%<Outputs_Real_Body(block, system)>  
      %endif
    %endif
  }

%endfunction %% CoreOutputsCode


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Function: Outputs_Scalar_Body
%function Outputs_Scalar_Body(block, system) Output
  %%
  int_T ch;
  for(ch=0; ch < u_nChans; ch++) {
    *y++ = *u++;
  }
  %% 
%endfunction  %%   Outputs_Scalar_Body
  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Function: Outputs_Cplx_Body
%function Outputs_Cplx_Body(block, system) Output
  %%
  { 
    int_T buff_end = 2 * u_nSamps - 1;
    int_T ch;

    for(ch=0; ch < u_nChans; ch++) {

      int_T         i;
      for (i=0; i<u_nSamps; i++) {
        dw_Buffer[i] = dw_Buffer[buff_end-i] = *u++;
      }
    
      dspfft(2*u_nSamps, dw_Buffer);
    
      /* Multiply by weights: */
      {
        int_T         i;
        for (i=0; i<u_nSamps; i++) {
	  y->re = 0.5*CMULT_RE(dw_Buffer[i], dw_Weights[i]);
	  (y++)->im = 0.5*CMULT_IM(dw_Buffer[i], dw_Weights[i]);
        }
      }    
    } /* End loop over channels */
  }
%endfunction %% Outputs_Cplx_Body


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Function: Outputs_Real_Body
%function Outputs_Real_Body(block, system) Output
  %%
  %assign OUTPORT = 0
  %assign cplx_y  = (LibBlockOutputSignalIsComplex(OUTPORT) != 0)
  %assign dtype   = cplx_y ? "creal_T" : "real_T"

  creal_T *yout = (creal_T *)y;
  int_T    N    = u_nSamps;
  int_T    N2   = N >> 1;
  int_T    N4   = N2 >> 1;
  real_T   N2ele;         /* Store the N2 element */
  int_T    ch;

  for(ch=0; ch<u_nChans; ch++) {
  
      if (N == 2) {
        yout->re = *u++;
        yout->im = *u++;
      
      } else { 
        int_T    i;
    
        for (i=0; i < N4; i++) {
          yout[i].re     = *(u + 4*i);
          yout[i].im     = *(u + 4*i + 2);
	  yout[N4+i].re  = *(u - 4*i + N-1);
	  yout[N4+i].im  = *(u - 4*i + N-1 - 2);
	}
	u += N; /* Increment to next frame */
      }
  
      /* Compute length N/2 FFT: */
      dspfft(N2, yout);    
  
      {
        real_T  theta = -8 * atan(1.0) / N;
        creal_T W     = {1.0, 0.0};
	creal_T twid;
	int_T   i;
    
	/* Complex exponential: cos+jsin */
	twid.re = cos(theta);
	twid.im = sin(theta);
	
	/* Store y[N2] element (always real) */
	N2ele = yout[0].re - yout[0].im;
	
	yout[0].re += yout[0].im; 
	yout[0].im  = 0.0;
	
	for (i = 1; i < N4; i++) {
	  creal_T a, b, c, d;
	  
	  a = yout[i];
	  b = yout[N2-i];
	  
	  c.re = 0.5 * (a.re + b.re);
	  c.im = 0.5 * (a.im - b.im);
	  
	  d.re =  0.5 * (a.im + b.im);
	  d.im = -0.5 * (a.re - b.re);
	  
	  {
	    creal_T ctemp;
	    
	    /* W *= twid */
	    ctemp.re = CMULT_RE(W, twid);
	    ctemp.im = CMULT_IM(W, twid);
	    W = ctemp;
	    
	    /* Precompute W*d: */
	    ctemp.re = CMULT_RE(W, d);
	    ctemp.im = CMULT_IM(W, d);
	    
	    /* y[i] = c + W * d */
	    yout[i].re =  c.re + ctemp.re;
	    yout[i].im =  c.im + ctemp.im;
	    
	    /* y[N2 - i] = conj(c - W*d) */
	    yout[N2 - i].re =  c.re - ctemp.re;
	    yout[N2 - i].im = -c.im + ctemp.im;
	  }
	}
	yout[N4].im = -yout[N4].im;
      }
      
      /* Multiply by weights: */
      {
 	/* Note: There's N2-1 values are conjugates of each other:
	*   for (i = 1; i < N2; i++) {
	*       yr[i + N2] =  yr[N2 - i];
	*       yi[i + N2] = -yi[N2 - i];
	*   }
	*/
	int_T    i;
	
	*y++ = CMULT_RE(yout[0],dw_Weights[0]);
	
	for (i=1; i<N2; i++) {
	  creal_T ctemp  = yout[i];
	  *y++ = CMULT_RE(ctemp,dw_Weights[i]);
	  *y++ = CMULT_XCONJ_RE(ctemp,dw_Weights[N-i]);	    
	}
	y -= N-1;
	
	{
	  real_T *ya = y+2;
	  real_T *yb = y+N-1;
	  int_T  j   = N2-1;
	  
	  while (j-- > 0) {
	    *yb = *ya;
	    for (i=0; i<(yb-ya); i++) {
	      ya[i] = ya[i+1];
	    }
	    ya++;
	    yb--;
	  }
	  y[N2] = N2ele * dw_Weights[N2].re;
	}
      }
      /* Set pointers to next frame. */
       yout += N2;   /* yout is complex so increment by N2 */
       y    += N;    /* y is real so increment by N */
       
    } /* End loop over channels */
    
%endfunction %% Outputs_Real_Body

%% [EOF] sdspdct3.tlc