function results = fft_autocepstrum_analysis(THM, fig_start, YS, nseg)
% Robust fundamental (harmonic spacing) detector for rotating machinery
%
% Computes BOTH:
%   (1) Averaged RAW cepstrum
%   (2) Averaged ENVELOPE cepstrum (Hilbert)
%
% Linear quefrency axis (seconds)
% Includes low-quefrency liftering, octave correction, and FFT fallback
%
% Inputs:
%   THM       = [time, signal]
%   fig_start = starting figure number
%   YS        = y-axis label
%   nseg      = number of averaging segments
%
% Outputs (struct):
%   results.f0
%   results.T0
%   results.method
%   results.cep_raw
%   results.cep_env
%   results.q
%   results.freq
%   results.fft_mag

fig = fig_start;

%% ---------------- Time History ----------------
t = THM(:,1);
t = t - t(1);
x = THM(:,2);

dt = mean(diff(t));
fs = 1/dt;
n  = length(x);
Trec = t(end) - t(1);

figure(fig)
plot(t, x), grid on
title('Time History')
xlabel('Time (sec)')
ylabel(YS)
fig = fig + 1;

%% ---------------- FFT ----------------
x0 = x - mean(x);
w_full = hann(n,'periodic');
xw = x0 .* w_full;

X = fft(xw);
NHS = floor(n/2);
freq = (0:NHS)'/(n*dt);
FFTmag = abs(X(1:NHS+1));

[~, iMax] = max(FFTmag);
f_peak = freq(iMax);

figure(fig)
plot(freq, FFTmag), grid on
title('FFT Magnitude Spectrum')
xlabel('Frequency (Hz)')
ylabel(YS)
fig = fig + 1;

%% ---------------- Envelope ----------------
env = abs(hilbert(x0));
env = env - mean(env);

%% ---------------- Averaged Cepstra ----------------
nseg = max(1, round(nseg));
segLen = floor(n / nseg);
segLen = max(segLen, 32);

w = hann(segLen,'periodic');

logMag_raw = zeros(segLen,1);
logMag_env = zeros(segLen,1);

for k = 1:nseg
    i1 = (k-1)*segLen + 1;
    i2 = i1 + segLen - 1;
    if i2 > n, break, end

    xr = x0(i1:i2) .* w;
    xe = env(i1:i2) .* w;

    logMag_raw = logMag_raw + log(abs(fft(xr)) + 1e-12);
    logMag_env = logMag_env + log(abs(fft(xe)) + 1e-12);
end

logMag_raw = logMag_raw / nseg;
logMag_env = logMag_env / nseg;

c_raw = real(ifft(logMag_raw));
c_env = real(ifft(logMag_env));

q = (0:segLen-1)' * dt;

%% ---------------- Cepstral Liftering ----------------
q_lifter_cut = 0.010;   % seconds
L = max(2, round(q_lifter_cut/dt));

c_raw(1:L) = 0;
c_env(1:L) = 0;

%% ---------------- Fundamental Search ----------------
f_min = 5;
f_max = min(200, 0.45*fs);

Tmin = 1/f_max;
Tmax = 1/f_min;

roi = (q >= Tmin) & (q <= Tmax);
q_roi = q(roi);

c_raw_s = smoothdata(c_raw(roi),'movmean',max(3,round(0.001/dt)));
c_env_s = smoothdata(c_env(roi),'movmean',max(3,round(0.001/dt)));

noise_raw = median(abs(c_raw_s - median(c_raw_s))) / 0.6745;
noise_env = median(abs(c_env_s - median(c_env_s))) / 0.6745;

minProm_raw = max(3*noise_raw, 0.05*max(abs(c_raw_s)));
minProm_env = max(3*noise_env, 0.05*max(abs(c_env_s)));

[pks_r, locs_r] = findpeaks(c_raw_s, q_roi, ...
    'MinPeakProminence', minProm_raw, 'MinPeakDistance', Tmin);

[pks_e, locs_e] = findpeaks(c_env_s, q_roi, ...
    'MinPeakProminence', minProm_env, 'MinPeakDistance', Tmin);

%% ---------------- Choose Cepstrum ----------------
use_env = ~isempty(pks_e) && (isempty(pks_r) || max(pks_e) > 1.2*max(pks_r));

if use_env
    [~,k] = max(pks_e);
    T0 = locs_e(k);
    method = "cepstrum-envelope";
    q_used = q_roi;
    c_used = c_env_s;
elseif ~isempty(pks_r)
    [~,k] = max(pks_r);
    T0 = locs_r(k);
    method = "cepstrum-raw";
    q_used = q_roi;
    c_used = c_raw_s;
else
    T0 = NaN;
    method = "cepstrum-none";
end

%% ---------------- Octave Correction ----------------
if ~isnan(T0)
    getC = @(T) interp1(q_used, c_used, T, 'linear', 0);
    A1 = getC(T0);
    A2 = getC(T0/2);

    if (T0/2 >= Tmin) && (T0/2 <= Tmax) && (A2 > 0.55*A1)
        T0 = T0/2;
        method = method + "-octave-corrected";
    end
end

f0 = 1/T0;

%% ---------------- FFT Fallback ----------------
if isnan(f0) || f0 < f_min || f0 > f_max
    f0 = f_peak;
    T0 = 1/f0;
    method = "fft-peak";
end

%% ---------------- Adaptive Quefrency Axis ----------------
qmax = min([6*T0, 0.25*Trec, Tmax]);
qmax = max(qmax, 3*T0);

%% ---------------- Cepstrum Plots ----------------
figure(fig)

figH = gcf;
pos = figH.Position;

newW = 1.1 * pos(3);
newH = 1.5 * pos(4);

screen = get(0,'ScreenSize');
newX = screen(1) + (screen(3) - newW)/2;
newY = screen(2) + (screen(4) - newH)/2;

figH.Position = [newX newY newW newH];

subplot(2,1,1)
plot(q, c_raw,'b'), hold on
xline(T0,'r--','LineWidth',1.2)
grid on
xlim([0 qmax])
title('Cepstrum – RAW signal')
xlabel('Quefrency (sec)')
ylabel('Amplitude')
yl = ylim;    % <-- capture y-limits

subplot(2,1,2)
plot(q, c_env,'b'), hold on
xline(T0,'r--','LineWidth',1.2)
grid on
xlim([0 qmax])
ylim(yl)      % <-- enforce same y-limits
title('Cepstrum – ENVELOPE (Hilbert)')
xlabel('Quefrency (sec)')
ylabel('Amplitude')

ax = gca;
fsz = ax.Title.FontSize;

sgtitle(sprintf('T0 = %.5f s   f0 = %.2f Hz   (%s)', ...
    T0, f0, method), ...
    'FontSize', fsz, 'FontWeight','bold');

fig = fig + 1;

%% ---------------- Package Results ----------------
results.freq     = freq;
results.fft_mag  = FFTmag;
results.f_peak   = f_peak;

results.q        = q;
results.cep_raw  = c_raw;
results.cep_env  = c_env;

results.T0       = T0;
results.f0       = f0;
results.method   = method;

fprintf('\nDetected fundamental spacing:\n');
fprintf('  f0 = %.3f Hz   (%s)\n', f0, method);
fprintf('  T0 = %.6f s\n', T0);

end
