function thres = findThreshold(y, FFTlength, f, nbits)

% find the frequencies in barks
b=hz2bark(f);

powernormalizationconstant = 90.302;

% calculate the absolute threshold of hearing
ATH=3.64.*(f./1000).^(-.8) - 6.5.*exp(-0.6.*(f./1000-3.3).^2) + 0.001.*(f./1000).^4;

% step 1: normalize the signal
%y=normalize(y, FFTlength, nbits);

% step 2: get power density spectrum
power=psd(hanning(length(y))'.*y, FFTlength, powernormalizationconstant);

% step 3: find tones
tone_psd=find_tones(power(1:FFTlength/2));

% step 4: find noise maskers within critical band
% start with 0-1
noise_psd=zeros(1,length(tone_psd));
lowbin=1;
highbin=max(find(b<1));

% remainder of critical bands can be done with loop
for band = 1:23,
   [noise_psd_at_loc, loc]=noise_masker(power, tone_psd, lowbin, highbin);
   if (loc ~= -1)
      noise_psd(floor(loc))=noise_psd_at_loc;
   end
   lowbin=highbin;
   highbin=max(find(b<(band+1)));
end

% step 5: remove maskers that are near each other
[tone_psd, noise_psd]=check_maskers(tone_psd, noise_psd, ATH, b);

% step 6: now that we have important maskers, calculate global threshold
thres=global_threshold(ATH, tone_psd, noise_psd, b);