function [f, trans, phase, H, info] = ...
    sine_sweep_transw_lsq(THM, a_resp, ijk, nnum, varargin)

% sine_sweep_transw_lsq  ver 3.0
%
% Adds excitation interpretation options:
%
%   'ExcitationType'
%        'force'
%        'measured_base'
%        'solver_base'
%
%   'BaseSign'
%        +1  (default)
%        -1  (flip base input sign)
%
% LSQ estimator itself is unchanged.
%

% ---------------- Options ----------------
p = inputParser;

p.addParameter('CyclesPerWindow', 5, @(x)isnumeric(x)&&isscalar(x)&&x>=2);
p.addParameter('UseHann', true, @islogical);
p.addParameter('MinInAmp', 1e-12, @(x)isnumeric(x)&&isscalar(x)&&x>0);
p.addParameter('Verbose', false, @islogical);

% residual gating
p.addParameter('RejectByResidual', true, @islogical);
p.addParameter('MaxRelResIn',  0.25, @(x)isnumeric(x)&&isscalar(x)&&x>0);
p.addParameter('MaxRelResOut', 0.25, @(x)isnumeric(x)&&isscalar(x)&&x>0);
p.addParameter('MinCoherence', 0.0,  @(x)isnumeric(x)&&isscalar(x)&&x>=0&&x<=1);

% sweep correction
p.addParameter('SweepCorrection', true, @islogical);
p.addParameter('SweepGamma', 1.0, @(x)isnumeric(x)&&isscalar(x)&&x>=0);
p.addParameter('SlopeFitHalfSpan', 6, @(x)isnumeric(x)&&isscalar(x)&&x>=2);

% NEW OPTIONS
p.addParameter('ExcitationType','measured_base',@ischar);
p.addParameter('BaseSign',1,@(x)isnumeric(x)&&isscalar(x));

p.parse(varargin{:});
opt = p.Results;

% ---------------- Basic checks ----------------
t = THM(:,1);
a = THM(:,2);
b = a_resp(:);

t = t(:); a = a(:);

if numel(b) ~= numel(a)
    error('a_resp must have same length as THM(:,2).');
end

dt = mean(diff(t));
if ~isfinite(dt) || dt<=0
    error('Invalid time vector in THM(:,1).');
end

% Apply optional base sign correction
if strcmpi(opt.ExcitationType,'measured_base') || ...
   strcmpi(opt.ExcitationType,'solver_base')
    a = opt.BaseSign * a;
end

% ---------------- Zero crossings ----------------
[tindex, tz] = pos_slope_zero_crossings(t, a);

if numel(tz) < (opt.CyclesPerWindow+2)
    error('Not enough zero crossings to form windows.');
end

% ---------------- Preallocate ----------------
maxK = floor((numel(tz) - opt.CyclesPerWindow - 1) / max(1,nnum)) + 1;
f_all = zeros(maxK,1);
H_all = zeros(maxK,1);

k = 0;
i = 1;

while true

    i2 = i + opt.CyclesPerWindow;
    if i2 > numel(tz), break; end

    ii = tindex(i);
    jj = tindex(i2);

    if jj <= ii+5, break; end

    delta = mean(diff(tz(i:i2)));
    if delta<=0, i=i+max(1,nnum); continue; end

    fk = 1/delta;
    w0 = 2*pi*fk;

    tw = t(ii:jj);
    aw = a(ii:jj);
    bw = b(ii:jj);

    [Xin,~,~,~]  = harmonic_lsq_complex(tw,aw,w0,opt.UseHann);
    [Yout,~,~,~] = harmonic_lsq_complex(tw,bw,w0,opt.UseHann);

    if abs(Xin) < opt.MinInAmp
        i=i+max(1,nnum);
        continue;
    end

    Hk = Yout / Xin;

    % ---------------- Interpretation ----------------
    switch lower(opt.ExcitationType)

        case 'force'
            % FRF already correct

        case 'measured_base'
            % absolute accel transmissibility udd / ydd
            % nothing else required

        case 'solver_base'
            % solver-generated response
            % ratio is still valid
            % but theoretical comparison must use discrete FRF
            % (flag only)

        otherwise
            error('Unknown ExcitationType');
    end

    k = k + 1;
    f_all(k) = fk;
    H_all(k) = Hk;

    i = i + max(1,nnum);

end

