% ////////////////////////////////////////////////////////////////////
% FILE,     HBF2LInit.m
% PROJECT,  HBF2L
% AUTHOR,   weitao.du
% DATE,     2007.Jun.11th
% VERSION,  0.2
% VER MOIDFICATION
% 0.2 - 0.1, modify the coefficent scale, now all coeffs sum should be
% 1, so the amplify is 1
% ////////////////////////////////////////////////////////////////////
% DECRIPTION
% 2 Level Half Band Fir Filter Design and Simulation, data initial
% ////////////////////////////////////////////////////////////////////

function simData = HBF2LInit(cfg);
  fprintf(1, 'HBF2LInit(), run\n');
  simData.NOUSE = 'NOUSE';
  simData = FilterInit(cfg, simData);
  simData = GenStim(cfg, simData);
  

end  % function HBF2LInit();

% ////////////////////////////////////////////////////////////////////
% FilterInit(), Init the filter data
% ////////////////////////////////////////////////////////////////////

function simData = FilterInit(cfg, simData)

  bwNoise   = cfg.bwNoise      ;
  WGNoise   = cfg.WGNoise      ;
  sine      = cfg.sine         ;
  stepUp    = cfg.stepUp       ;
  stimType  = cfg.stimType     ;
  stimLen   = cfg.stimLen      ;

  WLCfg = cfg.WLCfg;
                                
  inputWL   = WLCfg.inputWL    ;
  coeffWL   = WLCfg.coeffWL    ;
  outputWL  = WLCfg.outputWL   ;
                                
  K1        = cfg.K1           ;
  K2        = cfg.K2           ;
             
  passBand1 = cfg.passBand1    ;
  passBand2 = cfg.passBand2    ;

  tapHbf1   = 4*K1 -2          ; % for Matlab 'firhalfband()'
  tapHbf2   = 4*K2 -2          ;              

  % get coefficents 
  coeff1AllFlt = firhalfband(tapHbf1, passBand1) .';
  coeff2AllFlt = firhalfband(tapHbf2, passBand2) .';
  % the max Val of Coeff should be 0.5
  if(max(coeff1AllFlt) ~= 0.5)
    fprintf(1, 'WARNING, max val of coeff1AllFlt should be 0.5');
  end
  if(max(coeff2AllFlt) ~= 0.5)
    fprintf(1, 'WARNING, max val of coeff2AllFlt should be 0.5');
  end
  coeff1AllFlt = coeff1AllFlt*2;  
  coeff2AllFlt = coeff2AllFlt*2;
  
  % all taps Fixed Point Coeff
  maxCoeffVal  = 2^(coeffWL-1) - 1;     % coeffWL include the sign bit 
  coeff1AllFxp = floor(coeff1AllFlt*maxCoeffVal+0.5);
  coeff2AllFxp = floor(coeff2AllFlt*maxCoeffVal+0.5);
  
  % the condense coefficents, not include 0 and max value
  coeff1CdsFxp = coeff1AllFxp(1:2:2*K1);
  coeff2CdsFxp = coeff2AllFxp(1:2:2*K2);
  coeff1CdsFlt = coeff1AllFlt(1:2:2*K1);
  coeff2CdsFlt = coeff2AllFlt(1:2:2*K2);


  hHfbFxp12 = fvtool(coeff1AllFxp,1, coeff2AllFxp, 1,'Grid', 'on', 'Legend', 'on');
  legend(hHfbFxp12, 'LEVEL1 HBF', 'LEVEL2 HBF');

  inLv2Len = 2*(4*K1-1);
  inLv2 = zeros(inLv2Len,1);
  inLv2(2:2:inLv2Len) = coeff1AllFxp;
  % step response of total HBF
  out = conv(inLv2, coeff2AllFxp);
  hHfbFxpTOAL= fvtool(out, 1,'Grid', 'on', 'Legend', 'on');
  legend(hHfbFxpTOAL, 'TOTAL 2 LEVEL HBF');

  filterInit.coeff1AllFlt =  coeff1AllFlt ;
  filterInit.coeff2AllFlt =  coeff2AllFlt ;
  filterInit.coeff1AllFxp =  coeff1AllFxp ;
  filterInit.coeff2AllFxp =  coeff2AllFxp ;
  filterInit.coeff1CdsFlt =  coeff1CdsFlt ; 
  filterInit.coeff2CdsFlt =  coeff2CdsFlt ;  
