% helessfit computes the harmonic analysis throught least squares.
% Jorge Armando Laurel-Castillo
% University of Florida
% Febrary 25, 2014
function [meanS,amp,phase,r,s]=halessqfit(omega,time,sample,showFit,col,cc)
%input:
% omega     vector with the desired angular frequencies
% time      vector with the time
% sample    meaasured data 
% showFit   if show the data together the signal fit.
% output:
% amp       amplitud corresponding to the desired angular frequencies
% phase     phase(degrees) corresponding to the desired angular frequencies
m=length(omega);           %number of  harmonics
N=1+2*m;
A=zeros(N,N);
x=zeros(N,1);
b=zeros(N,1);
A(1,1)=length(sample);
b(1)=sum(sample);
for i=1:m                 %idxai=i*2;  idxbi=i*2+1;
    b(i*2)=sum(sample.*sin(omega(i).*time));
    b(i*2+1)=sum(sample.*cos(omega(i).*time));
    
    A(1,i*2)=sum(sin(omega(i).*time)); A(i*2,1)=A(1,i*2);
    A(1,i*2+1)=sum(cos(omega(i).*time)); A(i*2+1,1)=A(1,i*2+1);    
    for j=i:m             %idxaj=j*2;  idxbj=j*2+1;
       %%%% for ai
       A(i*2,j*2)=sum(sin(omega(i).*time).*sin(omega(j).*time)); 
       A(j*2,i*2)=A(i*2,j*2);
       A(i*2,j*2+1)=sum(sin(omega(i).*time).*cos(omega(j).*time));
       A(j*2+1,i*2)=A(i*2,j*2+1);
       %%%% for bi
       A(i*2+1,j*2)=sum(cos(omega(i).*time).*sin(omega(j).*time)); 
       A(j*2,i*2+1)=A(i*2+1,j*2);
       A(i*2+1,j*2+1)=sum(cos(omega(i).*time).*cos(omega(j).*time));
       A(j*2+1,i*2+1)=A(i*2+1,j*2+1);       
    end
end


x=A\b;

meanS=x(1);
aAndb=reshape(x(2:end,1),2,m)';
amp=sqrt( aAndb(:,1).^2 + aAndb(:,2).^2 );
phase=atan2(aAndb(:,2),aAndb(:,1));
tim=time;

if showFit==1
  nData=length(sample);
  ti=time(1);
  tf=time(nData);
  timeF=linspace(ti,tf,nData*3);
  sampleF=0*timeF+meanS;
  sampleP=0*tim+meanS;
  for i=1:m
    sampleF=sampleF+amp(i)*sin(omega(i)*timeF+phase(i));
    sampleP=sampleP+amp(i)*sin(omega(i)*tim+phase(i));
  end
%   plot(time,sample,col,'LineWidth',3);
%   ylabel('')
%   hold on;
%   plot(timeF,sampleF,'color',cc,'LineWidth',3);
%   hold off;
%   grid on
%   grid minor
%   set(gca,'FontSize',16)
end

phase=phase*180/pi;
rnu=sum((mean(sample)-sampleP).^2);
rde=sum((mean(sample)-sample).^2);
r=rnu/rde;

%skill calculation (s)
snu=sum(abs(sampleP-sample).^2);
sde=sum((abs(sampleP-mean(sample))+abs(sample-mean(sample))).^2);
s=1-snu/sde;
end