function [amp, badmask] = clean_amp(amp, varargin)
% [cleaned_amp, badmask] = clean_amp(amp, varargin);
% in order:
% - account for vmdas bug where bottom amp spikes
% - run a median filter in the horizontal (removes acoustic interference)
% - align so all profiles have the same vertical average (for one noisy beam)
% - return a mask with 1 where amp falls below noise floor or exceeds
%                 acoustic noise threshold (mfilt)
%
% ampfilt                number          (use [] to  disable) 3-point filter
%                                         10  default for AGC
% align                 0 or 1            align input profiles for easier 
%                                                 visual comparison
%                                         0 (default) do not align means
%                                         1 (make all beams have same 
%                                                vertical mean per profile)
% osbug              0 or 1                deal with vmdas bug 
%                                           where bottom amplitude spikes.
%                                             (replace bottom amp with next 
%                                             shallower value)
%                                          0  (default) 1 (enable)
%
% noisefloor       min(amp)                     noise floor per profile
%                                                  1xNprofs or constant
% nedit_threshold   []                make any amp greater than XX into NaNs:
%                                           XX = nedit_theshold + noisefloor
%                                     default:  set to [] to disable
%                                               
%                      


cfg.align = 0;
cfg.ampfilt = 10; %tolerance 
cfg.osbug = 0;
cfg.noisefloor = [];
cfg.nedit_threshold = 7; %roughly 3 dB 

cfg = fillstruct(cfg, varargin, 0);


if length(size(amp)) == 3
   [nbins, nbeams,nprofs]=size(amp);
   if cfg.osbug == 1
      amp(end,:,:) = amp(end-1,:,:);
   end
else
   [nbins, nprofs]=size(amp);
   nbeams = 1;
   if cfg.osbug == 1
      amp(end,:) = amp(end-1,:);
   end
end
badmask = zeros(nbins, nprofs);


vmult = ones(nbins,1);

if nbeams == 1
   cfg.align = 0;
end
   
if cfg.align == 1
   allm = zeros(nbins, 1);
end

%% invoke "median" filter (this block of code is transported from bigstarget)


if (~isempty(cfg.ampfilt))
   if length(cfg.ampfilt) == 2
      fprintf('ERROR in bigstarget: \n')
      fprintf('specify filter tolerance only\n')

      nanmask = [];
      return
   end
   if (nprofs >= 3)
      % Make a square version if it was 3-D.
      if nbeams > 1,
         aa = reshape(amp, nbins * nbeams, nprofs);
      else
         aa=amp;
      end

      % Average of profiles i-i and i+1.
      aamn = 0.5 * ([aa(:, 2:end), aa(:,end)] + ...
                    [aa(:, 1), aa(:,1:(end-1))]);
      % Find large positive outliers.
      badmask = (aa - aamn) > cfg.ampfilt;
      aa(badmask) = aamn(badmask);
      % Make it 3-D again if necessary.
      if nbeams > 1,
         amp = reshape(aa, nbins, nbeams, nprofs);
         badmask = reshape(badmask, nbins, nbeams, nprofs);
      else
         amp = aa;
      end
   end
end


if nbeams > 1
   for beami = 1:nbeams
      eval(sprintf('amp%d = squeeze(amp(:,%d,:));',beami, beami))
   end
else
   amp1 = amp;
end


if (cfg.align == 1) & (nprofs>1)
   for beami = 1:nbeams
      eval(sprintf('allm = allm + mstdgap(amp%d'')'';', beami))
      eval(sprintf('amp%dM =  repmat(mstdgap(amp%d'')'', 1,nprofs);',...
                   beami, beami))
   end
   allm = allm/nbeams;
   for beami = 1:nbeams
      eval(sprintf('amp%d = amp%d-(amp%dM-repmat(allm,1,nprofs));', ...
                   beami,beami, beami))
   end
end


%% noise floor
if isempty(cfg.noisefloor)
   noisefloor = 100000*ones(nbeams,nprofs);
   for beami = 1:nbeams
      for profi = 1:nprofs
         eval(sprintf('noisefloor(%d, %d) = min(amp%d(:,%d));',...
                      beami, profi, beami, profi))
      end
   end
elseif (length(cfg.noisefloor) == 1)
   noisefloor = cfg.noisefloor*ones(nbeams,nprofs);
end

% do not test equality or we'll lose too many in low amp regions
if ~isempty(cfg.nedit_threshold)
   if nbeams == 1
      bmask = zeros(size(amp));
      nfe  = cfg.nedit_threshold + vmult*noisefloor;
      bmask(amp < nfe) = 1;      
   else      
      bmask = [];
      for beami = 1:nbeams
         ca = squeeze(amp(:,beami, :));
         nfe  = cfg.nedit_threshold + vmult*noisefloor(beami,:);
         mask = zeros(size(ca));
         mask(ca < nfe) = 1;      
         bmask = cat(3, bmask, mask);
      end
   end

   if (nbeams > 1)
      bmask = permute(bmask, [1,3,2]);
   end   
   badmask = badmask | bmask;
end




