%% $RCSfile: sdspfft2.tlc,v $
%% $Revision: 1.2 $   
%% $Date: 2000/06/05 18:26:54 $
%%
%% Abstract: TLC Code generation for DSP Blockset FFT block.
%%           Generates in-line and function-call code as needed.
%%
%% Copyright (c) 1995-2000 by The MathWorks, Inc.
%%

%implements sdspfft2 "C"

%include "dsplib.tlc"

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

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


%% Function: Outputs ==========================================================
%%
%function Outputs(block, system) Output
  %%
  %% The list of variables to exclude must be comma separated and 
  %% stored in a single long string.  The variable names must match
  %% the convention of automatically generated argument names.
  %% Names are case sensative. 
  %%
  %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: CoreOutputsCode ====================================================
%% Abstract:
%%   Generate the core code for the block's operation
%%
%function CoreOutputsCode(block,system) Output
%%
%if !MergeComplexCode(block) && !block.CGInlineFcn 
    %% Complex and Real combined into one function
    if(u_cplx) {
      /* cast necessary because we are passed a void ptr */
      creal_T *u_tmp = (creal_T *)u;
      creal_T *y_tmp = (creal_T *)y;
      creal_T *u     = u_tmp;
      creal_T *y     = y_tmp;
      %<Outputs_Cplx_Body(block, system)>
    } else {
      real_T  *u_tmp = (real_T *)u;
      creal_T *y_tmp = (creal_T *)y;
      real_T  *u     = u_tmp;
      creal_T *y     = y_tmp;
      %<Outputs_Real_Body(block, system)>  
    }
%else
    %% Complex or Real input (whic is cplx[1])
    %% The output complexity is stored in cplx[0]
    %if cplx[1]	== 1	  
      %<Outputs_Cplx_Body(block, system)>
    %else
      %<Outputs_Real_Body(block, system)>  
    %endif
%endif

%endfunction %% CoreOutputsCode


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Function: Outputs_Cplx_Body
%function Outputs_Cplx_Body(block, system) Output
  %%
  {
    int_T f;
    for (f = 0; f++ < u_nChans; ) {
      %% 
      %%  Don't copy input to outputs if they share same buffer.
      %if !SharedIOBuffers
	int_T i;
	/* Copy inputs into outputs, since FFT operates in place: */
	for (i=0; i++ < u_nSamps;) {
	  *y++ = *u++;
	}
	y -= u_nSamps;
      %else
	/* Input and Output share same buffers */
      %endif
      dspfft(u_nSamps, y);
      y += u_nSamps;
    }
  }
  
%endfunction %% Outputs_Cplx_Body


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Function: Outputs_Real_Body
%function Outputs_Real_Body(block, system) Output
  %%
  {
    int_T f;
    for (f = 0; f++ < u_nChans; ) {
      
      /* Scalar */
      if (u_nSamps == 1) {
	y->re = *u++;
	(y++)->im = 0.0;
	
      } else {
	const int_T N2 = u_nSamps >> 1;  /* Length of complex FFT is half width */
	const int_T N4 = N2 >> 1;       /* Half length of complex FFT  */

	/* Copy inputs to outputs */
	/* A length-N real input is interpreted as a length-N/2 complex input*/
	int_T i;
	
	for (i=N2; i-- > 0;) {
	  /* Careful! The re and im parts in uptr may not be contiguous! */
	  y->re     = *u++;
	  (y++)->im = *u++;
	}
	y -= N2;
      %%
      %if SharedIOBuffers
	%% 
	%% Real input will NEVER share the same buffer as its output because 
	%% the output will be complex and thus different from the input type.
	%% However, we can get to this spot if the user selects to merge
	%% complex and real code into one function.
	%%
	%warning Real inputs cannot share the same buffer as Complex outputs
	/* Real inputs cannot share the same buffer as Complex outputs */
	%%
      %endif
      %%
	/* Compute length N/2 FFT: */
	dspfft(N2, y);
	
	{
	  creal_T      W     = {1.0, 0.0};
	  const real_T theta = -8 * atan(1.0) / u_nSamps;
	  creal_T      twid;
	  int_T        i;
	  
	  /* Complex exponential: cos+jsin */
	  twid.re = cos(theta);
	  twid.im = sin(theta);
	  
	  y[N2].re = y[0].re - y[0].im;
	  y[N2].im = 0.0;
	  
	  y[0].re += y[0].im; 
	  y[0].im  = 0.0;
	  
	  for (i = 1; i < N4; i++) {
	    creal_T a, b, c, d;
	    
	    a = y[i];
	    b = y[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 * (b.re - a.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 */
	      y[i].re = c.re + ctemp.re;
	      y[i].im = c.im + ctemp.im;
	      
	      /* y[N2 - i] = conj(c - W*d) */
	      y[N2 - i].re =  c.re - ctemp.re;
	      y[N2 - i].im = -c.im + ctemp.im;
	    }
	  }            
	  y[N4].im = -y[N4].im;
	  
	  /* y[i+N2] = conj(y[N2-i]) */
	  for (i = 1; i < N2; i++) {
	    y[i + N2].re =  y[N2 - i].re;
	    y[i + N2].im = -y[N2 - i].im;
	  }
	}
	/* Make sure uptr and y are set up for next channel */
	y += u_nSamps;
	
      } /* end of  N > 1 case */
    }
    
  }
%endfunction %% Outputs_Real_Body

%% [EOF] sdspfft.tlc