function  [L] = spl(S, Sn, f)

% L = spl(S, Sn, f)
%
% spl calculates the sound pressure level of the sample 
%                
%
%   in: S signal spectrum
% Sn  noise spectrum
% f frequency scale (default 0:(nf-1)/25000)
%
%  out: L sound pressure level

%   edited from flimit.m by Esa Ojala 12.8.1991

[n1, nf] = size(S);
if nargin<3; fs = 25000; f = (0:nf-1)/(nf-1)/2*fs; end
%
%  lim = f(max(find(S>p*max(S) & S>Sn )));
%
Ind = find(dbs(S)>-20 & S>Sn);
lind=length(Ind);
%  Ind= find(f(Ind)<6500);
%  Inds = sort(Ind);
%  if lind<=0; L=-50;end
%  if lind>=0;
Ptot=0;
for n = 3:lind
  P(n)=dbs(S(Ind(n)));
  P(n)=P(n)/10;
  Ptot=Ptot+10^P(n);
end
%Ptot
L=10*log10(Ptot);