function [Pnm_at_loc, loc] = noise_masker(P, Ptm, low, high)
% [power_spectral_density_of_noise_at_loc, location_within_critical_band]
%    = noise_maskers(power_spectral_density, power_spectral_density_tones,
%                    low_bin_of_critical_band, high_bin_of_critical_band)
%
% noise_masker takes a power spectrum and related tone spectrum,
% as well as the low and high ranges to deal with (so that you can
% later specify a critical bandwidth), and returns a noise masker
% value along with where it should be located.
%
% If loc=-1, there was NO noise masker in the region.

% this array will hold whether a given element is part of the
% noise.  Begin with ones.
noise_members=ones(1,high-low+1);

% go through tonal power spectrum
for k=low:high,
   % if there is a tone
   if (Ptm(k)>0),
      % check frequency location and determine neighborhood length
      if ((k>2) & (k<63))
         m=2;
      elseif((k>=63) & (k<127))
         m=3;
      elseif((k>=127) & (k<256))
         m=6;
      else
         m=0;
      end
      % set all members of the neighborhood to 0, which
      % removes them from the list of noise members
      for n=(k-low+1)-m:(k-low+1)+m,
         if (n > 0)
            noise_members(n)=0;
         end
      end
   end
end

% those indices in noise_members that have a 1
% are part of the noise masker.

% go through the noise locations and perform power summation,
% which requires going back to magnitude (dB doesn't add linearly.).

% if there are no noise members in the range, then get out
if (isempty(find(noise_members)))
   Pnm_at_loc=0;
   loc=-1;
else
   temp=0;
	for k=(low+find(noise_members)-1),
      temp = temp + dbinv(P(k));
	end

	% take summation and make it into dB
	Pnm_at_loc=10*log10(temp);

   loc=geomean(low+find(noise_members)-1);
end