function [S,D] = xsf(this,freq,varargin)
% xsf  Power spectrum and spectral density functions for VAR variables.
%
% Syntax
% =======
%
%     [S,D] = xsf(V,FREQ,...)
%
% Input arguments
% ================
%
% * `V` [ VAR ] - VAR object.
%
% * `FREQ` [ numeric ] - Vector of Frequencies at which the XSFs will be
% evaluated.
%
% Output arguments
% =================
%
% * `S` [ numeric ] - Power spectrum matrices.
%
% * `D` [ numeric ] - Spectral density matrices.
%
% Options
% ========
%
% * `'applyTo='` [ cellstr | char | *`':'`* ] - List of variables to which
% the `'filter='` will be applied; `':'` 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.
%
% * `'progress='` [ `true` | *`false`* ] - Display progress bar in the command
% window.
%
% Description
% ============
%
% The output matrices, `S` and `D`, are `N`-by-`N`-by-`K`, where `N` is the
% number of VAR variables and `K` is the number of frequencies (i.e. the
% length of the vector `freq`).
%
% The k-th page is the `S` matrix, i.e. `S(:,:,k)`, is the cross-spectrum
% matrix for the VAR variables at the k-th frequency. Similarly, the `k`-th
% page in `D`, i.e. `D(:,:,k)`, is the cross-density matrix.
%
% Example
% ========
%

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

% Parse input arguments.
IP = inputParser();
IP.addRequired('VAR',@(x) isa(x,'VAR'));
IP.addRequired('freq',@isnumeric);
IP.parse(this,freq);

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

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

ny = size(this.A,1);
nalt = size(this.A,3);
freq = freq(:)';
nfreq = length(freq);

[ans,filter,ans,applyto] = ...
    freqdom.applyfilteropt(this,opt,freq,ny,this.ynames); %#ok<NOANS,ASGLU>

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

S = nan(ny,ny,nfreq,nalt);
for ialt = 1 : nalt
    % Compute power spectrum function.
    S(:,:,:,ialt) = freqdom.xsfvar( ...
        this.A(:,:,ialt),this.Omega(:,:,ialt),freq,filter,applyto);
    if opt.progress
        update(progress,ialt/nalt);
    end
end
S = S / (2*pi);

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

end
