function [dg,dag,dah,g_Riemannian]=dg(f1,f2)

% function [dg,dag,dah,g_Riemannian]=dg(f1,f2) : dg geodesic spectral distance between f1 and f2
%
% dg : geodesic distance between power spectral densities f1 and f2
%
% dg is a (pseudo-)metric between spectral density functions which
%    is intrinsic and separates deterministic from nondeterministic spectra
%    (at least those which are log-L2 integrable.
%
% Then dag, dah are the following "divergences"
%
% dag : logarithm of the ratio of arithmetic over geometric means of f1/f2
% dah : logarithm of the ratio of arithmetic over harmonic means of f1/f2
%
% Then,
% g_Riemannian(f1,Delta) where Delta=f2=perturbation from f1,
%                       : a Riemannian metric quantifying degradation
%                         of predictive error variance when f1 is replaced
%                         by f1+Delta. It is the metric which gives rise to the geodesic
%                         distance dg.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N=length(f1); f=f1./f2; logf=log(f);

dg=sqrt(sum(logf.^2)/N-(sum(logf)/N)^2);
dag=log(sum(f)/N)-sum(logf)/N;
dah=log(sum(f)/N)+log(sum(1./f)/N);
dg=sum((1./f).^2)/N-(sum(1./f)/N)^2;