function [SNR_weighted2,SNR] = SNR_weighted2(S,N,fs,clipping);
%function [SNR_weighted2,SNR] = SNR_weighted2(S,N,fs);
%
%calculates a weighted SNR according to the band importance function
%
%INPUTS:
%	S				=	filtered clean speech signal
%	N				=  filtered noise signal
%	fs				=	sampling frequency
%						if fs <= 9 kHz, the 14 band 1/3 octave band 
%					             (ANSI S1.1-1986 standard) values	 are used
%						else the band importance values in the SII standard 
%						(ANSI S3.5-1997) are used
%	clipping		=	if '1' the SNRs are limited to the interval [0,30] dB
%						if '0' the SNRs are not limited to [0,30] dB


%Band importance function
r=2^(1/6);
if fs/2 > 4500;
   I=[83 95 150 289 440 578 653 711 818 844 882 898 868 844 771 527 364 185]*1e-4;
   F=[160 200 250 315 400 500 630 800 1000 1250 1600 2000 2500 3150 4000 5000 6300 8000];
   f2=F*r;
   n=sum(fs/2>f2);
	F=F(1:n);
   I=I(1:n);
   %disp([num2str(n) '/18 banden, max: ' num2str(sum(I)*100) '%']);
else
   I=[128 320 320 447 447 639 639 767 959 1182 1214 1086 1086 757]*1e-4;
   F=[200 250 315 400 500 630 800 1000 1250 1600 2000 2500 3150 4000];
   f2=F*r;
	n=sum(fs/2>f2);
	F=F(1:n);
   I=I(1:n);
   %disp([num2str(n) '/14 banden, max: ' num2str(sum(I)*100) '%']);
end

SNR=zeros(1,n);

%calculation of the SNR values in the different bands
for i=1:n
   [b,a]	 = oct3dsgn(F(i),fs,3);	
   Ep(i)=rmsdb(filter(b,a,S));
   
   %noise spectrum level -> should include the noncorrelated as well as the 
   %correlated noise (i.e. reverberation). Here, only the uncorrelated noise is
   %taken into account (-> only valid in case of small reverberation)!! 
   
   Np(i)=rmsdb(filter(b,a,N));
   SNR(i)=Ep(i)-Np(i);
   if clipping 
      SNR(i)=min(max(0,SNR(i)),30);
   end
end

SNR_weighted2=I/sum(I)*SNR';


%  100% rip-off from simon_doclo (found on the internet)
function y=rmsdb(x)
   % Root Mean Square
%	if x is a vector, y = rms( x ) returns the root mean square
%	value of the elements of x.
%
%	if x is a matrix, y = rms( x ) returns the root mean square
%	value of each column of x.

[r,c] = size( x );
if c == 1,
   y = 10 * log10( x'*x / length(x) );
elseif r ==1,
   y = 10 * log10( x*x' / length(x) );
else
   y = 10 * log10( sum( x .* x ) / length(x) );
end