function mod08_02
  % DT Notch Filter to Remove One of Two Sinusoids
  %
  % This Octave/MATLAB script was written to mimic a corresponding
  % Labview module (Module 8.2) 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/>.
  
  f1 = 120;     % wave 1 freq
  f2 = 200;     % wave 2 freq
  fn = 200;     % notch freq
  rp = 0.90;    % pole radius
  
  fs = 1000;    % sampling frequency
  
  simState('f1', f1);
  simState('f2', f2);
  simState('fn', fn);
  simState('rp', rp);
  simState('fs', fs);
  
  clf;
  createPanel();
  createPlot();
end

function createPanel
  pos = [0.66 0.48 .32 .52 ];
  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

  h = sliderControl(p, 'f1', [100 300], simState('f1'), ...
                    [1 5], [0 0.78 1 .2 ], @guiCallback);
  simState('f1h', h);
  h = sliderControl(p, 'notch frequency', [100 300], simState('fn'), ...
                    [1 5], [0 0.52 1 .2 ], @guiCallback);
  simState('fnh', h);
  h = sliderControl(p, 'f2', [100 300], simState('f2'), ...
                    [1 5], [0 0.27 1 .2 ], @guiCallback);
  simState('f2h', h);
  h = sliderControl(p, 'pole radius', [0.5 0.99], simState('rp'), ...
                    [0.05 0.10], [0.0 0.02 1 .2 ], @guiCallback);
  simState('rph', 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

function guiCallback(s, ~)
  g = simState();
  if (s == g.f1h)
    % f1 slider 
    v = round(get(s, 'value'));
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('f1', v);
  elseif (s == g.fnh)
    % notch slider 
    v = round(get(s, 'value'));
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('fn', v);
  elseif (s == g.f2h)
    % f2 slider 
    v = round(get(s, 'value'));
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('f2', v);
  elseif (s == g.rph)
    % pole radius slider 
    v = round(get(s, 'value')/0.01)*0.01;
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('rp', v);
  end
  createPlot();
end

function createPlot()
  %}
  s = simState();
  %%{
  %set(s.status, 'string', 'wait');

  % Generate input signal
  n = 1:50;
  nplot = 30:length(n);
  %nplot = n;
  W1 = 2*pi*s.f1/s.fs;
  W2 = 2*pi*s.f2/s.fs;  
  x = sin(W1*n) + sin(W2*n);

  % Generate filtered output
  WN = 2*pi*s.fn/s.fs;
  B=[1 -2*cos(WN) 1];
  A=[1 -2*s.rp*cos(WN) s.rp*s.rp];
  y = filter(B, A, x);

  % Compute spectra of input and output
  FX=2*abs(fftshift(fft(x, 1000)))/50;
  FY=2*abs(fftshift(fft(y, 1000)))/50;
  k = linspace(-1,1,1000);
  fplot = k*s.fs/2;
  
  % Compute the response of the filter
  EW = exp(1j*pi*k);
  H=polyval(B,EW)./polyval(A,EW);
  
  % Generate the plots
  p_hi = 0.32;
  p_wi = 0.22;
  p_xoff = [0.08 0.41 0.08 0.41 0.74];
  p_yoff = [0.61 0.61 0.12 0.12 0.10];
  sCnt = 1;
  subplot('position', [p_xoff(sCnt) p_yoff(sCnt)  p_wi p_hi]);
  stem(nplot, x(nplot));
  grid on;
  xlabel('Samples');
  ylabel('Value');
  h = title('Original Signal (time)');
  p = get(h, 'position');
  xl = xlim();
  set(h, 'position', [xl(1) p(2:3)], 'horizontalalignment', 'left');

  sCnt = sCnt + 1;
  subplot('position', [p_xoff(sCnt) p_yoff(sCnt)  p_wi p_hi]);
  semilogy(fplot, FX);
  %yl = ylim();
  ylim([1e-3 1.2]);
  xlim([-s.fs/2 s.fs/2]);
  grid on;
  hold on;
  semilogy(fplot, abs(H));
  hold off;
  xlabel('Frequency (Hz)');
  ylabel('Value');
  h = title('Original Signal (freq)');
  p = get(h, 'position');
  xl = xlim();
  set(h, 'position', [xl(1) p(2:3)], 'horizontalalignment', 'left');

  sCnt = sCnt + 1;
  subplot('position', [p_xoff(sCnt) p_yoff(sCnt)  p_wi p_hi]);
  stem(nplot, y(nplot));
  grid on;
  xlabel('Samples');
  ylabel('Value');
  h = title('Filtered Signal (time)');
  p = get(h, 'position');
  xl = xlim();
  set(h, 'position', [xl(1) p(2:3)], 'horizontalalignment', 'left');

  sCnt = sCnt + 1;
  subplot('position', [p_xoff(sCnt) p_yoff(sCnt)  p_wi p_hi]);
  semilogy(fplot, FY);
  %yl = ylim();
  ylim([1e-3 1.2]);
  xlim([-s.fs/2 s.fs/2]);
  grid on;
  xlabel('Frequency (Hz)');
  ylabel('Value');
  h = title('Filtered Signal (freq)');
  p = get(h, 'position');
  xl = xlim();
  set(h, 'position', [xl(1) p(2:3)], 'horizontalalignment', 'left');

  sCnt = sCnt + 1;
  subplot('position', [p_xoff(sCnt) p_yoff(sCnt)  p_wi p_hi]);
  % Plot unit circle
  angle = (-1:0.01:1)*pi;
  plot(cos(angle), sin(angle));
  grid on;
  axis([-1.2 1.2 -1.2 1.2], 'square');
  box('on');
  hold on;
  plot([cos(WN) cos(-WN)], [sin(WN) sin(-WN)], 'ok');
  plot(s.rp*[cos(WN) cos(-WN)], s.rp*[sin(WN) sin(-WN)], 'xr');
  hold off;
  xlabel('Real');
  ylabel('Imag');
  h = title('Pole-Zero Plot');
  p = get(h, 'position');
  xl = xlim();
  set(h, 'position', [xl(1) p(2:3)], 'horizontalalignment', 'left');

  %}
  set(s.cpanel, 'visible', 'on');
  %set(s.status, 'string', 'ready');
  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

