function mod06_02
  % Simulate separation of two trumpet notes.
  %
  % This Octave/MATLAB script was written to mimic a corresponding
  % Labview module (Module 6.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/>.
  %
  
  % Change these values to change the portion of the time waveform that
  % is displayed.  toff is the start time.  tdur is the duration
  toff = 100e-3;  tdur = 20e-3;

  % Filter order
  n = 9;

  % Fundamental Frequency
  fr = 436;
  
  % Decay constant for impulse response
  alpha = 100;
  
  % Save initial simulation state
  simState('toff', toff);
  simState('tdur', tdur); 
  simState('n', n);
  simState('fr', fr);
  simState('alpha', alpha);
  
  guiState('xup', true);
  guiState('yup', true);
  guiState('hup', true);
  guiState('xplay' , false);
  guiState('yplay' , false);
  guiState('s333h', -1);
  guiState('s334h', -1);
  guiState('s335h', -1);
  guiState('s336h', -1);

  clf;
  createPanel();
  createPlot();
end

function createPanel
  pos = [0.05 0.05 .9 .28 ];
  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
  %}

  uicontrol(p, 'style', 'text', ...
                'units', 'normalized', ...
                'string', 'A note = 440 Hz,   G note = 392 Hz', ...
                'position', [0.05  0.85 .6 .1]);
  %%{
  bg = uibuttongroup(p, 'position', [.05 .5 0.6 .2]);
  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]);
  uicontrol(bg, 'style', 'text', ...
                'units', 'normalized', ...
                'string', 'Play:', ...
                'position', [0.15  0.1 .2 .8]);
  
  guiState('bx', bx);
  guiState('by', by);
  
  h = sliderControl(p, 'Fund. Freq.', 300, ...
                      500, simState('fr'), [1 5], ...
                      [0.05 .2 0.6 .15], @guiCallback );
  guiState('frh', h);
  h = sliderControl(p, 'Order', 2, ...
                      10, simState('n'), [1 2], ...
                      [0.7 .55  .25 .15], @guiCallback );
  guiState('nh', h);
  h = sliderControl(p, 'alpha', 10, ...
                      400, simState('alpha'), [5 10], ...
                      [0.7 .25  .25 .15], @guiCallback );
  guiState('alphah', h);

  h = uicontrol(p, 'style', 'text', ...
                'units', 'normalized', ...
                'string', 'ready', ...
                'backgroundcolor', [.2 .2 .2], ...
                'foregroundcolor', [1 1 1], ...
                'position', [0.75  0.8 .15 .15]);
  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', [0 0 1 1]);
  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.frh)
    % fr slider (reson freq)
    v = round(get(s, 'value'));
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('fr', v);
    guiState('yup', true);
    guiState('hup', true);
  elseif (s == g.nh)
    % n slider  (fitler order)
    v = round(get(s, 'value'));
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('n', v);
    guiState('hup', true);
    guiState('yup', true);
  elseif (s == g.alphah)
    % alpha slider  
    v = round(get(s, 'value')/5)*5;
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('alpha', v);
    guiState('hup', true);
    guiState('yup', true);
  end
  
  createPlot();
end

function createPlot()

  %%{
  fs = 44100;
  % Keep x and y persistent so that they can be played without
  % having to be reread or recalculated
  persistent x 
  if (isempty(x))
      x = csvread('twotrumpetsGA.csv');
  end
  persistent y 
  if (isempty(y))
      y = 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 the output, y
  persistent FX FH

  %}
  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)');
  
    subplot(3,3,2);
    FX = fft(x);
    FXI = 2*abs(FX(fi))/L;
    semilogy(F/1000, FXI);
    xlim([0 5]);
    ylim([1e-6 1]);
    grid on;
    xlabel('freq (kHz)');
    ylabel('Amplitude');
    title('Original (x)');
  end
  %%{
  if(g.hup)
    guiState('hup', false);
    
    %{
    % Method 1 - Impulse Invariance
    % This is closest to the original Labview implementation
    % The Octave residue function gives a warning reguarding a singular
    % matrix for high order filters.
    tic
    N=poly(1i*2*pi*s.fr*[-s.n:-1 1:s.n]);
    D=poly(-s.alpha+1i*2*pi*s.fr*[-s.n:-1 1:s.n]);
    [R, P, K]=residue(N,D);
    h=real(R.'*exp(P*t));
    FH=1+fft(h)/fs;
    FHI=abs(FH(fi));
    [max(h) min(h) length(h)]
    toc
    %}
    
    %{
    % Method 2 - Modified Impulse Invariance
    % The Octave residue function performs better when working with a
    % normalized transfer function.  However, Octave still gives a warning 
    % for a 10th order filter.  This is slower than method 1.
    tic
    N=poly(1i*(-s.n:s.n));
    D=poly(-s.alpha/(2*pi*s.fr)+1i*(-s.n:s.n));
    % Omit zero/pole pair near origin
    %N=poly(1i*[-s.n:-1 1:s.n]);
    %D=poly(-s.alpha/(2*pi*s.fr)+1i*[-s.n:-1 1:s.n]);
    [R, P, K]=residue(N,D);
    % Denormalize R and P values
    h=real((2*pi*s.fr*R).'*exp(2*pi*s.fr*P*t));
    FH=1+fft(h)/fs;
    FHI=abs(FH(fi));
    [max(h) min(h) length(h)]
    toc
    %}
    
    %%{
    % Method 3
    % This method creates the frequency response by multiplying (in the
    % frequency domain) the individual (complex-conjugate pair) frequency
    % response terms.
    % This method is not only better conditioned than either of the
    % other methods (I tested up to 12th order filter), but much faster
    % too (in both Octave and MATLAB), a win-win.
    %tic
    % Uncomment the following line and comment out the one after it
    % to add a zero at the origin.  The other two methods have zeros
    % at the origin (but probably shouldn't)
    %FH = 1 - fft(s.alpha*exp(-s.alpha*t))/fs;
    FH = ones(size(t));
    for k = 1:s.n
      ghp = 2*s.alpha*exp(-s.alpha*t).*cos(k*2*pi*s.fr*t) - ...
            (s.alpha^2/(k*2*pi*s.fr))*exp(-s.alpha*t).*sin(k*2*pi*s.fr*t);
      FH = FH .* (1-fft(ghp)/(fs));
    end
    h = real(ifft(FH))*fs;
    % The large DC value swamps the rest of the impulse response.
    % Remove this value for plotting.  This is consistent with
    % the way the impulse response is shown in the other methods.
    h(1) = 0;
    FHI=abs(FH(fi));
    %[max(h) min(h) length(h)]
    %toc
    %}
    
    %%{
    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 Resp');
    else
      set(g.s3p1, 'ydata', h(1:round(L/15)));
    end  
    %}
    %%{
    if(~ishandle(g.s336h))
      sh = subplot(3,3,6);
      guiState('s336h', sh);
      lh = plot(F/1000, FHI);
      guiState('s6p1', lh);
      xlim([0 5]);
      ylim([1e-5 1]);
      grid on;
      xlabel('freq (kHz)');
      ylabel('Amplitude');
      title('Freq Resp');
    else
      set(g.s6p1, 'ydata', FHI);
    end  
    %}
  
  end

  if(g.yup)  
    %'y redraw'
    guiState('yup', false);
    y=real(ifft(FX.*FH));

    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('Filtered (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-6 1]);
      grid on;
      xlabel('freq (kHz)');
      ylabel('Amplitude');
      title('Filtered (y)');
    else
      set(g.s5p1, 'ydata', FYI);
    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

  %}
  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