function out = markov_psd_transition_closedform(P, X, f, Fs, remove_mean, nfft)
%MARKOV_PSD_TRANSITION_CLOSEDFORM
% PSD of a stationary Markov transition-output process:
%   y_t = X(s_t, s_{t-1})
% with column-stochastic transition matrix:
%   P(i,j)=Pr{s_t=i | s_{t-1}=j}
%
% Inputs
%   P           NxN, column-stochastic (columns sum to 1)
%   X           NxN, emission on transition j->i is X(i,j)
%   f           frequency grid
%               - if Fs empty: normalized cycles/sample (typically [0,0.5])
%               - if Fs given: Hz (typically [0,Fs/2])
%               If omitted/empty: generated as linspace(0,0.5,nfft) or linspace(0,Fs/2,nfft)
%   Fs          sampling rate in Hz (optional; default [])
%   remove_mean true removes DC impulse (subtract |mu|^2) (default true)
%   nfft        number of frequency samples if f omitted (default 4096)
%
% Output struct out:
%   out.pistat stationary distribution (P*pistat = pistat)
%   out.mu     mean E[y_t]
%   out.R0     R[0] = E[|y_t|^2]
%   out.f      frequency grid (same as input or generated)
%   out.S      PSD on that grid (nonnegative clipped)

    % Defaults
    if nargin < 4, Fs = []; end
    if nargin < 5 || isempty(remove_mean), remove_mean = true; end
    if nargin < 6 || isempty(nfft), nfft = 4096; end

    % If f omitted/empty, make a default one-sided grid
    if nargin < 3 || isempty(f)
        if isempty(Fs)
            f = linspace(0, 0.5, nfft).';
        else
            f = linspace(0, Fs/2, nfft).';
        end
    else
        f = f(:);
    end

    % Ensure numeric full matrices
    P = full(double(P));
    X = full(double(X));

    N = size(P,1);
    if ~isequal(size(P), [N N]) || ~isequal(size(X), [N N])
        error('P and X must both be NxN with the same N.');
    end

    % Normalize columns if slightly off
    cs = sum(P,1);
    if max(abs(cs-1)) > 1e-12
        P = P ./ cs;
    end
    P = max(P,0);

    one = ones(N,1);

    % Stationary distribution pistat: solve (I-P)pistat = 0, sum(pistat)=1
    A = eye(N) - P;
    A(end,:) = one.';
    b = zeros(N,1); b(end) = 1;
    pistat = A\b;
    pistat = real(pistat);
    pistat = max(pistat,0);
    pistat = pistat / sum(pistat);

    % B = X .* P (Hadamard product)
    B = X .* P;

    % Mean mu = E[y]
    mu = one.' * (B * pistat);

    % R0 = E[|y|^2]
    R0 = one.' * ( ((abs(X)).^2) .* P ) * pistat;

    % Optional mean removal (removes DC impulse)
    if remove_mean
        R0_eff = R0 - abs(mu)^2;
    else
        R0_eff = R0;
    end

    % t = B*pistat
    t = B * pistat;

    % Frequency to rad/sample
    if isempty(Fs)
        omega = 2*pi*f;           % normalized cycles/sample
    else
        omega = 2*pi*(f/Fs);      % Hz -> rad/sample
    end

    S = zeros(size(omega));
    I = eye(N);

    % Sum_{k>=1} R[k] e^{-j omega k} = z * 1^T B (I - zP)^{-1} t , z=e^{-jomega}
    for m = 1:numel(omega)
        z = exp(-1j*omega(m));
        u = (I - z*P) \ t;
        tail = z * (one.' * (B * u));   % scalar
        S(m) = real(R0_eff + 2*real(tail));
    end

    % Clip tiny negative numerical noise
    S = max(S, 0);

    out = struct();
    out.pistat = pistat;
    out.mu     = mu;
    out.R0     = R0;
    out.f      = f;
    out.S      = S;
end
