function [m, mf] = freqMask(f,A,Fs,delta)
%FREQMASK computes the masking created by frequency f with amplitude A (in db) 
%     [m, mf] = freqMask(f,A,Fs,delta), f is the frequency in Hz, A is the SPL 
%     (Sound Pressure Level) of frequency f in db, Fs is the sampling rate of the input signal 
%     and delta is the frequency resolution (delta > 1), A>0, f>1, Fs>1000.
%     mf are the masked frequencies and m is the corresponding masked level
%     in db
 
b = freq2bark(f);   % converts frequency to bark units where masking effects are 
                    % more easily modeled

% Compute masking effect with appropriate resolution
bw = 0:delta/100:floor(freq2bark(Fs/2));

lambda = 27;        % masking towards lower frequencies is modeled as a line with slope
                    % 27db
    m = A-lambda.*(b-bw);
i=find(bw>b);

lambda = -lambda./(1+(b^0.25)/4)+max(0,A-40)*0.37;
m(i) = A-lambda.*(b-bw(i));
m(m<0)=0;
m = round(m);
mf = bark2freq(bw);
i = find((mf<(f+0.25*f))&(mf>(f-0.25*f)));
m(i)=0;