function IDENT = n4sid_fr(FRF, f, fs, nn)
% N4SID_ID_F performs modal parameter identification using the N4SID (Subspace State-Space System ID) method in the frequency domain.
% INPUTS:
%   FRF - Frequency Response Function data
%   f - Frequency vector corresponding to FRF data
%   fs - Sampling frequency
%   nn - Vector of model orders to test
% OUTPUT:
%   IDENT - Structure containing identified modal parameters (frequency, damping ratio, mode shape) for each model order

% Convert FRF data to an IDFRD object for use with the System Identification Toolbox
data = idfrd(FRF, f, 1/fs);

for i = 1:length(nn)
    % Loop over each specified model order
    % Perform N4SID on the frequency domain data with the current model order
    sys = n4sid(data, nn(i), 'Display', 'off');

    % Extract system matrices A and C for modal analysis
    A = sys.A;                                % System matrix A
    C = sys.C;                                % Output matrix C

    % Compute the eigenvalues and eigenvectors of matrix A
    [V, D] = eig(A);                          % Solve eigenvalue problem for matrix A
    deltat = 1 / fs;                          % Sampling period

    % Extract eigenvalues and compute natural frequencies (w) and damping ratios (z)
    Di = diag(D);                             % Extract eigenvalues
    [w, idx] = sort(abs(log(Di)) / deltat);   % Compute frequencies from eigenvalues and sort
    z = -real(log(Di(idx))) ./ abs(log(Di(idx))); % Compute damping ratios

    % Reorder modal shapes based on sorted frequencies
    PHI = C(:, idx) * V(:, idx);              % Compute and reorder modal shapes
    mod = abs(PHI);                           % Compute magnitude of modal shapes
    mod = mod ./ max(mod, [], 1);             % Normalize modal shapes by their maximum values

    % Compute and adjust phase for each mode shape
    ph = angle(PHI);                          % Compute phase of modal shapes
    ph1 = ph - ph(1, :);                      % Adjust phase by subtracting the phase of the first element

    % Assemble result matrix with natural frequencies (converted to Hz), damping ratios, and normalized modal shapes
    DUMMY = [w'; z'; mod .* sign(cos(ph1))];

    % Store identified parameters, skipping every other column to avoid NaNs
    IDENT(i).ident = DUMMY(:, 2:2:size(DUMMY, 2));

    % Store the current model order in the output structure
    IDENT(i).order = nn(i);
end
end