function QMask = ThreshHearing(Fs,M)
%THESHHEARING computes the thresholds in human hearing (in db) for pure tones
%varaying in frequency 
%       QMask = ThreshHearing(Fs,M), Fs is the sampling frequency, M number
%       of zone, QMask threshold of hearing per zone

delta = 20;   % minimum frequency resolution in Hz
% Vector of frequencies at which the threshold of hearing is computed
f = delta:delta:Fs/2;
ht = round(5+3.64*(f/1000).^(-0.8)-6.5*exp(-0.6*(f/1000-3.3).*(f/1000-3.3))+0.001*(f/1000).*(f/1000).*(f/1000).*(f/1000));   

% Compute the threshold of hearing per zone

for i=1:M,
     j= find((f<=(i*Fs/(2*M)))&(f>((i-1)*Fs/(2*M))));
     if length(j),
         QMask(i)=min(ht(j));
     else
         QMask(i) = QMask(i-1)+20;
     end
end