function mod06_03
  % Simulate noise removal.
  %
  % This Octave/MATLAB script was written to mimic a corresponding
  % Labview module (Module 6.3) 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;

  % Resonater Filter order
  nres = 10;
  % Resonator Fundamental Frequency
  fr = 494;
  % Resonator Decay constant for impulse response
  alpha = 10;
  
  % Butterworth cutoff
  fc = 1000;
  % Butterworth order 
  nbut = 10;
  
  % Threshold
  T = 0.0015;
  
  % Noise level
  Ns = 0.2;
  
  % Save initial simulation state
  simState('toff', toff);
  simState('tdur', tdur); 
  simState('nres', nres);
  simState('fr', fr);
  simState('alpha', alpha);
  simState('fc', fc);
  simState('nbut', nbut);
  simState('T', T);
  simState('Ns', Ns);
  simState('filter', 0);
  
  guiState('xup', true);
  guiState('yup', true);
  guiState('zup', true);
  guiState('hup', true);
  guiState('xplay' , false);
  guiState('yplay' , false);
  guiState('zplay' , false);
  guiState('s4p1', -1);
  guiState('s5p1', -1);
  guiState('s5p2', -1);
  guiState('s7p1', -1);
  guiState('s8p1', -1);

  clf;
  createPanel();
  createPlot();
end

