function [stateLevels, ...
          readoutTransitions, ...
          readoutAberrations, ...
          readoutCycles, ...
          measureLines] = BilevelMetrics(x, t, ...
                                         enableAutoLevel, ...
                                         stateLevels, ...
                                         pctTol, ...
                                         pctRefLevels, ...
                                         settlingSeekDuration, ...
                                         xLim, ...
                                         enableAberrations, ...
                                         enableCycles, ...
                                         enableMeasureLines, ...
                                         numMeasLines) 
% Time Scope Bilevel Measurements callback routine
%     x                     - input voltage samples
%     t                     - sample times
%     enableAutoLevel       - true if auto-leveling is desired
%     stateLevels           - low and high logic levels
%     pctTol                - percent state level tolerance (usu 2%)
%     pctRefLevels          - percent reference levels (usu 10% 50% 90%)
%     settlingSeekDuration  - duration over which to compute settling
%     xLim                  - Time Scope axes borders
%     enableAberrations     - enables overshoot/undershoot/settling
%     enableCycles          - enables pulse width, period, etc.
%     enableMeasureLines    - enables annotation (for use in debug)
%     numMeasLines          - expected number of output lines
%
%  This file is for internal use only.

%   Copyright 2013 The MathWorks, Inc.
%#codegen
  
  readoutTransitions = NaN(1, 9);
  readoutAberrations = NaN(1, 8);
  readoutCycles = NaN(1, 8);
  coder.varsize('measureLines.xData');
  coder.varsize('measureLines.yData');
  measureLines = repmat(struct('xData', [], 'yData', []), 1, numMeasLines);
  
  if enableAutoLevel
    stateLevels = computeStateLevels(x);
  end
    
  readoutTransitions(1, [2 1]) = stateLevels;
  
  sl = stateLevels;
  pslt = pctTol;
  prl = pctRefLevels;
  A = diff(sl) / 100;
  us = sl(2);               % upper state
  ls = sl(1);               % lower state
  usub = sl(2) + A * pslt;  % upper bound of upper state
  uslb = sl(2) - A * pslt;  % lower bound of upper state
  lsub = sl(1) + A * pslt;  % upper bound of lower state
  lslb = sl(1) - A * pslt;  % lower bound of lower state
  url = sl(1) + A * prl(3); % upper  reference level
  mrl = sl(1) + A * prl(2); % middle reference level
  lrl = sl(1) + A * prl(1); % lower  reference level

  if enableMeasureLines
    measureLines(1, 1:7) = getReferenceLevels(xLim, usub, uslb, url, mrl, lrl, lsub, lslb);
  end      
          
  [tm, iPre, iPost] = signal.internal.getTransitions(x, t, uslb, lsub, url, mrl, lrl);
  n = numel(iPre);

  if n>0
    % find all mid crossings within zoom region
    idx = find(xLim(1)<=tm.MiddleCross & tm.MiddleCross<=xLim(2));

    % assume all transitions are completely computed (at first)
    isComplete = true;
    
    if enableAberrations
      sd = settlingSeekDuration;
      prem  = signal.internal.getPreshoots( x, t, us, ls, uslb, lsub, 3, tm, iPre, iPost);
      [postm, isPostComplete] = signal.internal.getPostshoots(x, t, us, ls, uslb, lsub, 3, tm, iPre, iPost, false);
      [sm, isSettleComplete]  = signal.internal.getSettling(  x, t, us, ls, pslt, sd, tm, iPre, iPost, false);
      
      % indicate if an aberration cannot be computed for any given transition      
      isComplete = isPostComplete && isSettleComplete;
      
      readoutAberrations(1, :) = getReadoutAberrations(xLim, tm, prem, postm, sm, isComplete);
      if enableMeasureLines
        % obtain the undershoot, overshoot, and settling points
        measureLines(1, 11:15) = getAberrations(prem, postm, sm, isComplete);
      end          
    end

    % now get the transition metrics readout panel info
    readoutTransitions(1, :) = getReadoutTransitions(idx, sl, tm, isComplete);
    
    if enableMeasureLines
      measureLines(1, 8:10) = getReferenceCrossings(sl, tm, prl);
    end
    
    if enableCycles
      readoutCycles = getReadoutCycles(idx, tm, isComplete);
    end
  end
  
end

% -- auxiliary functions --
function stateLevels = computeStateLevels(x)
  lower = min(x);
  upper = max(x);
  histogram = signal.internal.getHistogram(x, 100, lower, upper);
  stateLevels = signal.internal.getLevelsByHistogram(histogram, ...
                    lower, upper, 'Histogram mode');
