function [rad,per] = xsf2phase(S,varargin)
% xs2phase  Convert power spectrum or spectral density matrices to phase shift.
%
% Syntax
% =======
%
%     RAD = xsf2phase(S,...)
%     [RAD,PER] = xsf2phase(S,FREQ,...)
%
% Input arguments
% ================
%
% * `S` [ numeric ] - Power spectrum or spectral density matrices computed
% by the `xsf` function.
%
% Output arguments
% =================
%
% * `RAD` [ numeric ] - Phase shift in radians.
%
% * `PER` [ numeric ] - Phase shift in periods.
%
% Options
% ========
%
% * `'unwrap='` [ `true` | *`false`* ] - Correct phase angles to produce
% smoother phase vector.
%
% Description
% ============
%
% Positive numbers of `RAD` and `PER` mean leads, negative numbers lags.
% Note that this is unlike e.g. the definition by Harvey (1993) where
% positive numbers of phase shifts mean lags.
%
% Example
% ========
%

% -IRIS Toolbox.
% -Copyright (c) 2007-2012 Jaromir Benes.

if ~isempty(varargin) && isnumeric(varargin{1})
    freq = varargin{1};
    varargin(1) = [];
else
    freq = [];
end

opt = passvalopt('freqdom.xsf2phase',varargin{:});

%**************************************************************************

rad = atan2(imag(S),real(S));

if opt.unwrap
    rad = unwrap(rad,[],3);
end

per = [];
if nargout == 1 || isempty(freq)
    return
end

nfreq = length(freq);
per = rad;
realsmall = getrealsmall();
for i = 1 : nfreq
    if abs(freq(i)) <= realsmall
        per(:,:,i,:) = NaN;
    else
        per(:,:,i,:) = per(:,:,i,:) / freq(i);
    end
end

end