function [a_pos, a_neg] = arbit_engine_alt(a1,a2,b1,b2,b3,yy,f,res)
% arbit_engine_alt.m
%
% Stable batch-recursive Smallwood filter implementation.
% Vectorized over nspec filters, looping only over time.
%
% INPUTS:
%   a1,a2,b1,b2,b3 : nspec×1 Smallwood coefficients
%   yy             : Nt×1 input time history
%   f              : nspec×1 natural frequencies (only f(1) used if res==1)
%   res            : 1 = apply zero padding (ring-down), 0 = no padding
%
% OUTPUTS:
%   a_pos          : nspec×1 positive response peaks
%   a_neg          : nspec×1 negative response peaks (magnitude)

    % ---- Shape & type enforcement ----
    yy = double(yy(:));

    a1 = double(a1(:));
    a2 = double(a2(:));
    b1 = double(b1(:));
    b2 = double(b2(:));
    b3 = double(b3(:));

    nspec = length(a1);

    % Basic length consistency check
    if ~(length(b1)==nspec && length(b2)==nspec && length(b3)==nspec && length(a2)==nspec)
        error('Coefficient length mismatch: a1=%d, a2=%d, b1=%d, b2=%d, b3=%d.', ...
            length(a1), length(a2), length(b1), length(b2), length(b3));
    end

    % ---- Optional zero padding (ring-down) ----
    if res == 1
        Nres = round(0.75 / f(1));   % same logic as ws_srs
        yy   = [yy ; zeros(Nres,1)];
    end

    Nt = length(yy);

    % ---- Initialize states & response storage ----
    y1 = zeros(nspec,1);     % y[n-1] for all filters
    y2 = zeros(nspec,1);     % y[n-2] for all filters

    a_resp = zeros(Nt, nspec);   % time × nspec

    % ---- Main batch-recursive loop ----
    for n = 1:Nt

        a0 = yy(n);

        if n >= 2
            a1s = yy(n-1);
        else
            a1s = 0;
        end

        if n >= 3
            a2s = yy(n-2);
        else
            a2s = 0;
        end

        % Vectorized feed-forward term for all nspec filters
        % u = b1*x[n] + b2*x[n-1] + b3*x[n-2]
        u = b1 * a0 + b2 * a1s + b3 * a2s;   % nspec×1

        % Vectorized recursive term:
        % y[n] = u + a1.*y[n-1] + a2.*y[n-2]
        % (signs are assumed already baked into a1,a2 from Smallwood)
        y = u + a1 .* y1 + a2 .* y2;        % nspec×1

        a_resp(n,:) = y.';   % store row

        % Shift states
        y2 = y1;
        y1 = y;
    end

    % ---- Peak extraction ----
    ymax = max(a_resp, [], 1).';   % nspec×1
    ymin = min(a_resp, [], 1).';   % nspec×1

    a_pos = abs(ymax);
    a_neg = abs(ymin);
end
