Matrix Pencil method

[bo_fr...@my-deja.com]

I want to approximate a time series with a sum of complex exponentials,
i. e. finding the poles of the system whos output I am analysing.

I have red about the matrix pencil method [1] and it seems promising.
This reference actually includes some computer code.
[1] "Using the Matrix Pencil Method to Estimate the Parameters of a Sum
of Complex Exponentials", T K Sarkar and O Pereira, IEEE Antennas and
Propagation Magazine, vol 37, no 1, Feb 1995.

Q: Is there anyone out there having implemented the Matrix Pencil
method in Matlab code?

Sent via Deja.com http://www.deja.com/
Before you buy.

[Matti Taskinen]

Try this:

function [b, a, f] = mpencil(y, x, tol, L)
% function [B, A, F] = mpencil(Y, X, [TOL, L])
%MPENCIL Fits a sum of exponential functions to data by a Matrix Pencil method
%
% B = MPENCIL(Y,X) approximates Y(i) by a weighted sum of
% exponentials sum_j=1^p A(j)*exp(B(j)*X(i)). The X(i) are
% assumed equally spaced.
%
% The number of exponential is defined with the third optional
% parameter TOL. If TOL < 1 (and > 0) the number of terms is the same
% as the number of such singular values (of the data) that are larger
% than TOL*biggest singular value. If TOL >= 1 and an integer, the
% number of exponential terms is TOL. The default TOL value is 1e-3.
%
% The last optional parameter L is a parameter for the Matrix Pencil
% method. Values between N/3 and N/2 are preferred where N is the
% number of data (Y).
% The default L is the ceil(1/2*(ceil(N/3)+floor(N/2))).
%
% The A(i) and the fitted values F(i) are optionally returned by
% [B,A,F] = MPENCIL(Y,X).

% References:
% Sarkar, T.K. and Pereira, O.
% Using the Matrix Pencil Method to Estimate the Parameters of a Sum
% of Complex Exponentials.
% IEEE Antennas and Propagation Magazine, Vol. 37, No 1, Feb. 1995

% Matti Taskinen, Rolf Nevanlinna Institute, University of Helsinki
% 06.03.1996. Latest revision 01.04.1996.

% Check the vector sizes:
N = size(y, 1);
if size(x, 1) ~= N
error('length(Y) should be length(X)');
end

if size(y, 2) ~= 1
error('Y should be column vector');
end

if size(x, 2) ~= 1
error('X should be column vector');
end

% Fill in the missing parameters:
if nargin < 4
L1 = ceil(1/3 * N);
L2 = floor(1/2 * N);
L = ceil((L1 + L2) / 2);
if nargin < 3
tol = 1e-3;
end
end

% Check the parameter values:
M_given = round(tol) == tol;
if (tol < 0) | ((~M_given) & (tol > 1))
error('TOL should be either >= 1 and an integer or < 1 and > 0');
end

if L > N/2
error('L shoud be < N/2');
end

if M_given & (L < tol)
error('TOL should be <= L');
end

% X are assumed to be equally spaced:
T = diff(x(1:2));

Y = zeros(N-L, L+1);
ind = 0:N-L-1;
for j = 1:L+1
Y(:,j) = y(ind+j);
end

[U,S,V] = svd(Y,0); % Economy size!

if M_given
% The number of exponential terms may be given in TOL:
M = tol;
else
% Otherwise figure it out from the singular values:
D = diag(S);
for M = 1:length(D)-1
if (abs(D(M+1)/D(1)) <= tol)
break;
end
end
end

SM = S(:,1:M);

VM = V(:,1:M);
V1 = VM(1:L,:);
V2 = VM(2:L+1,:);

Y1 = U*SM*V1';
Y2 = U*SM*V2';

A = pinv(Y1)*Y2;
z = eig(A);
z = z(1:M);
b = 1/T * log(z);

if nargout > 1
Z = exp(x*b.');
a = Z\y;

if nargout > 2
f = Z*a;
end
end

[Garrett Barter]

I can attest that this code works well. One point of improvement is to sort the eigenvalues:
z = sort(eig(A), 'descend');
instead of,
z = eig(A);

[sandeepb...@gmail.com]

Can you describe what is x, y etc. What the variables represent.

[nagar...@gmail.com]

On Tuesday, December 12, 2017 at 10:38:14 PM UTC-8, sandeepb...@gmail.com wrote:
> Can you describe what is x, y etc. What the variables represent.

Dear Sandeep sir,

This only for you if you are willing to transfer $1000 to my bank account then I will provide the code for matrix pencil method which is applicable to your research.

[Hm...@mail.aub.edu]

does anyone have a more clear code plz