function [threshold, start] = mask_threshold(type, j, P, b)
% [masking_threshold, starting_bin] = mask_threshold(type, location_of_masker_bin,
%                                     power_spectral_density_at_j, bark_spectrum)
%
% mask_threshold returns an array of the masked threshold in dB SPL that
% results around a mask located at a frequency bin (i.e., in 
% discrete terms).  It also returns a starting index for this threshold,
% which is discussed later.
%
% The user should also supply the power spectral density and the related Bark
% spectrum so that all calculations can be made. Note also that two
% different threshold are possible, so the user should specify:
%
%  type = 0      threshold = NOISE threshold
%  type = 1      threshold = TONE  threshold
%
% This thresholding is determined in a range from -3 to +8 Barks
% from the mask. (This is why a bark spectrum is needed.)  In case you
% would like to overlay different thresholds, you need to know where each
% one actually starts.  Thus, the starting bin for the threshold is also
% returned.

% determine where masker is in barks
maskerloc=b(j);

% set up range of the resulting function
% in barks
low=maskerloc-3;
high=maskerloc+8;
% in discrete bins
lowbin=max(find(b<low));
if (isempty(lowbin))
   lowbin = 1;
end
highbin=max(find(b<high));

% calculate spreading function
SF = spreading_function(j, P, lowbin, highbin, b);

% which type are we dealing with?  The equations are slightly different.
if (type==0)
   % calculate noise threshold
   threshold=P-.175*b(j)+SF-2.025;
else
   % calculate tone threshold
   threshold=P-.275*b(j)+SF-6.025;
end

% finally, note that the lowest value in threshold corresponds
% to the frequency bin at lowbin.
start=lowbin;

% -----------------------------------------------------------------
function spread = spreading_function(masker_bin, power, low, high, b)
% spreading_around_masker = spreading_function(masker_bin, masker_power,
%                                              low_bin, high_bin, bark_spectrum)
%
% spreading_function takes a masker_bin and related power to determine
% the spreading function that surrounds it.  This function indicates
% the minimum level that nearby frequencies must attain before they
% are detectable by the human ear.  Along with the power in the 
% masker and other data, this can be used to determine whether a 
% signal is audible over a masker.

masker_bark=b(masker_bin);

for i=low:high,
   maskee_bark=b(i);
   deltaz=maskee_bark-masker_bark;
   
   if ((deltaz>=-3.5) & (deltaz<-1))
      spread(i-low+1)=17*deltaz-0.4*power+11;
   elseif ((deltaz>=-1) & (deltaz<0))
      spread(i-low+1)=(0.4*power+6)*deltaz;
   elseif ((deltaz>=0) & (deltaz<1))
      spread(i-low+1)=-17*deltaz;
   elseif ((deltaz>=1) & (deltaz<8.5))
      spread(i-low+1)=(0.15*power-17)*deltaz-0.15*power;
   end
end