f_all = f_all(1:k);
H_all = H_all(1:k);

% -------- Sweep correction (unchanged) --------
Hcorr = H_all;

if opt.SweepCorrection && numel(f_all)>=6
    [k_ln] = local_log_slope( ...
        linspace(0,1,numel(f_all))', f_all, opt.SlopeFitHalfSpan);
    omega = 2*pi*f_all;
    omega_dot = k_ln .* omega;
    corr = 1 + 1i * opt.SweepGamma .* (omega_dot ./ max(eps,omega.^2));
    Hcorr = H_all ./ corr;
end

% -------- Outputs --------
f = f_all(:);
H = Hcorr(:);
trans = abs(H);

if ijk==3
    phase = rad2deg(angle(H));
else
    phase = [];
end

% info
info = struct;
info.ExcitationType = opt.ExcitationType;
info.BaseSign = opt.BaseSign;
info.H_raw = H_all;
info.H = H;

if opt.Verbose
    fprintf('Accepted %d frequency points\n',numel(f));
    fprintf('ExcitationType = %s\n',opt.ExcitationType);
end

end

%% =====================================================================
% Positive slope zero crossing times (with linear interpolation)
%% =====================================================================
function [tindex, tz] = pos_slope_zero_crossings(t, a)

nt = numel(a);
tindex = zeros(nt,1);
tz = zeros(nt,1);

k = 0;
for i = 2:nt
    if (a(i) > a(i-1)) && (a(i) > 0) && (a(i-1) <= 0)
        dy = a(i) - a(i-1);
        dx = t(i) - t(i-1);

        if dx <= 0 || dy == 0
            continue;
        end

        m = dy/dx;
        bq = a(i-1) - m*t(i-1);
        xz = -bq/m;

        k = k + 1;
        tindex(k) = i;
        tz(k) = xz;
    end
end

tindex = tindex(1:k);
tz     = tz(1:k);

end

%% =====================================================================
% Harmonic least-squares projection to get complex amplitude
%
% Model: y(t) = dc + Ac*cos(wt) + As*sin(wt)
% Complex amplitude (cos-reference): Y = Ac - j*As
%% =====================================================================
function [Yc, dc, rms_res, yhat] = harmonic_lsq_complex(t, y, w0, useHann)

t = t(:);
y = y(:);

N = numel(t);
if N < 6
    Yc = 0;
    dc = mean(y);
    rms_res = NaN;
    yhat = dc*ones(size(y));
    return
end

t0 = t(1);
tt = t - t0;

c = cos(w0*tt);
s = sin(w0*tt);

M = [ones(N,1), c, s];

if useHann
    w = hann_local(N);
    Mw = M .* w;
    yw = y .* w;
    x = Mw \ yw;
else
    x = M \ y;
end

dc = x(1);
Ac = x(2);
As = x(3);

Yc = Ac - 1i*As;

yhat = M*x;
r = y - yhat;
rms_res = sqrt(mean(r.^2));

end

function w = hann_local(N)
n = (0:N-1).';
if N == 1
    w = 1;
else
    w = 0.5*(1 - cos(2*pi*n/(N-1)));
end
end

%% =====================================================================
% Local slope fit for LOG sweep: k = d(ln f)/dt using small neighborhood
%% =====================================================================
function k_ln = local_log_slope(tmid, f, halfSpan)

tmid = tmid(:);
f = f(:);

N = numel(f);
k_ln = NaN(N,1);

lf = log(max(f, eps));

for i = 1:N
    i1 = max(1, i-halfSpan);
    i2 = min(N, i+halfSpan);

    tt = tmid(i1:i2);
    yy = lf(i1:i2);

    if numel(tt) < 3
        continue;
    end

    % linear fit yy ≈ a + k*tt
    tt0 = tt - mean(tt);
    yy0 = yy - mean(yy);

    den = sum(tt0.^2);
    if den <= 0
        continue;
    end
    k = sum(tt0.*yy0)/den;     % slope

    k_ln(i) = k;
end

% mild cleanup: replace NaNs at ends with nearest finite
k_ln = fill_ends_with_nearest(k_ln);

end

function x = fill_ends_with_nearest(x)
x = x(:);
idx = find(isfinite(x));
if isempty(idx), return; end
x(1:idx(1)-1) = x(idx(1));
x(idx(end)+1:end) = x(idx(end));
end