function createPanel
  pos = [0.65 0.05 .30 .9 ];
  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 0.93 1 .07]);
  b = uicontrol(bg, 'style', 'radiobutton', ...
                'units', 'normalized', ...
                'callback', @guiCallback, ...
                'value', 1, ...
                'position', [0.2 0.4 .1 .4]);
  guiState('rbutton', b);
  b = uicontrol(bg, 'style', 'radiobutton', ...
                'units', 'normalized', ...
                'callback', @guiCallback, ...
                'value', 0, ...
                'position', [0.5 0.4 .1 .4]);
  guiState('bbutton', b);
  b = uicontrol(bg, 'style', 'radiobutton', ...
                'units', 'normalized', ...
                'callback', @guiCallback, ...
                'value', 0, ...
                'position', [0.8 0.4 .1 .4]);
  guiState('tbutton', b);

  % Noise Level
  h = sliderControl(p, 'Noise', 0, ...
                      1, simState('Ns'), [0.1 0.2], ...
                      [0.05 .32 0.9 .07], @guiCallback );
  guiState('Nsh', h);
  
  % Play buttons
  bg = uibuttongroup(p, 'position', [.0 .1 1 .15]);
  bx = uicontrol(bg, 'style', 'pushbutton', ...
                'units', 'normalized', ...
                'string', 'x', ...
                'callback', @guiCallback, ...
                'position', [0.25 0.4 .1 .6]);
  by = uicontrol(bg, 'style', 'pushbutton', ...
                'units', 'normalized', ...
                'string', 'y', ...
                'callback', @guiCallback, ...
                'position', [0.5 0.4 .1 .6]);
  bz = uicontrol(bg, 'style', 'pushbutton', ...
                'units', 'normalized', ...
                'string', 'z', ...
                'callback', @guiCallback, ...
                'position', [0.75 0.4 .1 .6]);
  uicontrol(bg, 'style', 'text', ...
                'units', 'normalized', ...
                'string', 'Play', ...
                'position', [0.3  0.1 .5 .3]);
  
  guiState('bx', bx);
  guiState('by', by);
  guiState('bz', bz);

  h = uicontrol(p, 'style', 'text', ...
                'units', 'normalized', ...
                'string', 'ready', ...
                'backgroundcolor', [.2 .2 .2], ...
                'foregroundcolor', [1 1 1], ...
                'position', [0.25  0.0 .5 .1]);
  guiState('status', h);

  % Resonator Panel
  tp = uipanel (p, 'position', [0 .4 1 .53], ...
                'units', 'normalized');
  guiState('rpanel', tp);
  uicontrol(tp, 'style', 'text', 'units', 'normalized',...
           'position', [0.25 0.9 0.5 .1], 'String', 'Resonator');

  h = sliderControl(tp, 'Fund. Freq.', 300, ...
                      1000, simState('fr'), [1 10], ...
                      [0.05 .75 0.9 .12], @guiCallback );
  guiState('frh', h);
  h = sliderControl(tp, 'Alpha', 10, ...
                      100, simState('alpha'), [5 10], ...
                      [0.05 .5 0.9 .12], @guiCallback );
  guiState('alphah', h);
  h = sliderControl(tp, 'Order', 2, ...
                      10, simState('nres'), [1 2], ...
                      [0.05 .25 0.9 .12], @guiCallback );
  guiState('nresh', h);
           
           
  % Butterworth Panel
  tp = uipanel (p, 'position', [0 .4 1 .53], ...
                'units', 'normalized', 'visible', 'off');
  guiState('bpanel', tp);
  uicontrol(tp, 'style', 'text', 'units', 'normalized',...
           'position', [0.25 0.9 0.5 .1], 'String', 'Butterworth');

  h = sliderControl(tp, 'Cut. Freq.', 1000, ...
                      2500, simState('fc'), [50 200], ...
                      [0.05 .75 0.9 .12], @guiCallback );
  guiState('fch', h);
  h = sliderControl(tp, 'Order', 2, ...
                      10, simState('nbut'), [1 2], ...
                      [0.05 .25 0.9 .12], @guiCallback );
  guiState('nbuth', h);

  % Threshold Panel
  tp = uipanel (p, 'position', [0 .4 1 .53], ...
                'units', 'normalized', 'visible', 'off');
  guiState('tpanel', tp);
  uicontrol(tp, 'style', 'text', 'units', 'normalized',...
           'position', [0.25 0.9 0.5 .1], 'String', 'Threshold');
  h = sliderControl(tp, 'Thresh', -3, ...
                      -1, round(log10(simState('T'))/0.1)*0.1, [.1 .2], ...
                      [0.05 .5 0.9 .12], @guiCallback );
  set(get(h, 'userdata'), 'string', num2str(simState('T')));

  guiState('Th', 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.rbutton)
    % Resonator button
    % Selecting a button causes the previously selected button to also respond
    % Ignore the deselection event
    % Switch between three control panels based on button selection
    if(get(s,'value') == 1)
      set(g.bpanel,'visible','off');
      set(g.tpanel,'visible','off');
      set(g.rpanel,'visible','on');
      simState('filter', 0);
      guiState('hup', true);
      guiState('zup', true);
    end
  elseif (s == g.bbutton)
    % Butterworth button
    if(get(s,'value') == 1)
      set(g.rpanel,'visible','off');
      set(g.tpanel,'visible','off');
      set(g.bpanel,'visible','on');
      simState('filter', 1);
      guiState('hup', true);
      guiState('zup', true);
    end
  elseif (s == g.tbutton)
    % Threshold button
    if(get(s,'value') == 1)
      set(g.rpanel,'visible','off');
      set(g.bpanel,'visible','off');
      set(g.tpanel,'visible','on');
      simState('filter', 2);
      guiState('hup', true);
      guiState('zup', true);
    end
  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('hup', true);
    guiState('zup', true);
  elseif (s == g.nresh)
    % resonator n slider  (filter order)
    v = round(get(s, 'value'));
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('nres', v);
    guiState('hup', true);
    guiState('zup', true);
  elseif (s == g.alphah)
    % resonator 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('zup', true);
  elseif (s == g.fch)
    % Butter fc slider (cutoff freq)
    v = round(get(s, 'value'));
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('fc', v);
    guiState('hup', true);
    guiState('zup', true);
  elseif (s == g.nbuth)
    % Butter n slider  (filter order)
    v = round(get(s, 'value'));
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('nbut', v);
    guiState('hup', true);
    guiState('zup', true);
  elseif (s == g.Th)
    % threshold slider  
    v = round(get(s, 'value')/0.1)*0.1;
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(10^v,'%6.4f'));
    simState('T', 10^v);
    guiState('hup', true);
    guiState('zup', true);
  elseif (s == g.Nsh)
    % Noise level slider  
    v = round(get(s, 'value')/0.1)*0.1;
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('Ns', v);
    guiState('yup', true);
    guiState('zup', true);
  elseif (s == g.bx)
    guiState('xplay', true);
  elseif (s == g.by)
    guiState('yplay', true);
  elseif (s == g.bz)
    guiState('zplay', 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('trumpet.csv');
  end
  persistent y 
  if (isempty(y))
      randn('seed', 1);
      %rng(0, 'twister');
      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 the output, z
  persistent FX FY 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-5 1]);
    grid on;
    xlabel('freq (kHz)');
    ylabel('Amplitude');
    title('Original (x)');
  end
  
  if(g.yup)  
    %'y redraw'
    guiState('yup', false);
    y = x + s.Ns * randn(1, L);

    if (~ishandle(g.s4p1))
      subplot(3,3,4);
      lh = plot(tplot, y(nplot));
      guiState('s4p1', lh);
      grid on;
      xlim([s.toff s.toff+s.tdur]*1000);
      xlabel('time (ms)');
      ylabel('Amplitude');
      title('Noisy (y)');
    else
      % Line has been previously drawn, update ydata
      set(g.s4p1, 'ydata', y(nplot));
    end

    FY = fft(y);
    FYI = 2*abs(FY(fi))/L;
    if (~ishandle(g.s5p1))
      subplot(3, 3, 5)
      lh = semilogy(F/1000, FYI);
      guiState('s5p1', lh);
      xlim([0 5]);
      ylim([1e-5 1]);
      grid on;
      xlabel('freq (kHz)');
      ylabel('Amplitude');
      title('Noisy (y)');
    else
      % Line has been previously drawn, update ydata
      set(g.s5p1, 'ydata', FYI);
    end
  end

  if(g.hup)  
    %'y redraw'
    guiState('hup', false);
    
    if(s.filter == 0)
      % Resonator
      %%{
      % Method 3
      % This method creates the frequency response by multiplying 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
      FH = ones(size(t));
      for k = 1:s.nres
        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
      % FH is multi-bandstop, subtract from 1 to convert to multi-bandpass
      FH = 1 - FH;
      %h = real(ifft(FH))*fs;
      FHI=abs(FH(fi));
      %toc
      %}
    elseif (s.filter == 1)
      % Butterworth
      odd = (mod(s.nbut, 2) == 1);
      if (odd)
        % # of complex pole pairs
        nc = (s.nbut-1)/2;
      else
          nc = (s.nbut)/2;
      end
      
      wc = 2*pi*s.fc;
      FH = ones(size(t));
      w = (0:(length(FH)-1))/length(FH)*2*pi*fs;
      for k = 1:nc
        phi_k = (pi/s.nbut)*(k - s.nbut/2 - 0.5);
        HF = (wc^2)./((1i*w).^2 + 2*cos(phi_k)*1i*w*wc + wc^2);
        FH = FH .* HF;
      end
      if (odd)
        HF = wc./(1i*w + wc);
        FH = FH.*HF;
      end
      nf = length(FH);
      % The 
      FH = [FH(1:(nf/2+1)) conj(fliplr(FH(2:(nf/2))))];
      FHI=abs(FH(fi));
    elseif (s.filter == 2)
      % Threshhold
      FHI = s.T*ones(size(F));
    end

    if (~ishandle(g.s5p2))
      subplot(3, 3, 5)
      hold on;
      lh = semilogy(F/1000, FHI);
      guiState('s5p2', lh);
      hold off
    else
      % Line has been previously drawn, update ydata
      set(g.s5p2, 'ydata', FHI);
    end
  end

  %%{
  if(g.zup)
    guiState('zup', false);
    
    if(s.filter < 2)
      z = real(ifft(FY.*FH));
      FZ = fft(z);
      FZI = 2*abs(FZ(fi))/L;
    else
      FZ = FY;
      FZ(2*abs(FY)<s.T*L)=0;
      z=real(ifft(FZ));
      FZI = 2*abs(FZ(fi))/L;
      FZI(FZI <= 0) = 1e-6;
    end
    if(~ishandle(g.s7p1))
      subplot(3,3,7);
      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

    if(~ishandle(g.s8p1))
      subplot(3,3,8);
      lh = semilogy(F/1000, FZI);
      guiState('s8p1', lh);
      xlim([0 5]);
      ylim([1e-5 1]);
      grid on;
      xlabel('freq (kHz)');
      ylabel('Amplitude');
      title('Filtered (z)');
    else
      set(g.s8p1, 'ydata', FZI);
    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