Finding parameters of exponentially decaying sinusoids (Matrix Pencil Method) The matrix pencil method is an algorithm which can be used to find the individual exponential decaying sinusoids' parameters (frequency, amplitude, decay factor and initial phase) in a signal consisting of multiple such signals added. I am trying to implement the algorithm. The algorithm can be found in the paper from this link or this link.

In order to test the algorithm, I created a synthetic signal composed of four exponentially decaying sinusoids generated as follows:

fs=2205;
t=0:1/fs:249/fs;
f(1)=80;
f(2)=120;
f(3)=250;
f(4)=560;
a(1)=.4;
a(2)=1;
a(3)=0.89;
a(4)=.65;
d(1)=70;
d(2)=50;
d(3)=90;
d(4)=80;
for i=1:4
    x(i,:)=a(i)*exp(-d(i)*t).*cos(2*pi*f(i)*t);
end
y=x(1,:)+x(2,:)+x(3,:)+x(4,:);


I then feed this signal to the algorithm described in the paper as follows: 

function [f d] = mpencil(y)

%construct hankel matrix
N = size(y,2);
L1 = ceil(1/3 * N);
L2 = floor(2/3 * N);
L = ceil((L1 + L2) / 2);

fs=2205;
for i=1:1:(N-L)
    Y(i,:)=y(i:(i+L));
end

Y1=Y(:,1:L);
Y2=Y(:,2:(L+1));

[U,S,V] = svd(Y);
D=diag(S);
tol=1e-3;
m=0;
l=length(D);
for i=1:l
    if( abs(D(i)/D(1)) &gt;= tol)
        m=m+1;
    end
end
Ss=S(:,1:m);
Vnew=V(:,1:m);
a=size(Vnew,1);
Vs1=Vnew(1:(a-1),:);
Vs2=Vnew(2:end,:);
Y1=U*Ss*(Vs1');
Y2=U*Ss*(Vs2');
D_fil=(pinv(Y1))*Y2;
z = eig(D_fil);

l=length(z);
for i=1:2:l
    f((i+1)/2)= (angle(z(i))*fs)/(2*pi);
    d((i+1)/2)=-real(z(i))*fs;
end


In the output from the above code, I am correctly getting the four constituent frequency components but am not getting their decaying factors. If anybody has prior experience with this algorithm or has some understanding about why this discrepancy might be there, I would be very grateful for your help. I have tried rewriting the code from a scratch multiple times but it has been of no help, giving the same results. 

Any help would be highly appreciated.
In his book, Numerical Methods that Work, Forman S. Acton writes, page 253:


  For it is  well known that an exponential equation of this type [$y=Ae^{-at}+Be^{-bt}$] in which all four parameters are to be fitted is extremely ill conditionned. That is, there are many combinations of $(a,b,A,B)$ that will fit most exact data quite well indeed (will you believe four significant figures? and when experimental noise is thrown into into the pot, the entire operation becomes hopeless. But those with Faith in Science...


A four exponential sum seems quite ill-posed to me, yet there have been recent advances using sparsity assumptions or more involved optimization algorithms. Recent examples can be found for instance in Parameter estimation of monomial-exponential sums or Nonlinear Approximation by Sums of Exponentials and Translates, or citations thereof (I have not had the opportunity to test them by myself.