function mod08_05
  % 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 = 23e-3;  tdur = 5e-3;

  % Fundamental Frequency
  fr = 440;
  
  % Pole radius
  pr = 0.9;
  
  % Save initial simulation state
  simState('toff', toff);
  simState('tdur', tdur); 
  simState('fr', fr);
  simState('pr', pr);
  
  simState('xplay' , false);
  simState('yplay' , false);

  clf;
  createPanel();
  createPlot();
end

function createPanel
  pos = [0.7 0.02 .29 .96 ];
  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, 'Fundamental Frequency', [430, 510], ...
                      simState('fr'), [1 5], ...
                      [0.05 .75 0.9 .25], @guiCallback );
  simState('frh', h);
  h = sliderControl(p, 'Pole Radius', [0.1, 0.99], ...
                      simState('pr'), [0.05 0.1], ...
                      [0.05 .50  .9 .25], @guiCallback );
  simState('prh', h);

  bx = uicontrol(p, 'style', 'pushbutton', ...
                'units', 'normalized', ...
                'string', 'Play x', ...
                'callback', @guiCallback, ...
                'position', [0.1 0.3 .8 .15]);
  by = uicontrol(p, 'style', 'pushbutton', ...
                'units', 'normalized', ...
                'string', 'Play y', ...
                'callback', @guiCallback, ...
                'position', [0.1 0.1 .8 .15]);
  simState('bx', bx);
  simState('by', by);

  h = uicontrol(p, 'style', 'text', ...
                'units', 'normalized', ...
                'string', 'ready', ...
                'backgroundcolor', [.2 .2 .2], ...
                'foregroundcolor', [1 1 1], ...
                'position', [0.35  0.02 .3 .05]);
  simState('status', 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');
  drawnow();
end

function guiCallback(s, ~)
  g = simState();
  if (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);
  elseif (s == g.prh)
    % pole radius slider  
    v = round(get(s, 'value')/0.01)*0.01;
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('pr', v);
  elseif (s == g.bx)
    simState('xplay', true);
  elseif (s == g.by)
    simState('yplay', 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('twotrumpetsAB.csv');
  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
  
  fplot_max = 5000;
  persistent fi
  if (isempty(fi))
      %fi = 1:L/2;
      Nplot = round(L*fplot_max/fs);
      fi = 1:Nplot;
  end
  persistent F
  if (isempty(F))
      F =(fi-1)*fs/L;
  end
  
  % These are needed to compute the output, y
  persistent FX

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

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Calculate all of the data
  FX = 2*abs(fft(x))/L;

  n=floor(fs/s.fr);
  B=[1 zeros(1,n) -1];
  A=[1 zeros(1,n) -s.pr^n];
  %W=2*pi/L*(0:L/2);
  W = 2*pi*F/fs;
  EW=exp(1i*W);
  HA=abs(polyval(B,EW)./polyval(A,EW));
  HA(HA < 1e-6) = 1e-6;
  y=filter(B,A,x);
  FY=2*abs(fft(y))/L;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Now create the plots
  dw = 0.05; plotw = 0.25;  
  dh = 0.05; ploth = 0.29; 
  plotrowy = [ploth+4.2*dh  dh] + 3*dh;
  plotcolx = [dw            plotw+2.7*dw   ] + 0.65*dw;
  
  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;
  
  % Create invisible plot for text labels at bottom
  subplot('position', [0.65*dw  0 plotcolx(2)+plotw 3*dh], 'visible', 'off');
  text(0.25, 0.40, 'A note = 440 Hz', ...
                  'fontsize', 16, ...
                  'fontweight', 'bold', ...
                  'horizontalalignment', 'center');
  text(0.75, 0.40, 'B note = 494 Hz', ...
                  'fontsize', 16, ...
                  'fontweight', 'bold', ...
                  'horizontalalignment', 'center');
                  
  % Plot section of original signal
  subplot('position', [plotcolx(1) plotrowy(1) plotw ploth]);
  plot(tplot, x(nplot));
  grid on;
  xlim([s.toff s.toff+s.tdur]*1000);
  xlabel('time (ms)');
  ylabel('Value');
  h = title('Original (x)');
  p = get(h, 'position');
  xl = xlim();
  set(h, 'position', [xl(1) p(2:3)], 'horizontalalignment', 'left');

  % Plot spectrum of original signal
  subplot('position', [plotcolx(2) plotrowy(1) plotw ploth]);
  %semilogy(FX);
  semilogy(F/1000, FX(fi));
  hold on;
  semilogy(F/1000, abs(HA), 'r');
  hold off;
  xlim([0 fplot_max/1000]);
  ylim([1e-4 2]);
  grid on;
  xlabel('Frequency (kHz)');
  ylabel('Value');
  h = title('Original (x)');
  p = get(h, 'position');
  xl = xlim();
  set(h, 'position', [xl(1) p(2:3)], 'horizontalalignment', 'left');
  
  % Plot section of filtered signal
  subplot('position', [plotcolx(1) plotrowy(2) plotw ploth]);
  plot(tplot, y(nplot));
  grid on;
  xlim([s.toff s.toff+s.tdur]*1000);
  xlabel('time (ms)');
  ylabel('Value');
  h = title('Filtered (y)');
  p = get(h, 'position');
  xl = xlim();
  set(h, 'position', [xl(1) p(2:3)], 'horizontalalignment', 'left');

  % Plot spectrum of filtered signal
  subplot('position', [plotcolx(2) plotrowy(2) plotw ploth]);
  %semilogy(FX);
  semilogy(F/1000, FY(fi));
  xlim([0 fplot_max/1000]);
  ylim([1e-4 2]);
  grid on;
  xlabel('Frequency (kHz)');
  ylabel('Value');
  h = title('Filtered (y)');
  p = get(h, 'position');
  xl = xlim();
  set(h, 'position', [xl(1) p(2:3)], 'horizontalalignment', 'left');

  if(s.xplay)
    xplayer = audioplayer(x, fs);
    playblocking(xplayer);
    simState('xplay', false);
  end

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

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