end

function m = safemean(x)
  if numel(x) > 0
    m = mean(x);
  else
    m = NaN;
  end
end

function c = vvIdx(A,index)
% handle Vector-Vector Indexing 
%
% In vector-vector indexing, you use one vector as an index into another
% vector. When either vector is variable sized, you might get a run-time
% error during code generation. Consider the index expression A(B). The
% general rule for indexing is that size(A(B)) == size(B). However, when
% both A and B are vectors, MATLAB applies a special rule: use the
% orientation of A as the orientation of the output. For example, if
% size(A) == [1 5] and size(B) == [3 1], then size(A(B)) == [1 3].
% 
% In this situation, if the code generation software detects that both A
% and B are vectors at compile time, it applies the special rule and gives
% the same result as MATLAB. However, if either A or B is a variable-size
% matrix (has shape ?x?) at compile time, the code generation software
% applies only the general indexing rule. Then, if both A and B become
% vectors at run time, the code generation software reports a run-time
% error when you run the MEX function.
  Avec = A(:);
  c = Avec(index(:));
end

function readout = getReadoutTransitions(idx, sl, tm, isComplete)
  readout = NaN(1,9);
  
  % get state levels and waveform amplitude
  readout(1) = sl(2);
  readout(2) = sl(1);
  readout(3) = diff(sl);
  
  % skip last transition if its aberration metrics could not be computed
  if ~isComplete
    idx = idx(1:end-1,1);
  end

  % find all rising transitions within zoom region
  pIdx = vvIdx(idx,tm.Polarity(idx)>0);
  readout(4) = numel(pIdx);
  readout(5) = safemean(vvIdx(tm.Duration,pIdx));
  readout(6) = safemean(vvIdx(tm.SlewRate,pIdx));

  % find all falling transitions within zoom region
  nIdx = vvIdx(idx,tm.Polarity(idx)<0);
  readout(7) = numel(nIdx);
  readout(8) = safemean(vvIdx(tm.Duration,nIdx));
  readout(9) = safemean(vvIdx(tm.SlewRate,nIdx));
end

function results = inBoundsPos(isComplete, xLim, polarity, metric, instant)
  % filter incoming metrics and instants based upon whether or not 
  % they are within bounds of xLim and have correct positive polarity

  polarity = polarity(:);
  metric = metric(:);
  instant = instant(:);
  
  if ~isComplete
    polarity = polarity(1:end-1,1);
    instant = instant(1:end-1,1);
    metric = metric(1:end-1,1);
  end
  
  isPositive = polarity > 0;
  isInBounds = xLim(1) <= instant & instant <= xLim(2);
  results = metric(isPositive & isInBounds);
end

function results = inBoundsNeg(isComplete, xLim, polarity, metric, instant)
  % filter incoming metrics and instants based upon whether or not 
  % they are within bounds of xLim and have correct negative polarity
  polarity = polarity(:);
  metric = metric(:);
  instant = instant(:);
  if ~isComplete
    polarity = polarity(1:end-1,1);
    instant = instant(1:end-1,1);
    metric = metric(1:end-1,1);
  end
  
  isnegative = polarity < 0;
  isInBounds = xLim(1) <= instant & instant <= xLim(2);
  results = metric(isnegative & isInBounds);
end

function readout = getReadoutAberrations(xLim, tm, prem, postm, sm, isComplete)
  % pre-initialize readout
  readout = NaN(1,8);

  % aberrations corresponding to positive transitions
  readout(1) = safemean(inBoundsPos(isComplete, xLim, tm.Polarity, prem.Undershoot,  prem.UndershootInstant));
  readout(2) = safemean(inBoundsPos(isComplete, xLim, tm.Polarity, postm.Overshoot,  postm.OvershootInstant));
  readout(3) = safemean(inBoundsPos(isComplete, xLim, tm.Polarity, postm.Undershoot, postm.UndershootInstant));
  readout(4) = safemean(inBoundsPos(isComplete, xLim, tm.Polarity, sm.Duration,      sm.Instant));

  % aberrations corresponding to negative transitions
  readout(5) = safemean(inBoundsNeg(isComplete, xLim, tm.Polarity, prem.Overshoot,   prem.OvershootInstant));
  readout(6) = safemean(inBoundsNeg(isComplete, xLim, tm.Polarity, postm.Overshoot,  postm.OvershootInstant));
  readout(7) = safemean(inBoundsNeg(isComplete, xLim, tm.Polarity, postm.Undershoot, postm.UndershootInstant));
  readout(8) = safemean(inBoundsNeg(isComplete, xLim, tm.Polarity, sm.Duration,      sm.Instant));
