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.
%
%   [...] = 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.
    if dobp
        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);
        [b,a] = butter(6, Wn, 'bandpass');
        x = filtfilt(b, a, 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