function lines=music(y,K,fs) %y: signal, K: #spectral lines, fs: sampling frequency

N=5*K+1; % ??
Ryy=mycorr(y,N); %N: window size

[V,D]=eig(Ryy);
[d,ind]=sort(diag(D), 'descend'); %sort the eigenvalues
Ds=D(ind,ind);
Vs=V(:,ind); %reordered eigenvectors

G_N_NK=Vs(:,[K+1,length(Ryy)]); %take the N-K last (sorted) eigenvectors

omegas=-pi:2*pi/100:pi; %the spectral lines are in the interval [-pi,pi]

A=zeros(1,length(omegas));
for k=1:length(omegas)
   eN1w=[exp(-1i*[0:length(Ryy)-1]*omegas(k))].'; % .': tranpose
   valw=(eN1w')*(G_N_NK)*(G_N_NK')*(eN1w); % ': Hermitian
   A(k)=1/valw; %easier to find max than min
end

maxes_idx=findmax(A,K);
lines=omegas(maxes_idx);
lines=lines.*(fs/(2*pi));

end
