%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ModelRun.m    : radix 2 fft model run
% author        : duweitao
% date          : 2008.May.25
% version       : 0.1
% modify        : initial version
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function simData = ModelRun(cfg, simData)
  model1Data = Model1Run(cfg, simData);
  model2Data = Model2Run(cfg, simData);
  
  simData.model1Data = model1Data;
  simData.model2Data = model2Data;
end % ModelRun()
% -----------------------------------------------------
% function Model1Run(), this model only verifying
% the basic FFT concepts, not including the  bitwise 
% operation and TWIDDLE FACTORs symmetrical and look-up
% table for output re-mapping
% -----------------------------------------------------
function model1Data = Model1Run(cfg, simData)
  fprintf(1, '# Model1Run() \n');
% load from struct obj
  random_T = cfg.random_T   ;  
  sine_T   = cfg.sine_T     ;  
  stimType = cfg.stimType   ;  
  freqNorm = cfg.freqNorm   ;  
  Amp      = cfg.Amp        ;  
  N        = cfg.N          ;  

  inputData = simData.inputData;
  
  nLv    = log2(N)         ;% number of levels
  nBfL   = N/2             ;% number of butterflies
  idxLv  = [0:nLv-1]       ;% index of level
  idxBfL = [0:nBfL-1]      ;% index of butterfly
  nBfG   = 2.^(nLv-idxLv-1);% number of butterflies in one group
  nGpL   = 2.^idxLv        ;% number of groups in one level

  % row index is idxBfL+1 col index is idxLv + 1 and element value
  % is idxGpL value or idxBfG value
  idxGpL = zeros(nBfL,nLv) ;
  idxBfG = zeros(nBfL,nLv) ;
  idxWn  = zeros(nBfL,nLv) ; % index of twiddle factor
  addrRd0 = zeros(nBfL, nLv);
  addrRd1 = zeros(nBfL, nLv);
% 2D data buffer help debugging   
  dataBuf = zeros(N, nLv+1);    
  dataBuf(:,1) = inputData;
% twiddle factor table
  Wn = zeros(nBfL,1);

  for(idx1 = 1: nLv) 
    idxGpL(:, idx1) = floor(idxBfL ./ nBfG(idx1)) .' ;
    idxBfG(:, idx1) = mod(idxBfL, nBfG(idx1)) .' ;
  end

  % get the read address of every butterfly 
  for(idx1 = 1:nLv)     % by level
    for(idx2 = 1:nBfL)  % by butterfly index
      addrRd0(idx2, idx1) = idxBfG(idx2, idx1) + idxGpL(idx2, idx1)*nBfG(idx1)*2; 
      addrRd1(idx2, idx1) = addrRd0(idx2, idx1) + nBfG(idx1); 
      idxWn(idx2, idx1)   = idxBfG(idx2, idx1)*(2^(idx1-1));
    end
  end
  % get twiddle factors
  for(idx2 = 1:nBfL)  % total nBfL twiddle factors
    theta = 2*pi*(idx2-1)/N;
    Wn(idx2) = exp(-j*theta);
  end

  % get the bit-reorder map
  outMap = zeros(N,1);
  for(idx1 = 0:N-1)
    for(idx2 = 0:nLv-1)
      bitV = bitget(idx1, idx2 +1);
      outMap(idx1 +1) = bitset(outMap(idx1 +1), nLv-1-idx2 +1,  bitV);
    end
  end

  % bufferfly 
  for(idx1 = 1:nLv)     % by level
    for(idx2 = 1:nBfL)  % by butterfly index
      % get butterfly data
      bfIn0  = dataBuf(addrRd0(idx2, idx1)+1, idx1);
      bfIn1  = dataBuf(addrRd1(idx2, idx1)+1, idx1);
      bfTW   = Wn(idxWn(idx2, idx1)+1);
      % computing
      bfOut0 = bfIn0 + bfIn1;
      bfOut1 = (bfIn0 - bfIn1)*bfTW;
      % write out
      dataBuf(addrRd0(idx2, idx1)+1, idx1+1) = bfOut0; 
      dataBuf(addrRd1(idx2, idx1)+1, idx1+1) = bfOut1; 
    end
  end

  % reorder the output
  output = zeros(N,1);
  for(idx1 = 0:N-1) 
    output(outMap(idx1 +1) +1) = dataBuf(idx1 +1, nLv +1);
  end

  model1Data.nLv    = nLv    ;
  model1Data.nBfL   = nBfL   ;
  model1Data.idxLv  = idxLv  ;
  model1Data.idxBfL = idxBfL ;
  model1Data.nBfG   = nBfG   ;
  model1Data.nGpL   = nGpL   ;
  model1Data.idxGpL = idxGpL ;
  model1Data.idxBfG = idxBfG ;

  model1Data.idxWn   = idxWn   ;         
  model1Data.addrRd0 = addrRd0 ;     
  model1Data.addrRd1 = addrRd1 ;     
  model1Data.dataBuf = dataBuf ;  
  model1Data.Wn      = Wn      ;   
  model1Data.mapData = outMap  ;
  model1Data.output  = output  ;

  matOut(:,1) = fft(inputData);
  diff = matOut - output;
  % check the average error is OK
  if(sum(abs(diff))/N < 1E-10)
    fprintf(1, '# model1 run OK\n');
  else
    fprintf(1, '# model1 run ERROR\n');
  end
