clear;
t_f=100;
t=0:1/t_f:1-1/t_f;
x=2*sin(2*pi*10*t)+sin(2*pi*20*t)+randn(size(t));
N=length(x);
M=round(N/3);   
K=2;
L=N-M;
%%%%%%%%%%%求R%%%%%%%%%%%%
R=0; 
for i=1:L
    xx=x((i):(i+M-1));
    b=rot90(xx,3);
    cy=b*b';
    R=R+cy;
end
R=R/(N-M+1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w1=-pi:2*pi/N:pi-2*pi/N;
for n=1:N
   for k=1:M
      a(k)=exp(-j*(k-1)*w1(n));
   end
   c=a.';
   d=c;
   bb(n)=c'*inv(R)*c;
end
cc=M./(bb+conj(bb))*2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dd=cc;
for i=1:K
    [maxroot,jj]=max(abs(dd));
     freq(i)=-1/2+(jj-1)/(2*N);
    dd(jj)=0;
end
q=freq*t_f;
%%%%繪圖%%%
w=t_f*(-(2*N)/4:(2*N)/4-1)/(N);
subplot(2,1,1);plot(x);
title('signal picture');
xlabel('t');ylabel('x value');
subplot(2,1,2);plot(w,abs(cc));
title('Min-norm spectral method picture');
xlabel('w(hz)');ylabel('pbt value');