function nsegsrr = nsegsrr(sd,shat,flen,P)
% NSEGSNR - Computes the segmental SNR with direct-path normalization
%   
%   Usage: 
%       [nsegsrr] = nsegsrr(s,x,flen);
%   
%   Inputs:
%       sd          reference signal (e.g. clean speech)
%       shat        signal under evaluation (e.g. dereverberated speech)
%       flen        frame length (smples)
%
%   Outputs:
%       nsegsrr     normalized segmental SRR
%
%   Reference:
%       [1] P.A. Naylor, N.D. Gaubitch and E.A.P Habets,
%           'Signal-based performance evaluation of dereverberation algorithms',
%           Journal of Electrical and Computer Engineering
%           (united with Research Letters in Signal Processing),
%           Volume 2010, 2010.

%
%
%**************************************************************************
% Author: E.A.P. Habets and N.D. Gaubitch
% Version: $Id: ltassbci.m 269 2011-05-31 09:15:01Z ndg $
%**************************************************************************

%% Ensure both signals are of equal length
sd = sd(:);
shat = shat(:);
if (size(sd,1)~=size(shat,1))
    len = min(sd,shat);
    sd = sd(1:len);
    shat = shat(1:len);
else
    len = length(sd);
end

%% Check input parameters
if nargin < 3
    flen = min(3*8000,len); % Assuming sample frequency is 8000 Hz;
end
if nargin < 4
    P = 10; % LP order
end

%% Whiten signals
[arsd,~,k]  = lpcauto(sd,P,len);
esd = lpcifilt(sd,arsd,k(:,1));

[arshat,~,k]  = lpcauto(shat,P,len);
eshat = lpcifilt(shat,arshat,k(:,1));

%% Compute segsrr using whitened signals
finc = round(.25*flen);
sdf = enframe(sd,flen,finc);
esdf = enframe(esd,flen,finc);
shatf = enframe(shat,flen,finc);
eshatf = enframe(eshat,flen,finc);

% Direct-path normalization per frame
ahatf = sum(sdf.*shatf,2)./sum(sdf*sdf.').';
eshatfdpn = eshatf./repmat(ahatf,1,flen);

% Compute srr per frame
srr_frame = sum(esdf.^2,2)./(sum((eshatfdpn-esdf).^2,2) + eps);  

% Lower-bound on the SRR per frame
low_bound = 10^(max(10*log10(srr_frame)-50)/10);
wrong = find(srr_frame < low_bound); 
srr_frame(wrong) = ones(length(wrong),1)*low_bound;

% Compute mean nsegsrr in dB
nsegsrr = 10*mean(log10(srr_frame));