% filterInit.hHfbFxp      =  hHfbFxp      ;

  simData.filterInit      =  filterInit   ;

end  % fucntion FiterInit()

% ////////////////////////////////////////////////////////////////////
% GenStim(), generate the stimulation data
% ////////////////////////////////////////////////////////////////////

function simData = GenStim(cfg, simData)
  inputBw       = cfg.inputBw      ;
  inputFs       = cfg.inputFs      ;
  WGNoise       = cfg.WGNoise      ;
  bwNoise       = cfg.bwNoise      ;
  sine          = cfg.sine         ;
  sineAry       = cfg.sineAry      ;
  stepUp        = cfg.stepUp       ;
  stimType      = cfg.stimType     ;
  Vp            = cfg.Vp           ;        
  sineFreq      = cfg.sineFreq     ;
  sineDc        = cfg.sineDc       ;

  baseAryFreq   = cfg.baseAryFreq  ;
  sineAryNum    = cfg.sineAryNum   ;
  stimLen       = cfg.stimLen      ;
  loadDbgStim   = cfg.loadDbgStim  ;
  DBGFILENAME   = cfg.DBGFILENAME  ;

  WLCfg     = cfg.WLCfg;
  
  inputWL   = WLCfg.inputWL    ;
  coeffWL   = WLCfg.coeffWL    ;
  outputWL  = WLCfg.outputWL   ;

  K1        = cfg.K1           ;
  K2        = cfg.K2           ;
             
  passBand1 = cfg.passBand1    ;
  passBand2 = cfg.passBand2    ;

  inputData = 0;
  if(loadDbgStim) 
    fprintf(1, 'HBF2LInit(), DBGUG MODE load Debug file "%s"\n', DBGFILENAME);
    if(~exist(DBGFILENAME));
      fprintf(1, 'HBF2LInit(), FILE "%s" not exist\n', DBGFILENAME);
      return;
    end
    loadData    = load(DBGFILENAME);
    inputData   = loadData.inputData;
  else
    if(stimType == bwNoise) 
      inputData = BandNoise(stimLen, inputBw, inputFs );
      inputData = inputData * Vp;
    end
    if(stimType == WGNoise) 

       dBw          = 1;
       inputData    = wgn(stimLen, 1, dBw);
       maxOut       = max(abs(inputData));
       inputData    = inputData/(maxOut);
       inputData    = inputData * Vp;
    end
    if(stimType == sine) 
      inputData = 0 ;
      dc = sineDc;
      f0 = sineFreq; 
      inputData = sin(2*pi*f0*[0:stimLen-1]/inputFs);
      inputData = inputData * Vp;
      inputData = inputData + dc;
    end
    if(stimType == sineAry) 
      inputData = 0 ;
      dc = sineDc;
      for(idx = 1:sineAryNum)
        inputData = inputData + sin(2*pi*baseAryFreq*idx*[0:stimLen-1]/inputFs);
      end
      inputData = inputData/max(abs(inputData));
      inputData = inputData * Vp;
      inputData = inputData + dc;
    end
    if(stimType == stepUp)
      M     = -4; 
      inputData = mod([M:M:stimLen*M], K1*M) / (2^(inputWL-1));
    end

  end % else (loadDbgStim) 

  simData.inputData = inputData;

end % function GenStim()

% ////////////////////////////////////////////////////////////////////
% function BandNoise(), generate band noise signal, that 
% amplitude rescale to [-1, 1] 
% ////////////////////////////////////////////////////////////////////

function out = BandNoise(len, bw, fs)
  dBw = 1;
  out = wgn(len, 1, dBw);
  
  passBw    = bw;
  stopBw    = bw + bw/10;
  % use kasier window design fir
  fcuts = [passBw, stopBw];
  % gain of pass band, stop band
  mags = [60  0];
  % ripple of pass band, stop band
  devs = [1E-2 2E-2];
  [n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,fs);
  % hh is fir coefficents  
  hh = fir1(n,Wn,ftype,kaiser(n+1,beta));
  
  out = conv(out, hh);
  out = out(1:len).';
  
  maxOut = max(abs(out));
  out = out/(maxOut);
  

end
