function FF=samplesompi(x,m)

% samplesompi.m

% Save as 'samplesompi.m' on the same directory with 'testsompi.m'.

% input

% x=[x1 x2 ... xN] :Data vector; line vector

% m : AR order; scalar

% output

% FF=[f g A phai]

% f: frequency vector [f1 f2 .. f2m]'

% g: growing rate [g1 g2 .. g2m]'

% A: initial amplitude [A1 A2 .. A2m]'

% phai: initial phase [phai1 phai2 .. phai2m]'

N=size(x);

N=(N(1,2)-1)/2;

Y=x(2*m+1:2*N+1);

for j=2:2*m+1

Y=[Y;x(2*m+2-j:2*N+2-j)];

end

P=Y*Y'/(2*N-2*m+1);

% solving normal equation (lambda equation)

[V,D]=eig(P);

R=diag(D);

[lambda,I]=min(R);

a=V(:,I);

% to obtain namiso z

z=roots(a');

% f-g plot

g=log(abs(z));

f=angle(z)/2/pi;

[f g];

% to obtain initial amplitude and initial phase

zz=z';

Z=z';

for i=2:2*N+1

zz=zz.*z';

Z=[Z;zz];

end

c=pinv(Z)*x';

A=abs(c);

phai=angle(c);

u=Z*c;

FF=[f g A phai];

% drawing the graphs

figure(1)

plot(f,g,'+')

xlabel('Frequency /Hz');

ylabel('Gradient /sec-1');

grid on

figure(2)

px=[f f]'; py=[g g]';

pz=[zeros(size(A)) A]';

plot3(px,py,pz,'r')

grid on

xlabel('f /Hz'); ylabel('g /sec-1')

zlabel('Initial Amplitude')

% to caluculate model time series

%u=Z*c;

figure(3)

t=1:2*N+1;

plot(t,x,'b-',t,real(u),'r:')

grid on

xlabel('Time /sec');

ylabel('Amplitude');