end

function results = cycleMetrics(firstIdx,cross)
  if ~isempty(firstIdx) && numel(cross)>firstIdx(1)+1
    idx = firstIdx(1);
    
    % count the number of final transitions
    numPulses  = numel(cross(idx+1:2:end,1));

    % subtract the final transitions from the initial transitions
    avgWidth   = safemean(cross(idx+1:2:end,1)-cross(idx:2:end-1,1));

    % dutycycle = (final - initial) ./ (next - initial) * 100%
    avgDutyCyc = 100*safemean((cross(idx+1:2:end-1,1)-cross(idx:2:end-2,1)) ...
                              ./(cross(idx+2:2:end,1)-cross(idx:2:end-2,1)));
    % combine for assignment
    results = [numPulses avgWidth avgDutyCyc];
  else
    results = [0 NaN NaN];
  end
end

function readout = getReadoutCycles(idx, tm, isComplete)
  % pre-initialize
  readout = NaN(1, 8);

  % don't report last transition if incomplete
  if isComplete
    midCross = vvIdx(tm.MiddleCross,idx);
  else
    midCross = vvIdx(tm.MiddleCross,idx(1:end-1,1));
  end

  % find average period of mid-crossings of polarities matching
  % the first transition found inside the zoom region
  readout(1) = safemean(diff(midCross(1:2:end,1)));
  readout(2) = 1 ./ safemean(diff(midCross(1:2:end,1)));
    
  %pulse polarity > 0
  pp = find(vvIdx(tm.Polarity,idx)>0, 1, 'first');
  readout(3:5) = cycleMetrics(pp,midCross);

  %pulse polarity < 0
  pp = find(vvIdx(tm.Polarity,idx)<0, 1, 'first');
  readout(6:8) = cycleMetrics(pp,midCross);
end

function points = getReferenceLevels(xLim, usub, uslb, url, mrl, lrl, lsub, lslb)
  coder.varsize('points.xData');
  coder.varsize('points.yData');
  points = repmat(struct('xData', [], 'yData', []), 1, 7);
  v = ones(size(xLim));
  for i=(1:7)
    points(i).xData = xLim;
  end
  points(1).yData = v * usub;
  points(2).yData = v * uslb;
  points(3).yData = v * url;
  points(4).yData = v * mrl;
  points(5).yData = v * lrl;
  points(6).yData = v * lsub;
  points(7).yData = v * lslb;
end

function points = getReferenceCrossings(sl, tm, prl)
  coder.varsize('points.xData');
  coder.varsize('points.yData');
  points = repmat(struct('xData', [], 'yData', []), 1, 3);
  A = diff(sl)/100;
  wa = A * ones(size(tm.MiddleCross));
  points(1).xData = tm.LowerCross;
  points(2).xData = tm.MiddleCross;
  points(3).xData = tm.UpperCross;
  points(1).yData = sl(1) + wa*prl(1);
  points(2).yData = sl(1) + wa*prl(2);
  points(3).yData = sl(1) + wa*prl(3);
end

function points = getAberrations(prem, postm, sm, isComplete)
  coder.varsize('points.xData');
  coder.varsize('points.yData');
  points = repmat(struct('xData', [], 'yData', []), 1, 5);
  
  if isComplete  
    points(1).xData = prem.UndershootInstant;
    points(2).xData = prem.OvershootInstant;
    points(3).xData = postm.UndershootInstant;
    points(4).xData = postm.OvershootInstant;
    points(5).xData = sm.Instant;
    points(1).yData = prem.UndershootLevel;
    points(2).yData = prem.OvershootLevel;
    points(3).yData = postm.UndershootLevel;
    points(4).yData = postm.OvershootLevel;
    points(5).yData = sm.Level;
  else % don't include last transition
    points(1).xData = prem.UndershootInstant(1:end-1,1);
    points(2).xData = prem.OvershootInstant(1:end-1,1);
    points(3).xData = postm.UndershootInstant(1:end-1,1);
    points(4).xData = postm.OvershootInstant(1:end-1,1);
    points(5).xData = sm.Instant(1:end-1,1);
    points(1).yData = prem.UndershootLevel(1:end-1,1);
    points(2).yData = prem.OvershootLevel(1:end-1,1);
    points(3).yData = postm.UndershootLevel(1:end-1,1);
    points(4).yData = postm.OvershootLevel(1:end-1,1);
    points(5).yData = sm.Level(1:end-1,1);
  end
end