end % Model1Run
%===============================================================================

% -----------------------------------------------------
% function Model2Run(), this model using the  bitwise 
% operation and TWIDDLE FACTORs symmetrical and look-up
% table for output re-mapping
% -----------------------------------------------------
function model2Data = Model2Run(cfg, simData)
  fprintf(1, '# Model2Run() \n');
% load from struct obj
  random_T = cfg.random_T   ;  
  sine_T   = cfg.sine_T     ;  
  stimType = cfg.stimType   ;  
  freqNorm = cfg.freqNorm   ;  
  Amp      = cfg.Amp        ;  
  N        = cfg.N          ;  

  inputData = simData.inputData;
  
  nLv    = log2(N)         ;% number of levels
  nBfL   = N/2             ;% number of butterflies
  idxLv  = [0:nLv-1]       ;% index of level
  idxBfL = [0:nBfL-1]      ;% index of butterfly
  nBfG   = 2.^(nLv-idxLv-1);% number of butterflies in one group
  nGpL   = 2.^idxLv        ;% number of groups in one level

  % row index is idxBfL+1 col index is idxLv + 1 and element value
  % is idxGpL value or idxBfG value
  idxGpL = zeros(nBfL,nLv) ;
  idxBfG = zeros(nBfL,nLv) ;
  idxWn  = zeros(nBfL,nLv) ; % index of twiddle factor
  addrRd0 = zeros(nBfL, nLv);
  addrRd1 = zeros(nBfL, nLv);
% 2D data buffer help debugging   
  dataBuf = zeros(N, nLv+1);    
  dataBuf(:,1) = inputData;
