function [mask, F] = SimMask(x,Fs,M)
%SIMMASK computes the overall masking created by the frequencies of signal x in 
%a set of M zones 
%     [mask, F] = SimMask(x,Fs,M), x is the original signal, Fs is the sampling 
%     rate of the input signal, M is the number of zones in which the 
%     spectral masking is computed
%     mask is the overall mask created in the M zones
 
delta = 20;   % minimum frequency resolution in Hz

L = 1024;
% Compute masking effect with appropriate resolution
W = fft(x,L);    % FFT 1024 points
W=W(1:L/2);
W=abs(W)/L;

SPL = round(120+20*log10(W)); % Sound Pressure Level

spl = reshape(SPL,16,L/32);
mf = [0:L/2-1]*Fs/L;        % We consider that masking effects create only the frequencies 
j = find(max(spl)>60)*16-9;             % with power higher than 60db

le=length(j);
[mask,F] = freqMask(mf(j(1)),SPL(j(1)),Fs,delta);
% mask(j(1)-3:j(1)+3)=0;
for i=2:le,
    [m,F1] = freqMask(mf(j(i)),SPL(j(i)),Fs,delta);
  %  m(j(i)-3:j(i)+3)=0;
    mask = [mask; m];
end
mask=sum(mask);
Mask = zeros(1,M);
K = round(3*M/4);
for i = 1:K,
    j = find((F<(i*Fs/(2*M)))&(F>((i-1)*Fs/(2*M))));
    Mask(i) = min(mask(i));
end
mask = Mask;