function [S,D,LIST,FREQ] = xsf(THIS,FREQ,varargin)
% xsf  Power spectrum and spectral density of model variables.
%
% Syntax
% =======
%
%     [S,D,LIST] = xsf(M,FREQ,...)
%     [S,D,LIST,FREQ] = xsf(M,NFREQ,...)
%
% Input arguments
% ================
%
% * `M` [ model ] - Model object.
%
% * `FREQ` [ numeric ] - Vector of frequencies at which the XSFs will be
% evaluated.
%
% * `NFREQ` [ numeric ] - Total number of requested frequencies; the
% frequencies will be evenly spread between 0 and pi.
%
% Output arguments
% =================
%
% * `S` [ namedmat | numeric ] - Power spectrum matrices.
%
% * `D` [ namedmat | numeric ] - Spectral density matrices.
%
% * `LIST` [ cellstr ] - List of variable in order of appearance in rows
% and columns of `S` and `D`.
%
% * `FREQ` [ numeric ] - Vector of frequencies at which the XSFs has been
% evaluated.
%
% Options
% ========
%
% * `'applyTo='` [ cellstr | char | *`Inf`* ] - List of variables to which
% the option `'filter='` will be applied; `Inf` means all variables.
%
% * `'filter='` [ char  | *empty* ] - Linear filter that is applied to
% variables specified by 'applyto'.
%
% * `'nFreq='` [ numeric | *`256`* ] - Number of equally spaced frequencies
% over which the 'filter' is numerically integrated.
%
% * `'output='` [ *`'namedmat'`* | `'numeric'` ] - Output matrices `S` and `F`
% will be either namedmat objects or plain numeric arrays; if the option
% `'select='` is used, `'output='` is always a namedmat object.
%
% * `'progress='` [ `true` | *`false`* ] - Display progress bar on in the
% command window.
%
% * `'select='` [ cellstr | char ] - Select only a subset of variables for
% which the XSF matrices will be returned.
%
% Description
% ============
%
% Example
% ========
%

% -IRIS Toolbox.
% -Copyright 2007-2012 Jaromir Benes.

opt = passvalopt('model.xsf',varargin{:});

if isnumericscalar(FREQ) && FREQ == round(FREQ)
    nfreq = FREQ;
    FREQ = linspace(0,pi,nfreq);
else
    FREQ = FREQ(:).';
    nfreq = length(FREQ);
end

if ischar(opt.select)
    opt.select = regexp(opt.select,'\w+','match');
end

isselect = iscellstr(opt.select);
isnamedmat = strcmpi(opt.output,'namedmat') || isselect;

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

ny = length(THIS.solutionid{1});
nx = length(THIS.solutionid{2});
nalt = size(THIS.Assign,3);

[ans,filter,ans,applyto] = freqdom.applyfilteropt( ...
    THIS,opt,FREQ,ny+nx,THIS.name(THIS.nametype <= 2)); %#ok<NOANS,ASGLU>

if opt.progress
    progress = progressbar('IRIS VAR.xsf progress');
end

S = nan([ny+nx,ny+nx,nfreq,nalt]);
[flag,index] = isnan(THIS,'solution');
for ialt = find(~index)
    [T,R,K,Z,H,D,U,Omega] = mysspace(THIS,ialt,false);
    S(:,:,:,ialt) = freqdom.xsf(T,R,K,Z,H,D,U,Omega,FREQ,filter,applyto);
    if opt.progress
        update(progress,ialt/sum(~index));
    end
end
S = S / (2*pi);

% Solution not available.
if flag
    utils.warning('model', ...
        '#Solution_not_available',sprintf(' #%g',find(index)));
end

% List of variables in rows and columns of `S` and `D`.
LIST = [THIS.solutionvector{1:2}];

% Convert power spectrum to spectral density.
if nargout > 1
    C = acf(THIS);
    D = freqdom.psf2sdf(S,C);
end

% Convert output matrices to namedmat objects.
if isnamedmat
    S = namedmat(S,LIST,LIST);
    if exist('D','var')
        D = namedmat(D,LIST,LIST);
    end
end

% Select variables. For backward compatibility only.
% Use SELECT function afterwards instead.
if isselect
    [S,index] = select(S,opt.select);
    if exist('D','var')
        D = D(index{1},index{2},:,:);
    end
end

end