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)) >= 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.

• You should take Real[Ln(z)], not Real(z). Ln() being natural logarithm. – user24568 Nov 1 '16 at 22:38

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.

• You have considered exponential equation. The question is about complex exponentials (exponentially decaying sinusoids). Infact, Matrix Method and Prony Analysis are two of the most preferred methods to find parameters of such signals. – avr Jan 29 '16 at 12:50

With more sample points, fs = 5000, Scipy least_squares is pretty good with 5 % multiplicative noise, and not bad even with 10 %: Params: tmax 0.11  fs 5000  *noise 5 % laplace
a [0.42 0.99 0.92 0.65] = A + [0.024 -0.013 0.028 -0.0028]
d [71 49 92 78]        = D + [1.4 -0.72 2.3 -2.1]
f [80 120 250 560]     = F + [0.096 -0.068 -0.22 0.0064]

Params: tmax 0.11  fs 5000  *noise 10 % laplace
a [1.2 0.58 0.39 0.66] = A + [0.81 -0.42 -0.5 0.0059]
d [64 115 78 77]       = D + [-5.7 65 -12 -2.7]
f [119 248 251 560]    = F + [39 128 0.59 0.013]

Notes:

As initial x0 I just took the average A, average D, average F. Don't know how sensitive the results are to x0, tmax, fs, noise.

What errors do we want to minimize here -- relative, or absolute ? Early t / big y or later t / tiny y ? Thorny questions, beyond what the OP asked. Huber loss, a "dead band" around errors near 0, looks reasonable, and scipy least_squares can do that; see robust regression in the scipy cookbook.

A trick to avoid multiple minima in optimizing $f()$ with permutable arguments, e.g. $f(x,y,z) = f(z,y,x) = ...$, is to constrain $x \le y \le z$. Here, I didn't bother; A D F come out permuted.

Harmonic inversion aka the filter diagonalization method, FDM, looks powerful. But the code is ~ 1000 lines of dense C, and Mandelshtam's papers are over my head. If anyone knows of other programs or papers on harminv / FDM, or real test cases, please let me know.