function mod06_01
  % Demonstrate Sinusoidal Interference Removal
  %
  % This Octave/MATLAB script was written to mimic a corresponding
  % Labview module (Module 6.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/>.
  
  % Change these values to change the portion of the time waveform that
  % is displayed.  toff is the start time.  tdur is the duration
  % These should be GUI controls, but they don't fit in the current
  % design
  toff = 100e-3;  tdur = 20e-3;

  % Interference level
  Ai = 0.3; 
  % Interference and notch frequencies
  fi = 800; fn = 800; 

  % Decay constant for impulse response
  alpha = 200;
  
  % Save initial simulation state
  simState('Ai', Ai); simState('fi', fi); simState('fn', fn);
  simState('alpha', alpha);
  simState('toff', toff);
  simState('tdur', tdur);
  
  guiState('xup', true);
  guiState('yup', true);
  guiState('zup', true);
  guiState('hup', true);
  guiState('xplay' , false);
  guiState('yplay' , false);
  guiState('zplay' , false);
  guiState('s333h', -1);
  guiState('s334h', -1);
  guiState('s335h', -1);
  guiState('s337h', -1);
  guiState('s338h', -1);

  clf;
  createPanel();
  createPlot();
end

function createPanel
  pos = [.67 0.1 .30 .50 ];
  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
  
  bg = uibuttongroup(p, 'position', [.0 .88 1 .12]);
  bx = uicontrol(bg, 'style', 'pushbutton', ...
                'units', 'normalized', ...
                'string', 'x', ...
                'callback', @guiCallback, ...
                'position', [0.45 0.1 .1 .8]);
  by = uicontrol(bg, 'style', 'pushbutton', ...
                'units', 'normalized', ...
                'string', 'y', ...
                'callback', @guiCallback, ...
                'position', [0.65 0.1 .1 .8]);
  bz = uicontrol(bg, 'style', 'pushbutton', ...
                'units', 'normalized', ...
                'string', 'z', ...
                'callback', @guiCallback, ...
                'position', [0.85 0.1 .1 .8]);
  uicontrol(bg, 'style', 'text', ...
                'units', 'normalized', ...
                'string', 'Play:', ...
                'position', [0.05  0.1 .2 .8]);
  
  guiState('bx', bx);
  guiState('by', by);
  guiState('bz', bz);
  
  p = uipanel (p, 'position', [0 0 1 .88 ]);
  h = sliderControl(p, 'Ai', 0, ...
                      1, simState('Ai'), [.05 .2], ...
                      [0.1 .9 1 .1], @guiCallback );
  guiState('Aih', h);
  h = sliderControl(p, 'fi', 600, ...
                      2000, simState('fi'), [10 50], ...
                      [0.1 .7 1 .1], @guiCallback );
  guiState('fih', h);
  h = sliderControl(p, 'fn', 600, ...
                      2000, simState('fn'), [10 50], ...
                      [0.1 .5 1 .1], @guiCallback );
  guiState('fnh', h);
  h = sliderControl(p, 'alpha', 10, ...
                      500, simState('alpha'), [10 50], ...
                      [0.1 .3 1 .1], @guiCallback );
  guiState('ah', h);

  h = uicontrol(p, 'style', 'text', ...
                'units', 'normalized', ...
                'string', 'ready', ...
                'backgroundcolor', [.2 .2 .2], ...
                'foregroundcolor', [1 1 1], ...
                'position', [0.35  0.02 .3 .1]);
  guiState('status', h);
end

% Reset GUI panel after window size change
function resetControlPanel(s, e, 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', [1 1 0 0]);
  set(p, 'position', pos);
  set(p, 'visible', 'on');
end

function guiCallback(s, ~)
  g = guiState();
  if (s == g.bx)
    guiState('xplay', true);
  elseif (s == g.by)
    guiState('yplay', true);
  elseif (s == g.bz)
    guiState('zplay', true);
  elseif (s == g.Aih)
    % Ai slider (inter amplitude)
    v = round(get(s, 'value')/0.05)*0.05;
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('Ai', v);
    guiState('yup', true);
    guiState('zup', true);
  elseif (s == g.fih)
    % fi slider  (interference freq)
    v = round(get(s, 'value')/10)*10;
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('fi', v);
    guiState('yup', true);
    guiState('zup', true);
  elseif (s == g.fnh)
    % fn slider  (notch freq)
    v = round(get(s, 'value')/10)*10;
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('fn', v);
    guiState('hup', true);
    guiState('zup', true);
  elseif (s == g.ah)
    % alpha slider  (decay)
    v = round(get(s, 'value')/10)*10;
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('alpha', v);
    guiState('hup', true);
    guiState('zup', true);
  end
  
  createPlot();
end

function createPlot()

  %%{
  fs = 44100;
  % Keep x, y, and z persistent so that they can be played without
  % having to be recalculated
  persistent x 
  if (isempty(x))
      x = csvread('trumpet.csv');
  end
  persistent y 
  if(isempty(y))
      y = zeros(size(x));
  end
  persistent z 
  if (isempty(z))
      z = zeros(size(x));
  end
  
  % These are persistent for performance reasons.
  % They only need to be computed once and never change.
  persistent L 
  if (isempty(L))
      L = length(x);
  end
  persistent t 
  if (isempty(t))
      t = (0:L-1)/fs;
  end
  persistent fi
  if (isempty(fi))
      fi =(1:4000);
  end
  persistent F
  if (isempty(F))
      F =(fi-1)*fs/L;
  end
  
  % These are needed to compute z and FZ
  persistent FH FY

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

  i_sta = round(s.toff*fs)+1;
  i_end = i_sta + round(s.tdur*fs);
  
  nplot = i_sta:i_end;
  tplot = (nplot-1)*1000/fs;
  
  if(g.xup)
    guiState('xup', false);
    subplot(3,3,1);
    plot(tplot, x(nplot));
    grid on;
    xlim([s.toff s.toff+s.tdur]*1000);
    xlabel('time (ms)');
    ylabel('Amplitude');
    title('Original (x)');
  
    FX = fft(x);
    FXI = 2*abs(FX(fi))/L;
    subplot(3,3,2);
    semilogy(F/1000, FXI);
    xlim([0 5]);
    ylim([1e-5 1]);
    grid on;
    xlabel('freq (kHz)');
    ylabel('Amplitude');
    title('Original (x)');
  end
  
  if(g.hup)
    guiState('hup', false);
    alpha2=s.alpha*s.alpha;
    wo=2*pi*s.fn;
    A=sqrt(4*alpha2+alpha2*alpha2/wo/wo);
    theta=atan(s.alpha/wo/2);
    h=-A*exp(-s.alpha*t).*cos(wo*t+theta);
    FH=fft(h)/fs;

    if(~ishandle(g.s333h))
      sh = subplot(3,3,3);
      guiState('s333h', sh);
      lh = plot(t(1:round(L/15))*1000, h(1:round(L/15)));
      guiState('s3p1', lh);
      grid on;
      %xlim([0 tplotmax]);
      xlabel('time (ms)');
      ylabel('Amplitude');
      title('Impulse Response');
    else
      set(g.s3p1, 'ydata', h(1:round(L/15)));
    end  
  end

  if(g.yup)  
    % 'redraw'
    guiState('yup', false);
    y=x+s.Ai *cos(2*pi*s.fi*t);  
    
    if(~ishandle(g.s334h))
      sh = subplot(3,3,4);
      guiState('s334h', sh);    
      lh = plot(tplot, y(nplot));
      guiState('s4p1', lh);
      grid on;
      xlim([s.toff s.toff+s.tdur]*1000);
      xlabel('time (ms)');
      ylabel('Amplitude');
      title('Orig + Inter (y)');
    else
      set(g.s4p1, 'ydata', y(nplot));
    end  
    
    FY = fft(y);
    FYI = 2*abs(FY(fi))/L;
    if(~ishandle(g.s335h))
      sh = subplot(3, 3, 5);
      guiState('s335h', sh);    
      lh = semilogy(F/1000, FYI);
      guiState('s5p1', lh);
      xlim([0 5]);
      ylim([1e-5 1]);
      grid on;
      xlabel('freq (kHz)');
      ylabel('Amplitude');
      title('Orig + Inter (y)');
    else
      set(g.s5p1, 'ydata', FYI);
    end  
  end

  if(g.zup)
    guiState('zup', false);
    z=real(ifft(FY.*FH));
    z=y+z;

    if(~ishandle(g.s337h))
      sh = subplot(3, 3, 7);
      guiState('s337h', sh);    
      lh = plot(tplot, z(nplot));
      guiState('s7p1', lh);
      grid on;
      xlim([s.toff s.toff+s.tdur]*1000);
      xlabel('time (ms)');
      ylabel('Amplitude');
      title('Filtered (z)');
    else
      set(g.s7p1, 'ydata', z(nplot));
    end  
    
    FZ = fft(z);
    FZI = 2*abs(FZ(fi))/L;
    FHI=abs(1+FH(fi));
    if(~ishandle(g.s338h))
      sh = subplot(3, 3, 8);
      guiState('s338h', sh);    
      lh = semilogy(F/1000, FZI);
      guiState('s8p1', lh);
      hold on;
      lh = semilogy(F/1000, FHI, 'r');
      guiState('s8p2', lh);
      hold off;
      xlim([0 5]);
      ylim([1e-5 1]);
      grid on;
      xlabel('freq (kHz)');
      ylabel('Amplitude');
      title('Filtered (z)');
    else
      set(g.s8p1, 'ydata', FZI);
      set(g.s8p2, 'ydata', FHI);
    end  
  end
  
  if(g.xplay)
    xplayer = audioplayer(x, fs);
    playblocking(xplayer);
    guiState('xplay', false);
  end

  if(g.yplay)
    yplayer = audioplayer(y, fs);
    playblocking(yplayer);
    guiState('yplay', false);
  end

  if(g.zplay)
    zplayer = audioplayer(z, fs);
    playblocking(zplayer);
    guiState('zplay', false);
  end

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