function S = xsf(T,R,K,Z,H,D,U,Omega,freq,filter,applyto) %#ok<INUSL>
% xsf  Power spectrum generating function for general state space.
%
% Backend IRIS function.
% No help provided.

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

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

isfilter = ~isempty(filter) && any(applyto);

realsmall = getrealsmall();
ny = size(Z,1);
[nx,nb] = size(T);
nf = nx - nb;
nunit = sum(all(abs(eye(nb) - T(nf+1:end,:)) < realsmall,1));
ne = size(R,2);

Tf = T(1:nf,:);
Rf = R(1:nf,1:ne);
Ra = R(nf+1:end,:);
Ta12 = T(nf+1:nf+nunit,nunit+1:end);
Ta22 = T(nf+nunit+1:end,nunit+1:end);
Sigmaaa = Ra*Omega*transpose(Ra);
Sigmaff = Rf*Omega*transpose(Rf);
Sigmayy = H*Omega*transpose(H);
Sigmafa = Rf*Omega*transpose(Ra);

freq = freq(:)';
nfreq = length(freq);
S = zeros(ny+nf+nb,ny+nf+nb,nfreq);
s = zeros(ny+nf+nb,ny+nf+nb);

status = warning();
warning('off'); %#ok<WNOFF>
for i = 1 : nfreq
    lambda = freq(i);
    if isfilter && filter(i) == 0 && all(applyto) && lambda ~= 0
        continue
    end
    ee = exp(-1i*lambda);
    % F = eye(nf+nb) -  [zeros(nf+nb,nf),T]*exp(-1i*lambda);
    % xxx = F \ Sigmax / ctranspose(F);
    s(ny+1:end,ny+1:end) = doinv();
    s(1:ny,1:ny) = Z*s(ny+nf+1:end,ny+nf+1:end)*transpose(Z) + Sigmayy;
    s(1:ny,ny+1:end) = Z*s(ny+nf+1:end,ny+1:end);
    s(ny+1:end,1:ny) = s(ny+1:end,ny+nf+1:end)*transpose(Z);
    if lambda == 0
        % diffuse y
        index = find(any(abs(Z(:,1:nunit)) > realsmall,2));
        s(index,:) = Inf;
        s(:,index) = Inf;
    end
    if lambda == 0
        % Diffuse xf.
        index = find(any(abs(Tf(:,1:nunit)) > realsmall,2));
        s(ny+index,:) = Inf;
        s(:,ny+index) = Inf;
    end
    if ~isempty(U)
        s(ny+nf+1:end,:) = U*s(ny+nf+1:end,:);
        s(:,ny+nf+1:end) = s(:,ny+nf+1:end)*transpose(U);
        if lambda == 0
            % Diffuse xb.
            index = find(any(abs(U(:,1:nunit)) > realsmall,2));
            s(ny+nf+index,:) = Inf;
            s(:,ny+nf+index) = Inf;
        end
    end
    if isfilter
        s(applyto,:) = filter(i)*s(applyto,:);
        s(:,applyto) = s(:,applyto)*conj(filter(i));
    end
    S(:,:,i) = s;
end
warning(status);

% Skip dividing S by 2*pi.

if ~isfilter
    for i = 1 : size(S,1)
        S(i,i,:) = real(S(i,i,:));
    end
end

%**************************************************************************
    function [Sxx,Saa] = doinv()
        A = Tf*ee;
        % B = inv(eye(nb) - Ta*ee) = inv([A11,A12;0,A22]);
        % A11 = eye(nunit) - eye(unit)*ee
        % A12 = -Ta12*ee
        % A22 = eye(nb-nunit) - Ta22*ee
        B22 = inv(eye(nb-nunit) - Ta22*ee);
        if lambda == 0
            B11 = zeros(nunit);
            B12 = zeros(nunit,nb-nunit);
        else
            d = 1/(1-ee);
            B11 = diag(d(ones(1,nunit)));
            B12 = d*Ta12*B22*ee; %#ok<MINV>
        end
        B = [B11,B12;zeros(nb-nunit,nunit),B22];
        Saa = B*Sigmaaa*ctranspose(B);
        Sfa = Sigmafa*ctranspose(B) + A*Saa;
        X = A*B*transpose(Sigmafa);
        Sff = Sigmaff + X + ctranspose(X) + Tf*Saa*transpose(Tf);
        Sxx = [Sff,Sfa;ctranspose(Sfa),Saa];
    end
% doinv().

end
