function [lambda, e] = lsce(h, dt, N)
%LSCE Least Squares Complex Exponential method for EMA
% [LAMBDA, E] = LSCE(H, DT, N) fits N poles using the impulse responses
% in the rows of the matrix H with a time increment DT. The poles are
% returned in the vector LAMBDA, and the least square error is returned
% in E.

% EMA Toolbox
% Peter Fiala, 2008. feb.

Ni = size(h,1); % Number of IRs
Nt = size(h,2); % Number of samples per IR
% assembling the covariance matrix
R = zeros(Ni*(Nt-2*N),2*N); % covariance matrix
H = zeros(Ni*(Nt-2*N),1); % output vector
k = 1;
for i = 1 : Ni
for n = 2*N+1 : Nt
R(k,:) = h(i,n-1:-1:n-2*N);
H(k,:) = -h(i,n);
k = k + 1;
end
end
H(isnan(H)) = 0;
H(isinf(H)) = 0;
a = R \ H; % Solution in LMS sense
z = roots([1 a.']);
lambda = log(z)/dt;
lambda = comconpair(lambda); % keep complex conjugate pairs
lambda = lambda(real(lambda) < 0); % with negative real part
e = norm(R * a - H)^2; % The LS error
end