function [env, freq, mag, xbp] = envelope_fft(x, fs, varargin)
%ENVELOPE_FFT  Envelope (demodulation) + FFT, with optional band-pass pre-filter.
%
%   [ENV, FREQ, MAG] = ENVELOPE_FFT(X, FS) band-pass filters X (optional),
%   extracts the envelope via the Hilbert transform, then returns the
%   single-sided amplitude spectrum MAG of that envelope versus FREQ (Hz).
%
%   [ENV, FREQ, MAG, XBP] = ENVELOPE_FFT(...) also returns the band-pass
%   filtered signal XBP when filtering is on; XBP = 0 when filtering is off.
%
%   [...] = ENVELOPE_FFT(..., 'BANDPASS', TF) turns the band-pass
%   preprocessor ON (TF = true) or OFF (false, default). When ON, X is
%   filtered with a 6th-order Butterworth applied via FILTFILT (forward-
%   backward) for zero phase before the envelope is computed. The filter
%   is designed in second-order-section (SOS) form for numerical
%   robustness.
%
%   [...] = ENVELOPE_FFT(..., 'BAND', [LO HI]) sets the band-pass
%   passband in Hz (only used when 'BANDPASS' is true). Bracket the carrier
%   resonance band. Default: [3000 5000].
%
%   [...] = ENVELOPE_FFT(..., 'WINDOW', W) windows the envelope
%   before the FFT: 'hann' (default) or 'none'.
%
%   [...] = ENVELOPE_FFT(..., 'FMAX', F) trims the returned
%   spectrum to F Hz. Default: full spectrum to Nyquist.
%
%   Inputs:
%       X   - signal (vector); band-pass filtered or raw
%       FS  - sample rate, Hz
%
%   Outputs:
%       ENV  - envelope signal (same length as X)
%       FREQ - frequency vector, Hz (single-sided)
%       MAG  - single-sided amplitude spectrum
%       XBP  - band-pass filtered signal if filtering on, else 0

% ---- parse inputs ----
    p = inputParser;
    p.addRequired('x',  @(v) isnumeric(v) && isvector(v));
    p.addRequired('fs', @(v) isnumeric(v) && isscalar(v) && v > 0);
    p.addParameter('bandpass', false, @(v) islogical(v) || isscalar(v));
    p.addParameter('band', [3000 5000], ...
        @(v) isnumeric(v) && numel(v) == 2 && v(1) > 0 && v(2) > v(1));
    p.addParameter('window', 'hann', ...
        @(s) any(strcmpi(s, {'hann','none'})));
    p.addParameter('fmax', [], @(v) isnumeric(v) && isscalar(v) && v > 0);
    p.parse(x, fs, varargin{:});

    dobp = logical(p.Results.bandpass);
    band = p.Results.band;
    win  = lower(p.Results.window);
    fmax = p.Results.fmax;

    x = double(x(:));
    N = numel(x);

% ---- optional band-pass preprocessor ----
% 6th-order Butterworth, forward-backward (filtfilt) for zero phase.
% Designed in second-order-section (SOS) form to avoid the ill-conditioned
% (b,a) transfer-function representation that triggers near-singular
% RCOND warnings in filtfilt.
if dobp
    if band(1) <= 0
        error('envelope_fft:band', ...
            'Lower band edge (%.1f Hz) must be greater than 0.', band(1));
    end
    if band(2) >= fs/2
        error('envelope_fft:band', ...
            'Upper band edge (%.1f Hz) must be below Nyquist (%.1f Hz).', ...
            band(2), fs/2);
    end
    Wn = band / (fs/2);
    [z,pp,k] = butter(6, Wn, 'bandpass');
    [sos,g]  = zp2sos(z, pp, k);
    x   = filtfilt(sos, g, x);
    xbp = x;        % filtered signal
else
    xbp = 0;        % no filtering
end

% ---- envelope via Hilbert transform ----
    env = abs(hilbert(x));

% ---- remove DC so the 0 Hz line doesn't swamp the defect peaks ----
    e = env - mean(env);

% ---- window ----
switch win
    case 'hann'
        w   = hann(N, 'periodic');
        acg = sum(w) / N;          % coherent gain for amplitude correction
    case 'none'
        w   = ones(N, 1);
        acg = 1;
end
    xw = e .* w;

% ---- single-sided amplitude spectrum ----
    X = fft(xw);
    X = X(1:floor(N/2)+1);
    mag = abs(X) / N / acg;
    mag(2:end-1) = 2 * mag(2:end-1);
    freq = (0:floor(N/2)).' * (fs / N);

% ---- trim to fmax ----
if ~isempty(fmax)
        keep = freq <= fmax;
        freq = freq(keep);
        mag  = mag(keep);
end

end