function mod08_07
  % Denoising a Periodic Signal by Thresholding
  %
  % This Octave/MATLAB script was written to mimic a corresponding
  % Labview module (Module 8.7) 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/>.
  %
  
  % Noise Level
  Ni = 5;
  
  % Threshold
  Thresh = 0.9;
  
  % Save initial simulation state
  simState('Ni', Ni);
  simState('Thresh', Thresh);
  simState('xup', true);
  simState('SNRf', '%6.3f');
  
  clf;
  createPanel();
  createPlot();
end

function createPanel
  pos = [0.71 0.02 .28 .96 ];
  %cpanel = uipanel ('position', pos, 'visible', 'off');
  cpanel = uipanel ('position', pos);
  simState('cpanel', cpanel);
  
  %%{
  %% 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, cpanel, pos});
  end
  %}

  h = sliderControl(cpanel, 'noise level', [0.1, 5], ...
                      simState('Ni'), [0.1 0.5], ...
                      [0.05 .55 0.9 .2], @guiCallback );
  simState('Nih', h);

  h = sliderControl(cpanel, 'threshold', [0.01, 1], ...
                      simState('Thresh'), [0.02 0.1], ...
                      [0.05 .30  .9 .2], @guiCallback );
  simState('Threshh', h);

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

  uicontrol(cpanel, 'style', 'text', ...
                'units', 'normalized', ...
                'string', 'SNR (dB)', ...
                'position', [0.25  0.85 .5 .05]);
  SNRdB = 10*log10(59/simState('Ni')^2);
  h = uicontrol(cpanel, 'style', 'text', ...
                'units', 'normalized', ...
                'string', num2str(SNRdB, simState('SNRf')), ...
                'horizontalalignment', 'center', ...
                'backgroundcolor', [0.7 0.7 0.7], ...
                'position', [0.25  0.8 .5 .05]);
  simState('SNRh', h);

  set(cpanel, 'visible', 'on');
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.Nih)
    % Ni slider (noise level)
    v = round(get(s, 'value')/0.1)*0.1;
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    SNRdB = 10*log10(59/v^2);
    set(g.SNRh, 'string', num2str(SNRdB, g.SNRf));
    simState('Ni', v);
    simState('xup', true);
  elseif (s == g.Threshh)
    % Threshold slider  
    v = round(get(s, 'value')/0.01)*0.01;
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('Thresh', v);
  end
  
  createPlot();
end

function createPlot()

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

  persistent x;
  persistent t;
  persistent X;
  persistent XMAG;
  persistent yl1;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Calculate the data

  fs = 1000; 
  
  % Only update x if the noise level changes
  if(s.xup)
    x_t = @(t) 1 + 4*sin(10*pi*t) + 10* cos(20*pi*t + 0.6435);
    t = (0:(fs-1))/fs;
    x = x_t(t) + s.Ni*randn(size(t));
    X = fft(x)/fs;
    XMAG = abs(X);
  end

  Y = X;
  Y(XMAG < s.Thresh) = 0;
  y = real(ifft(Y*fs));
  YMAG = abs(Y);


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Now create the plots
  dw = 0.05; plotw = 0.26;  
  dh = 0.05; ploth = 0.33; 
  plotrowy = [ploth+4*dh  dh] + 1.1*dh;
  plotcolx = [dw            plotw+2.6*dw] + 0.65*dw;
  
  % Plot the original signal
  if(s.xup)
    subplot('position', [plotcolx(1) plotrowy(1) plotw ploth]);
    plot(t, x);
    yl1 = ylim();
    grid on;
    xlabel('time (s)');
    ylabel('Value');
    title('Noisy Signal (time)');

    simState('xup', false);
  end

  % Plot FFT of the noisy signal
  fi = 1:21;
  F = (fi - 1)*fs/length(X);
  subplot('position', [plotcolx(2) plotrowy(1) plotw ploth]);
  h = semilogy(F, XMAG(fi), 'o');
  % Plot the stems
  lcolor = get(h, 'color');
  yl2 = ylim();
  hold on
  for k = 1:length(F)
    line([F(k) F(k)], [yl2(1) XMAG(k)], 'color', lcolor);
  end
  % Plot the threshold line
  xl = xlim();
  line([xl(1) xl(2)], [s.Thresh s.Thresh], 'color', 'r'); 
  hold off;
  grid on;
  xlabel('frequency (Hz)');
  ylabel('Value');
  title('Filtered Signal (freq)');

  % Plot the filtered signal
  subplot('position', [plotcolx(1) plotrowy(2) plotw ploth]);
  plot(t, y);
  ylim(yl1);
  grid on;
  xlabel('time (s)');
  ylabel('Value');
  title('Noisy Signal (time)');
  
  % Plot FFT of the filtered signal
  subplot('position', [plotcolx(2) plotrowy(2) plotw ploth]);
  % If any YMAG values are zero we get an error in the log plot
  % Clip the minimum value.
  YMAG(YMAG < yl2(1)) = yl2(1)/100;
  h = semilogy(F, YMAG(fi), 'o');
  % Plot the stems
  lcolor = get(h, 'color');
  ylim(yl2);
  hold on
  for k = 1:length(F)
    line([F(k) F(k)], [yl2(1) YMAG(k)], 'color', lcolor);
  end
  hold off;
  grid on;
  xlabel('frequency (Hz)');
  ylabel('Value');
  title('Noisy Signal (freq)');

  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
