function mod02_01
  % Illustrate convolution of two causal exponentials
  %
  % This Octave/MATLAB script was written to mimic a corresponding
  % Labview module (Module 2.1) that accompanies the free text 
  % "Signals and Systems: Theory and Applications" by Ulaby and Yagle
  % (ss2.eecs.umich.edu).
  %
  % In almost all cases the compuation code used in these scripts
  % is the same as that used in the Labview module.  Credit for the
  % computation code should go to the textbook authors.
  
  % These scripts have been tested under Octave 4.2.2, Octave 4.4.0
  % and MATLAB R2015a.  The Octave GUI implementation is still a
  % little buggy.  Some of strange plotting code used here attempts
  % to work around a few Octave bugs.

  % Copyright (C) 2018 Anthony Richardson <richardson.tony@gmail.com>
  %
  % This program is free software; you can redistribute it and/or modify it under
  % the terms of the GNU General Public License as published by the Free Software
  % Foundation; either version 3 of the License, or (at your option) any later
  % version.
  %
  % This program is distributed in the hope that it will be useful, but WITHOUT
  % ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  % FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  % details.
  %
  % You should have received a copy of the GNU General Public License along with
  % this program; if not, see <http://www.gnu.org/licenses/>.
  %

  % TMAX is the plot duration of the input
  % FS is the sampling frequency for both the plot and the
  % numerical convolution
  TMAX = 5;  FS = 100;
  
  % Initial simulation values
  alpha1 = 1.5;  alpha2 = 2;
  extmode = false;

  % Save all of these values so that the can be recalled later
  % All simulation state is saved and retrieved via the simState function
  simState('TMAX', TMAX);
  simState('FS', FS);
  simState('alpha1', alpha1);
  simState('alpha2', alpha2);
  simState('extmode', extmode);
  guiState('s221h', -1);
  guiState('s222h', -1);
  guiState('s223h', -1);
  
  clf;
  % Create the GUI control panel
  createGUI();
  % Create the initial set of plots
  createPlots();
end  

function createGUI
  % Create the UI elements (a checkbox and two sliders)
  pos = [.575 .11 .33 .325];
  p = uipanel ('position', pos, 'visible', 'off');
  guiState('cpanel', p);
  
  %% BUG workaround - may be able to remove this in a
  %% future version of octave
  if (exist('OCTAVE_VERSION', 'builtin'))
    % The panel position seems to get messed up when the
    % window is resized in Octave. This was an attempt to fix it.
    %set(gcf, 'sizechangedfcn', @resetControlPanel);
    % That didn't work but this seems too.
    dellistener(gcf, 'position');
    addlistener (gcf, 'position', {@resetControlPanel, p, pos});
  end
  
  cb = uicontrol(p, 'style', 'checkbox',...
        'units', 'normalized', 'string', 'Extended Mode',...
        'value', simState('extmode'), 'callback',@guiCallback,...
        'position',[0.12 0.8 0.7 0.1]);
  guiState('cb', cb);
  h = sliderControl(p, 'Alpha 1', 1, 10, simState('alpha1'), [.1 1], [.1 .6 .9 0.1], @guiCallback);
  guiState('ha1', h);
  h = sliderControl(p, 'Alpha 2', 1, 10, simState('alpha2'),  [.1 1], [.1 .3 .9 0.1], @guiCallback);
  guiState('ha2', h);
end

% The panel position needs to be reset when the window is moved
% or resized.  This is to work around a bug in Octave
function resetControlPanel(s, e, p, pos)
  set(p, 'visible', 'off');
  set(p, 'position', [0 0 1 1]);
  set(p, 'position', pos);
  set(p, 'visible', 'on');
end

function guiCallback(s, e)
  % Only ha1 and ha2 use this callback
  if (s == guiState('ha1'))
    v = round(get(s, 'value')*10)/10;
    set(s, 'value', v);
    %get(s, 'my_slider_val')
    %set(s,'my_slider_val'), 'string', num2str(v));
    set(get(s,'userdata'), 'string', num2str(v));
    simState('alpha1', v);
  elseif (s == guiState('ha2'))
    v = round(get(s, 'value')*10)/10;
    set(s, 'value', v);
    %set(get(s,'my_slider_val'), 'string', num2str(v));
    set(get(s,'userdata'), 'string', num2str(v));
    simState('alpha2', v);
  elseif (s == guiState('cb'))
    v = get(s, 'value');
    simState('extmode', v);
  end
  createPlots();
end

