function [Pr,fr,A,z0,A0,ofac]=vLSP(x, y, f, p0, iofac)
%
%  [Pn, f, A, z0, A0] = vLSP(x, y, f, p0)
%
%  This is a vectorised MATLAB implementation of the Lomb-Scargle's periodogram estimator.
%  Note that this implementation sacrifices space complexity to get a 
%  moderately fast algorithm which avoids loops over the data.
%  
% INPUT
% x - time vector
% y - amplitude vector
% f - frequency vector which define the frequency bands of interest
% p0 - false alarm probability for detection (default=0.05)
% iofac - oversampling factor (default=4)
%
% OUTPUT
% Pr - cell array of periodogram estimates for each frequency band
% fr - cell array of frequencies for each frequency band
% A - vector of integrals of the psd for each frequency band
% z0 - vector of thresholds on psd amplitude for detection of signals in frequency band
% A0 - vector of thresholds on integrals of psd for detection of signals in
% frequency band
% ofac - effective oversampling factor, different from input when the
% lowest frequency required is smaller than the one defined by iofac

% A Eleuteri 23/7/2008

if nargin<5
    iofac=4;
end
if nargin<4
    p0=0.05;
end

x=x(:);
y=y(:);

% center the process
z = y - mean(y);

f=sort(f); % list of frequencies which define the frequency bands of interest
fmax=f(end); % maximum frequency
fmin=f(1); % minimum frequency
T=max(x)-min(x);
ofac=max(iofac,1/(T*fmin)); % if the lowest frequency is not compatible with the default oversampling factor, redefine it
df=1/(T*ofac); % frequency resolution

if ((fmax-fmin)/df*length(x))>200000000
    error('The space complexity of this problem instance is too large. Try using the fLSPw algorithm.');
end


for n=1:(length(f)-1)
    
    fr{n}=[f(n):df:f(n+1)]';

    w=2*pi*fr{n};
    PP=2*x*w';
    TAU=atan2(sum(sin(PP)),sum(cos(PP)))'/2./w;
    tmT=ones(length(TAU),1)*x'-TAU*ones(1,length(x));
    WTMT=(w*ones(1,length(x))).*tmT;
    CW=cos(WTMT);
    SW=sin(WTMT);
    zC=(ones(length(TAU),1)*z').*CW;
    zS=(ones(length(TAU),1)*z').*SW;
    Pn=sum(zC').^2./sum(CW'.^2)+sum(zS').^2./sum(SW'.^2);
    id=(w<=0);
    if any(id)
        Pn(id)=sum(z.*t).^2/sum(t.^2);
    end

    % normalize the psd by the variance
    Pr{n}=Pn/2/var(y);

    % evaluate the total power in the band [f0,f1] by rectangular integration
    A(n)=df*sum(Pr{n});

    % detection statistics
    z0(n)=log(length(fr{n}))-log(p0); % approximate threshold for spectral density peaks. Assumes p0<<1
    A0(n)=gaminv(1-p0,length(fr{n}),df); % threshold for integrated power over band

end
