function out = dist_intel2(speech_out,speech_in,fs,N);
%
% function out = dist_intel2(speech_out,speech_in,fs,N);
%
% calculates the intelligibility weighted spectral distortion for stationary signals
% INPUTS: speech_out = output speech signal
%         speech_in  = input speech signal
%         fs         = sampling frequency
%         N          = number of frequency points 
% OUTPUTS: out       = intelligibility weighted distortion

%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];
   f1=F/r;
   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;
   f1=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

PSpeech_out = psd(speech_out,2*N);
PSpeech_in  = psd(speech_in,2*N);
freq        = (fs/2)*(0:1:N)/N;

Edist=zeros(n,1);
nrband=1;
f = 0;

while nrband <= n,
  teller=0;
  distP=0; 
  tf2 = f2(nrband); 
  f = max(f,ceil(f1(nrband)*N/(fs/2)+1));
  while freq(f) < tf2,
    teller = teller+1;
    distP = distP + abs(10*log10(PSpeech_out(f)./PSpeech_in(f))); 
    f = f+1;
  end  
  Edist(nrband) = distP/teller;
  nrband = nrband+1;
end

out = I/sum(I)*Edist;