function theta=mpm(u,d,N)
%N=length(u);
L=floor(N/2)-1;
Y=hankel(u(1:N-L),u(N-L:N));
[U,S,Vs]=svd(Y);
D=abs(real(S));
[Sn,k]=sort(diag(-D));
V=Vs(:,k);
V=V(:,1:d);
V1=V(1:length(V)-1,:);
V2=V(2:length(V),:);
V3=pinv(V1)*V2;
eigy=eig(V3);
f=angle(eigy)/(pi);
if(f<0)
    f=-f;
elseif(f>0)
    f=1-f;
end
theta=asin(f)*180/pi;