function createPlots()
  TMAX = simState('TMAX'); FS = simState('FS');
  % The t vector is only computed once for efficiency.
  persistent t;
  if(isempty(t))
    t = [0:1/FS:2*TMAX];
  end
  hideline = false;
  
  alpha1 = simState('alpha1');
  alpha2 = simState('alpha2');
  g = guiState();

  % Create numerical conv plot
  r = conv(exp(-alpha1*t),exp(-alpha2*t))/FS;
  % Plot r versus t in subplot(2,2,2)
  if(~ishandle(g.s222h))
    sh = subplot(2,2,2);
    guiState('s222h', sh);
    lh = plot(t, r(1:length(t)), 'linewidth', 4);
    guiState('s2p1', lh);
    ymaxv = max(r(:));
    ind = max(1,ceil(log2(ymaxv/0.05))+1);
    ymax = [0.05 0.1 0.2 0.4];
    ymaxt = ymax(ind);
  
    ylim([ 0 ymaxt]);
    set(gca, 'ytick', 0:ymaxt/10:ymaxt);
    set(gca, 'xtick', 0:1:10);
    % Ugly, but works in both MATLAB and Octave
    set (gca, 'yticklabel', ...
           arrayfun (@(x) sprintf ('%5.3f', x), ...
           get (gca, 'ytick'), 'uniformoutput', false)); 
    grid on;
    title('Convolution (Numeric)');
    xlabel('Time (s)');
    ylabel('Amplitude');
  else
    set(g.s2p1, 'ydata', r(1:length(t)));
  end  
    
  % Create theoretical conv plot
  if(alpha1 == alpha2)
    % In extended mode calculate correct response when alpha1 == alpha2
    r = t.*exp(-alpha1*t);
    if(~simState('extmode'))
      % In non-extended mode, hide the line
      hideline = true;
    end
  else
    r = 1/(alpha2-alpha1)*(exp(-alpha1*t) - exp(-alpha2*t));
  end

  % Plot r versus t in subplot(2,2,1)
  if(~ishandle(g.s221h))
    sh = subplot(2,2,1);
    guiState('s221h', sh);
    lh = plot(t, r, 'linewidth', 4);
    guiState('s1p1', lh);
    g.s1p1 = lh;
    ylim([ 0 ymaxt]);
    set(gca, 'ytick', 0:ymaxt/10:ymaxt);
    set(gca, 'xtick', 0:1:10);
    set(gca, 'yticklabel', ...
        arrayfun (@(x) sprintf ('%5.3f', x), ...
        get (gca, 'ytick'), 'uniformoutput', false)); 
    grid on;
    title('Convolution (Theoretical)');
    xlabel('Time (s)');
    ylabel('Amplitude');
  else
    set(g.s1p1, 'ydata', r);
  end  
  
  if(hideline)
      set(g.s1p1, 'visible', 'off');
  else
      set(g.s1p1, 'visible', 'on');
  end
  
  % Show the input signals
  h = exp(-alpha1*t(1:floor(end/2)));
  x = exp(-alpha2*t(1:floor(end/2)));
  % Plot h and x versus t in subplot(2,2,3)
  if(~ishandle(g.s223h))
    sh = subplot(2,2,3);
    guiState('s223h', sh);
    % Use only half the time vector
    lh = plot(t(1:floor(end/2)), h, 'linewidth', 4);
    guiState('s3p1', lh);
    hold on;
    lh = plot(t(1:floor(end/2)), x, 'linewidth', 4);
    guiState('s3p2', lh);
    hold off
    set(gca, 'xtick', 0:0.5:5);
    set(gca, 'ytick', 0:0.1:1);
    set (gca, 'yticklabel', ...
         arrayfun (@(x) sprintf ('%3.1f', x), ...
         get (gca, 'ytick'), 'uniformoutput', false)); 
    grid on;
    title('Input Signals');
    xlabel('Time (s)');
    ylabel('Amplitude');
  else
    set(g.s3p1, 'ydata', h);
    set(g.s3p2, 'ydata', x);
  end  
  legend(['\alpha_1 = ' sprintf('%4.1f', alpha1)],...
         ['\alpha_2 = ' sprintf('%4.1f', alpha2)]);
  set(guiState('cpanel'), 'visible', 'on');
end

function v = guiState(key, val)
  persistent state;
  if (nargin == 0)  
    v = state;
  elseif (nargin == 1)
    v = state.(key);
  elseif (nargin == 2)
    state.(key) = val;
    v = val;
  end
end

function v = simState(key, val)
  persistent state;
  if (nargin == 0)  
    % Return the entire struct
    v = state;
  elseif (nargin == 1)
    % Return a particular field from the struct
    v = state.(key);
  elseif (nargin == 2)
    % Set a key value (key must be a string)
    state.(key) = val;
    v = val;
  end
end

function h = sliderControl(parent, label, min, max, val, step, pos, cback)
  % sliderstep = [minstep maxstep]
  % where minstep is fraction of range on arrow click
  % and maxstep is fraction of range on trough click
  h = uicontrol(parent, 'style' ,'slider',...
           'units', 'normalized',...
           'position', [pos(1)+0.1*pos(3) pos(2) 0.7*pos(3) pos(4)],...
           'min', min, 'max', max,...
           'sliderstep', [step(1)/(max-min) step(2)/(max-min)],...
           'value', val, 'callback', cback);
  uicontrol(parent, 'style', 'text',...
           'units', 'normalized',...
           'position', [pos(1) pos(2) 0.1*pos(3) pos(4)],...
           'String', num2str(min));
  %addproperty('my_slider_min', h, 'handle', hc);
  uicontrol(parent, 'style', 'text',...
           'units', 'normalized',...
           'position', [pos(1)+0.8*pos(3) pos(2) 0.1*pos(3) pos(4)],...
           'String', num2str(max));         
  %addproperty('my_slider_max', h, 'handle', hc);
  uicontrol(parent, 'style', 'text',...
           'units', 'normalized',...
           'position', [pos(1)+0.15*pos(3) pos(2)-pos(4) 0.4*pos(3) pos(4)],...
           'horizontalalignment','left',...
           'String', label);        
  hc = uicontrol(parent, 'style', 'text',...
           'units', 'normalized',...
           'position', [pos(1)+.5*pos(3) pos(2)-pos(4) 0.3*pos(3) pos(4)],...
           'String', num2str(val));
  %addproperty('my_slider_val', h, 'handle', hc);
  set(h, 'userdata', hc);
end