function mod08_01
  % Filter Design by Pole Zero Placement
  %
  % This Octave/MATLAB script was written to mimic a corresponding
  % Labview module (Module 8.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/>.
  
  %%{
  % Initial Pole/zero Radius
  pz_r = [ 1.00 1.00 1.00 0.60 0.80 0.80 ];
  % Initial Pole/zero Angle (mult by pi to get angle in radians)
  pz_a = [ 0.50 0.75 1.00 0.00 0.25 0.50 ];
  % Initial pole/zero state (zero if true, pole if false)
  pz_s = [ true true true false false false ];
  % Initial pole/zero inclusion state (include if true, exclude if false)
  pz_i = [ true true true true true true ];
  %}
  
  % Initial gain term
  C = 1;
  %{
  % Initial Pole/zero Radius
  pz_r = [ 1.00 0.50 0.50 0.50 0.00 0.00 ];
  % Initial Pole/zero Angle (mult by pi to get angle in radians)
  pz_a = [ 139/180 60/180 1.00 0.00 0.25 0.50 ];
  % Initial pole/zero state (zero if true, pole if false)
  pz_s = [ true false false false false false ];
  % Initial pole/zero inclusion state (include if true, exclude if false)
  pz_i = [ true true false false false false ];
  %}
  
  % Save initial simulation state
  simState('pz_r', pz_r);
  simState('pz_a', pz_a); 
  simState('pz_s', pz_s);
  simState('pz_i', pz_i);
  simState('C', C);
  % Initial cursor angle
  simState('cur_ang', 0);
  simState('cur_ang_inc',0);
  simState('pzr_inc', 0);
  simState('pza_inc', 0);
  
  simState('pup', true);    % Update plot
  simState('cup', true);    % Update cursor position
  clf;
  createPanel();
  createPlot();
end

function createPanel
  s = simState();
  pos = [0.01 0.01 .98 .28 ];
  p = uipanel ('position', pos, 'visible', 'off');
  simState('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

  % subpanel geometry
  wi = 0.298; hi = 0.42;
  xoff = [0.0 0.3 0.6 0.0 0.3 0.6];
  yoff = [0.43 0.43 0.43 0.0 0.0 0.0];
  radius_handles = zeros(1, 6);
  angle_handles = zeros(1, 6);
  enable_handles = zeros(1, 6);
  state_handles = zeros(1, 6);
  for k = 1:6
    pk = uipanel (p, 'position', [xoff(k) yoff(k) wi hi]);
    
    % Enable toggle button
    if(s.pz_i(k))
      str = 'Enable';
      bgc = get(0, 'defaultuipanelbackgroundcolor');
    else
      str = 'Disable';
      bgc = [1 1 0.2];
    end
    h = uicontrol(pk, 'style', 'togglebutton', ...
                      'units', 'normalized', ...
                      'position', [0 0.5 0.25 0.5], ...
                      'callback', @pziCallback, ...
                      'string', str, ...
                      'backgroundcolor', bgc, ...
                      'value', s.pz_i(k));
    enable_handles(k) = h;
    
    % Pole/Zero toggle button
    if(s.pz_s(k))
      str = 'Zero';
      bgc = get(0, 'defaultuipanelbackgroundcolor');
    else
      str = 'Pole';
      bgc = [1 0.2 0.2];
    end
    h = uicontrol(pk, 'style', 'togglebutton', ...
                      'units', 'normalized', ...
                      'position', [0 0.0 0.25 0.5], ...
                      'callback', @pzsCallback, ...
                      'string', str, ...
                      'backgroundcolor', bgc, ...
                      'value', s.pz_s(k));
    state_handles(k) = h;
    
    % Radius control (with label and +/-
    uicontrol(pk, 'style', 'text', ...
                      'units', 'normalized', ...
                      'position', [0.275 0.5 0.3 0.5], ...
                      'string', 'Radius');
    h = uicontrol(pk, 'style', 'edit', ...
                      'units', 'normalized', ...
                      'position', [0.6 0.5 0.3 0.5], ...
                      'callback', @pzrCallback, ...
                      'string', num2str(s.pz_r(k),'%4.2f'));
    radius_handles(k) = h;
    hp = uicontrol(pk, 'style', 'pushbutton', ...
                      'units', 'normalized', ...
                      'position', [0.9 0.75 0.1 0.25], ...
                      'callback', @pzrpCallback, ...
                      'string', '+');
    set(hp, 'userdata', h);
    hm = uicontrol(pk, 'style', 'pushbutton', ...
                      'units', 'normalized', ...
                      'position', [0.9 0.5 0.1 0.25], ...
                      'callback', @pzrmCallback, ...
                      'string', '-');
    set(hm, 'userdata', h);
    set(h, 'userdata', [hp hm]);
    
    % Angle control (with label and +/-
    uicontrol(pk, 'style', 'text', ...
                      'units', 'normalized', ...
                      'position', [0.275 0.0 0.3 0.5], ...
                      'string', 'Ang/pi');
    h = uicontrol(pk, 'style', 'edit', ...
                      'units', 'normalized', ...
                      'position', [0.6 0.0 0.3 0.5], ...
                      'callback', @pzaCallback, ...
                      'string', num2str(s.pz_a(k),'%4.2f'));
    angle_handles(k) = h;
    hp = uicontrol(pk, 'style', 'pushbutton', ...
                      'units', 'normalized', ...
                      'position', [0.9 0.25 0.1 0.25], ...
                      'callback', @pzapCallback, ...
                      'string', '+');
    set(hp, 'userdata', h);
    hm = uicontrol(pk, 'style', 'pushbutton', ...
                      'units', 'normalized', ...
                      'position', [0.9 0.0 0.1 0.25], ...
                      'callback', @pzamCallback, ...
                      'string', '-');
    set(hm, 'userdata', h);
    set(h, 'userdata', [hp hm]);
  end
  simState('enable_handles', enable_handles);
  simState('state_handles', state_handles);
  simState('radius_handles', radius_handles);
  simState('angle_handles', angle_handles);
  
  % Cursor panel
  cp = uipanel (p, 'position', [.9 0 .1 1]);
  uicontrol(cp, 'style', 'text', ...
                    'units', 'normalized', ...
                    'position', [0.0 0.75 1 0.25], ...
                    'horizontalalignment', 'center', ...
                    'string', 'Cursor');
  h = uicontrol(cp, 'style', 'text', ...
                    'units', 'normalized', ...
                    'position', [0.0 0.0 1 0.25], ...
                    'string', num2str(0,'%4.2f'));
  simState('cyh', h);
  h = uicontrol(cp, 'style', 'edit', ...
                    'units', 'normalized', ...
                    'position', [0.0 0.25 1 0.25], ...
                      'callback', @cursCallback, ...
                      'string', num2str(0,'%4.2f'));
  simState('cxh', h);
  hp = uicontrol(cp, 'style', 'pushbutton', ...
                    'units', 'normalized', ...
                    'position', [0.0 0.5 0.5 0.25], ...
                    'callback', @cpCallback, ...
                    'string', '+');
  set(hp, 'userdata', h);
  hm = uicontrol(cp, 'style', 'pushbutton', ...
                    'units', 'normalized', ...
                    'position', [0.5 0.5 0.5 0.25], ...
                    'callback', @cmCallback, ...
                    'string', '-');
  set(hm, 'userdata', h);
  set(h, 'userdata', [hp hm]);

  %%{
  % Gain Term
  uicontrol(p,  'style', 'text', ...
                'units', 'normalized', ...
                'position', [0.35 0.85 .03 0.15], ...
                'horizontalalignment', 'center', ...
                'string', 'C');
  h = uicontrol(p, 'style', 'edit', ...
                'units', 'normalized', ...
                'position', [0.4 0.85 .2 0.15], ...
                'callback', @gCallback, ...
                'string', num2str(s.C,'%f'));
  simState('gh', h);
  %}
  
end

% Reset GUI panel after window size change
function resetControlPanel(~, ~, p, pos)
  % Just setting the panel position does not seem to work but
  % maximizing it and then resetting the position does.
  set(p, 'visible', 'off');
  set(p, 'position', [0 0 1 1]);
  set(p, 'position', pos);
  set(p, 'visible', 'on');
end

% Enable/Disable Toggle Button Callback
function pziCallback(s, ~)
  pz_i = simState('pz_i');
  enable_handles = simState('enable_handles');
  for k = 1:length(enable_handles)
    if ( s == enable_handles(k))
      if(pz_i(k))
        pz_i(k) = false;
        set(enable_handles(k), 'string', 'Disable');
        set(enable_handles(k), 'backgroundcolor', [1 1 0.2]);
      else
        pz_i(k) = true;
        set(enable_handles(k), 'string', 'Enable');
        bgc = get(0, 'defaultuipanelbackgroundcolor');
        set(enable_handles(k), 'backgroundcolor', bgc);
      end
    end
  end
  simState('pz_i', pz_i);
  simState('pup', true);
  createPlot();
end

% Pole/Zero Toggle Button Callback
function pzsCallback(s, ~)
  pz_s = simState('pz_s');
  state_handles = simState('state_handles');
  for k = 1:length(state_handles)
    if ( s == state_handles(k))
      if(pz_s(k))
        pz_s(k) = false;
        set(state_handles(k), 'string', 'Pole');
        set(state_handles(k), 'backgroundcolor', [1 0.2 0.2]);
      else
        pz_s(k) = true;
        set(state_handles(k), 'string', 'Zero');
        bgc = get(0, 'defaultuipanelbackgroundcolor');
        set(state_handles(k), 'backgroundcolor', bgc);
      end
    end
  end
  simState('pz_s', pz_s);
  simState('pup', true);
  createPlot();
end

% Pole/Zero Radius Change Callback
function pzrCallback(s, ~)
  pz_r = simState('pz_r');
  radius_handles = simState('radius_handles');
  for k = 1:length(radius_handles)
    if ( s == radius_handles(k))
      r = str2double(get(radius_handles(k), 'string'));
      if(isnan(r))
          r = 1;
      end
      delta = simState('pzr_inc');
      r = r + delta*0.1;
      simState('pzr_inc', 0);
      pz_r(k) = r;
      set(radius_handles(k), 'string', num2str(r,'%4.2f'))
    end
  end
  simState('pz_r', pz_r);
  simState('pup', true);
  createPlot();
end

% Pole/Zero Radius Increment Button Callback
function pzrpCallback(s, ~)
  simState('pzr_inc', +1);
  pzrCallback(get(s, 'userdata'));
end

% Pole/Zero Radius Decrement Button Callback
function pzrmCallback(s, ~)
  simState('pzr_inc', -1);
  pzrCallback(get(s, 'userdata'));
end

% Pole/Zero Angle Change Callback
function pzaCallback(s, ~)
  pz_a = simState('pz_a');
  angle_handles = simState('angle_handles');
  for k = 1:length(angle_handles)
    if ( s == angle_handles(k))
      a = str2double(get(angle_handles(k), 'string'));
      if(isnan(a))
          a = 0;
      end
      delta = simState('pza_inc');
      % Radius increment or decrement
      a = a + delta*0.1;
      pz_a(k) = a;
      simState('pza_inc', 0);
      set(angle_handles(k), 'string', num2str(a,'%4.2f'))
    end
  end
  simState('pz_a', pz_a);
  simState('pup', true);
  createPlot();
end

% Pole/Zero Angle Increment Button Callback
function pzapCallback(s, ~)
  simState('pza_inc', +1);
  pzaCallback(get(s, 'userdata'));
end

% Pole/Zero Angle Decrement Button Callback
function pzamCallback(s, ~)
  simState('pzr_inc', -1);
  pzaCallback(get(s, 'userdata'));
end

% Cursor Change Callback
function cursCallback(s, ~)
  a = str2double(get(s,'string'));
  if(isnan(a))
      a = 0;
  end
  simState('cur_ang', a);
  createPlot();
end

function cpCallback(~, ~)
  simState('cur_ang_inc',1);
  createPlot();
end

function cmCallback(~, ~)
  simState('cur_ang_inc',-1);
  createPlot();
end

% gain Change Callback
function gCallback(s, ~)
  a = str2double(get(s,'string'));
  if(isnan(a))
      a = 1;
  end
  simState('C', a);
  set(s, 'string', num2str(a,'%f'));
  createPlot();
end

function createPlot()
  s = simState();
  
  persistent omega;
  if (isempty(omega))
    omega = (0:0.002:1)*pi;
  end
  
  persistent H;

  % Update cursor position  
  indarr = find(s.cur_ang <= omega);
  if(isempty(indarr))
    ind = length(omega);
  else
    ind = indarr(1);
  end
  ind = ind + s.cur_ang_inc*5;
  simState('cur_ang_inc', 0);
  if(ind < 1)
    ind = 1;
  end
  if(ind > length(omega))
    ind = length(omega);
  end
  set(s.cxh, 'string', num2str(omega(ind),'%6.3f'));
  s.cur_ang = omega(ind);
  simState('cur_ang', s.cur_ang);

  % Update plots
  simState('xup', false);
  
  % Pole/Zero Map
  % Completely recreate plot for any change.  This could probably be more
  % efficient by not recreating plot for a cursor change, but it is
  % still relatively fast.
  angle = (-1:0.01:1)*pi;
  subplot('position', [5/52 25/64 4.8/13 8/16]);
  % Plot unit circle
  plot(cos(angle), sin(angle));
  grid on;
  axis([-1.2 1.2 -1.2 1.2], 'square');
  box('on');
  hold on;
  
  % Mark poles and zeros
  % Extract enabled poles and zeros
  pz_r = s.pz_r(s.pz_i);
  pz_a = s.pz_a(s.pz_i);
  pz_s = s.pz_s(s.pz_i);
  % Separate zeros and poles from enabled pole/zero array
  z_r = pz_r(pz_s);
  z_a = pz_a(pz_s);
  p_r = pz_r(~pz_s);
  p_a = pz_a(~pz_s);
  % Separate real zeros from complex zeros
  real_z = or(z_a == 0, z_a == 1);
  % Real zeros
  plot(z_r(real_z).*cos(z_a(real_z)*pi), z_r(real_z).*sin(z_a(real_z)*pi), 'ok');
  % Conj pair zeros
  plot(z_r(~real_z).*cos(z_a(~real_z)*pi), z_r(~real_z).*sin(z_a(~real_z)*pi), 'ok');
  plot(z_r(~real_z).*cos(-z_a(~real_z)*pi), z_r(~real_z).*sin(-z_a(~real_z)*pi), 'ok');
  % Separate real poles from complex poles
  real_p = or(p_a == 0, p_a == 1);
  % Real zeros
  plot(p_r(real_p).*cos(p_a(real_p)*pi), p_r(real_p).*sin(p_a(real_p)*pi), 'xr');
  % Conj pair zeros
  plot(p_r(~real_p).*cos(p_a(~real_p)*pi), p_r(~real_p).*sin(p_a(~real_p)*pi), 'xr');
  plot(p_r(~real_p).*cos(-p_a(~real_p)*pi), p_r(~real_p).*sin(-p_a(~real_p)*pi), 'xr');
  
  % Plot cursor on pole/zero map
  plot([0 cos(s.cur_ang)], [0 sin(s.cur_ang)], 'r');
  hold off;
  xlabel('Real');
  ylabel('Imaginary');
  h = title('Pole-Zero Plot');
  p = get(h, 'position');
  xl = xlim();
  set(h, 'position', [xl(1) p(2:3)], 'horizontalalignment', 'left');
  
  %%{
  % Plot the magnitude response
  H = s.C*ones(size(omega));
  for k = 1:length(pz_r)
    if(or(pz_a(k) == 0, pz_a(k) == 1))
      %Real pole/zero
      if(pz_s(k))
        % real zero
        H = H .* (exp(1i*omega) - pz_r(k)*cos(pz_a(k)*pi));
      else
        % real pole
        H = H ./ (exp(1i*omega) - pz_r(k)*cos(pz_a(k)*pi));
      end
    else
      if(pz_s(k))
        % conj pair of zeros
        %H = H .* (exp(1i*2*angle) - ...
        %          2*s.pz_r(k)*cos(s.pz_a(k)*pi)*exp(1i*angle) + ...
        %          s.pz_r(k)^2);
        H = H .* (exp(1i*omega) - pz_r(k)*exp(+1i*pz_a(k)*pi));
        H = H .* (exp(1i*omega) - pz_r(k)*exp(-1i*pz_a(k)*pi));
      else
        % conj pair of poles
        %H = H ./ (exp(1i*2*angle) - ...
        %          2*s.pz_r(k)*cos(s.pz_a(k)*pi)*exp(1i*angle) + ...
        %          s.pz_r(k)^2);
        H = H ./ (exp(1i*omega) - pz_r(k)*exp(+1i*pz_a(k)*pi));
        H = H ./ (exp(1i*omega) - pz_r(k)*exp(-1i*pz_a(k)*pi));
      end
    end
  end
  
  % Find the response at the cursor
  set(s.cyh, 'string', num2str(abs(H(ind)),'%6.3f'));
  % Plot the response
  subplot('position', [29.7/54 25/64 5/13 8/16]);
  plot(omega, abs(H), 'linewidth', 2);
  box('on')
  grid on;
  hold on;
  ylims = ylim();
  % Plot the cursor
  plot([omega(ind) omega(ind)],[ 0 ylims(2)], 'r');
  hold off;
  xlabel('\Omega (rad)');
  ylabel('Magnitude');
  h = title('Frequency Response');
  p = get(h, 'position');
  xl = xlim();
  set(h, 'position', [xl(1) p(2:3)], 'horizontalalignment', 'left');
  
  set(s.cpanel, 'visible', 'on');
  % Create invisible axes for labels
  axes('position', [0 0 1 1], 'visible', 'off');
  text(0.5, 0.97, 'Filter Design by Pole-Zero Placement', ...
                  'fontsize', 16, ...
                  'fontweight', 'bold', ...
                  'horizontalalignment', 'center');
  drawnow();
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
