function [hpsd,pwrc,bs,as]=gp_psd(ynorm,rg,fz)

% GP_PSD        Hyperparameter Initialisation using Power Spectrum in
% Frequency Domain.
%
% [HPSD,PWRC,BS,AS] = GP_PSD(ynorm,rg,fz)
%
% This function utilises the power spectrum density (PSD) to analyse
% time-series data in frequency domain, for the use of Gaussian process 
% prior model. The prior covariance function chosen to maximise the 
% log-likelihood function is given of the form,
%
%           C(zi,zj) = a*exp([zi-zj]'*d*[zi-zj]) + bI
%
% where {a, d, b} are hyperparameters of the Gaussian process, such that
% 'as' and 'bs' are approximation values of 'a' and 'b' for the 
% initialisation of the training or optimisation procedure.
%
% Inputs:
%       ynorm   = Normalised data, y, i.e. y-mean(y)
%       rg      = Frequency Range, for calculation using AVGPOWER
%       fz      = Sampling frequency, specified in Hz
%
% Output:
%       HPSD    = Power spectrum density, calculated using PSD
%       PWRC    = Power ratio of noise over data (with noise)
%       BS      = Approximation to hyperparameter 'b'
%       AS      = Approximation to hyperparameter 'a'
%
% The Power Spectrum Density (PSD) is calculated using "Periodogram"
% (Periodogram spectral estimator) approach.
%
%
% (C) Gaussian Process Toolbox 1.0
% (C) Copyright 2006-2007, Keith Neo <keith.neo@nuim.ie>
%                          http://www.hamilton.ie

h=spectrum.periodogram('rectangular');
hopts=psdopts(h);
set(hopts,'NFFT',2^12,'Fs',fz);

hpsd=psd(h,ynorm,hopts);
pwrc=avgpower(hpsd,rg)/avgpower(hpsd);
bs=pwrc/(1-pwrc);
as=(ynorm'*ynorm)/((1+bs)*size(ynorm,1));