function [a,Rx,fLS]=ARMA_LS(M,x,p,varargin)

%p=4;
[Q,N]=size(x);
%M=100;
%n=1:N;
%sampling=100;
%fls=zeros(sampling,2);
%for k=1:sampling
%     w=randn(size(n));
%     x = sqrt(20)*sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+w;
    Rx=xcorr(x,varargin{:});
    R=zeros(M,p);
    for i=1:M
        R(i,:)=Rx(p+i-1+N:-1:i+N);
    end
    b=-Rx(p+1+N:p+M+N)';
    a=(inv(R'*R))*R'*b;
    
    a=[1;a];   %add a0=1
    z=roots(a');
    f=atan(imag(z)./real(z))/2/pi;
    f=sort(f,'descend');
    fLS=[f(2);f(1)];   %the frequency obtained by LS method
%end