%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%     esprit   method    %%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
clc;
%%%%%%%%%%%%%         signal        %%%%%%%%%%%%%
t_f=100
t=0:1/t_f:1-1/t_f;
x=randn(1,length(t))+2*sin(2*pi*20.5*t)+10*sin(2*pi*48*t);
%%%%%%%%%%%%%%%     conditions     %%%%%%%%%%%%%%%%
N=length(x);
M=round(N/3);
L=N-M;
K=4;
%%%%%%%%%%%          y matrix         %%%%%%%%%%%%
for i=1:L
    xx=x((i+1):(i+M));
    b=rot90(xx,3);
    y(:,i)=b;
end
%%%%%%%%%%%          R matrix         %%%%%%%%%%%%
R=0;
for i=1:L
   cy=y(:,i)*(y(:,i))';
   R=R+cy;
end
R=R/N;
%%%%%%%%%%%%%%%% eigendecomposition %%%%%%%%%%%%%%%
[v,d]=eig(R);
G=v(:,1:M-K);
S=v(:,M-K+1:M);
GG=G*G';
S1=S(1:M-1,:);
S2=S(2:M,:);
%%%%%%%%%%%%%%%%%       fai       %%%%%%%%%%%%%%%%
fai=inv(S2'*S2)*S2'*S1;
[v,d]=eig(fai);
%%%%%%%%%%%%%%%%%%%%   theta   %%%%%%%%%%%%%%%%%%%%
for i=1:K
    freq(i)=phase(d(i,i))/2/pi;
end
freq
%%%%%%%%%%%%%%%%%%%%    繪圖   %%%%%%%%%%%%%%%%%%%%
subplot(2,1,1);
plot(x);
s=1;
title('signal picture');
xlabel('t');ylabel('x value');
subplot(2,1,2);stem(freq,abs(d),'c',':');
hold on;
subplot(2,1,2);plot(freq,s,'.');
hold off;
title('esprit method');
xlabel('F(hz)');ylabel('ample value')