% twiddle factor table
  Wn = zeros(nBfL,1);
  WnTable41 = zeros(N/4,1);   % this is synmetric twiddle table


  for(idx1 = 1: nLv) 
    idxGpL(:, idx1) = floor(idxBfL ./ nBfG(idx1)) .' ;
    idxBfG(:, idx1) = mod(idxBfL, nBfG(idx1)) .' ;
  end

  % get the read address of every butterfly 
  for(idx1 = 1:nLv)     % by level
    for(idx2 = 1:nBfL)  % by butterfly index
      addrRd0(idx2, idx1) = idxBfG(idx2, idx1) + idxGpL(idx2, idx1)*nBfG(idx1)*2; 
      addrRd1(idx2, idx1) = addrRd0(idx2, idx1) + nBfG(idx1); 
      idxWn(idx2, idx1)   = idxBfG(idx2, idx1)*(2^(idx1-1));
    end
  end
  % get twiddle factors
  % get symmetrical wave table
  for(idx1 = 0:N/4-1)
    WnTable41(idx1 +1) = -sin(2*pi/N*idx1);
  end
  % image part of Wn
  for(idx1 = 0:N/2-1)
    if(idx1 < N/4)
      ImWn = WnTable41(idx1 +1);
    else
      if(idx1 == N/4)
        ImWn = -1;
      else
        idxMap = N/2 - idx1;    % map to imag 0~N/4-1
        ImWn = WnTable41(idxMap+1);
      end % else(idx1)
    end % else(idx1)
    Wn(idx1+1) = ImWn*j;
  end
  % real part of Wn
  for(idx1 = 0:N/2-1)
    if(idx1 < N/4)
      idxMapIm = N/4 - idx1;
      ReWn = -imag(Wn(idxMapIm+1));
    else
      idxMapIm = mod(idx1,N/4);    % map to Imag 0~N/4-1
      ReWn = imag(Wn(idxMapIm+1));
    end % else(idx1)
    Wn(idx1+1) = ReWn + Wn(idx1+1);
  end

  % get the bit-reorder map
  idxMap8bits = zeros(256,1);
  for(idx1 = 0:255)
    for(idx2 = 0:7)
      bitV = bitget(idx1, idx2 +1);
      idxMap8bits(idx1 +1) = bitset(idxMap8bits(idx1 +1), 8-idx2 ,  bitV);
    end
  end

  outMap = zeros(N,1);
  for(idx1 = 0:N-1)
    idxHigh = bitand(idx1, hex2dec('ff00'));
    idxHigh = bitshift(idxHigh, -8);
    idxLow  = bitand(idx1, hex2dec('00ff'));
    revHigh = idxMap8bits(idxHigh+1);
    revLow  = idxMap8bits(idxLow+1);
    revIdx = bitor(bitshift(revLow,8), revHigh);
    revIdx = bitshift(revIdx, nLv-16);
    outMap(idx1+1) = revIdx;
  end

  % bufferfly 
  for(idx1 = 1:nLv)     % by level
    for(idx2 = 1:nBfL)  % by butterfly index
      % get butterfly data
      bfIn0  = dataBuf(addrRd0(idx2, idx1)+1, idx1);
      bfIn1  = dataBuf(addrRd1(idx2, idx1)+1, idx1);
      bfTW   = Wn(idxWn(idx2, idx1)+1);
      % computing
      bfOut0 = bfIn0 + bfIn1;
      bfOut1 = (bfIn0 - bfIn1)*bfTW;
      % write out
      dataBuf(addrRd0(idx2, idx1)+1, idx1+1) = bfOut0; 
      dataBuf(addrRd1(idx2, idx1)+1, idx1+1) = bfOut1; 
    end
  end

  % reorder the output
  output = zeros(N,1);
  for(idx1 = 0:N-1) 
    output(outMap(idx1 +1) +1) = dataBuf(idx1 +1, nLv +1);
  end

  model2Data.nLv    = nLv    ;
  model2Data.nBfL   = nBfL   ;
  model2Data.idxLv  = idxLv  ;
  model2Data.idxBfL = idxBfL ;
  model2Data.nBfG   = nBfG   ;
  model2Data.nGpL   = nGpL   ;
  model2Data.idxGpL = idxGpL ;
  model2Data.idxBfG = idxBfG ;

  model2Data.idxWn   = idxWn   ;         
  model2Data.addrRd0 = addrRd0 ;     
  model2Data.addrRd1 = addrRd1 ;     
  model2Data.dataBuf = dataBuf ;  
  model2Data.WnTable41 = WnTable41;
  model2Data.Wn      = Wn      ;   
  model2Data.mapData = outMap  ;
  model2Data.output  = output  ;

  matOut(:,1) = fft(inputData);
  diff = matOut - output;
  % check the average error is OK
  if(sum(abs(diff))/N < 1E-10)
    fprintf(1, '# model2 run OK\n');
  else
    fprintf(1, '# model2 run ERROR\n');
  end
end % Model2Run
%===============================================================================
  