function h = sliderControl(parent, label, limits, 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) pos(2)+0.33*pos(4) pos(3) 0.32*pos(4)],...
           'min', limits(1), 'max', limits(2),...
           'sliderstep', [step(1)/(limits(2)-limits(1)) step(2)/(limits(2)-limits(1))],...
           'value', val, 'callback', cback);
  uicontrol(parent, 'style', 'text',...
           'units', 'normalized',...
           'position', [pos(1) pos(2) 0.2*pos(3) 0.3*pos(4)],...
           'horizontalalignment', 'left', ...
           'String', num2str(limits(1)));
  uicontrol(parent, 'style', 'text',...
           'units', 'normalized',...
           'position', [pos(1)+0.8*pos(3) pos(2) 0.2*pos(3) 0.3*pos(4)],...
           'horizontalalignment', 'right', ...
           'String', num2str(limits(2)));         
  uicontrol(parent, 'style', 'text',...
           'units', 'normalized',...
           'position', [pos(1) pos(2)+0.67*pos(4) pos(3) 0.3*pos(4)],...
           'horizontalalignment','left',...
           'String', label);        
  hc = uicontrol(parent, 'style', 'text',...
           'units', 'normalized',...
           'position', [pos(1)+.3*pos(3) pos(2) 0.4*pos(3) 0.3*pos(4)],...
           'backgroundcolor', [0.8 0.8 0.8], ...
           'String', num2str(val));
  set(h, 'userdata', hc);
end
