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.