function A = seed_amplitudes_ls(f, Hc, fn, zeta, frf, Gdc, tau, Areal, Model)
% SEED_AMPLITUDES_LS  Per-mode amplitudes via local weighted least squares.
% Back-compat signatures:
%   A = seed_amplitudes_ls(f,Hc,fn,zeta)
%   A = seed_amplitudes_ls(f,Hc,fn,zeta,frf)
%   A = seed_amplitudes_ls(f,Hc,fn,zeta,frf,Gdc,tau,Areal,Model)

if nargin < 5 || isempty(frf),   frf   = ''; end
if nargin < 6 || isempty(Gdc),   Gdc   = iff(strcmp(frf,''), 1, 0); end
if nargin < 7 || isempty(tau),   tau   = 0; end
if nargin < 8 || isempty(Areal), Areal = true; end
if nargin < 9 || isempty(Model), Model = 'auto'; end

f=f(:); w=2*pi*f; Hc=Hc(:); if isreal(Hc), Hc=complex(Hc,0); end
fn=fn(:); zeta=zeta(:);
model = iff(strcmpi(Model,'auto'), iff(isempty(frf),'multiplicative','additive'), lower(Model));

switch model
    case 'multiplicative', y = (Hc./max(Gdc,eps)).*exp(1i*w*tau) - 1;
    otherwise,             y = (Hc - Gdc).*exp(1i*w*tau);
end

A = zeros(numel(fn),1);
for k=1:numel(fn)
    wn = 2*pi*fn(k); z = zeta(k);
    col = frf_kernel(w, wn, z, frf);
    bw  = max( 2*z*fn(k), max(1.0, 0.04*fn(k)) );
    Wsq = sqrt(1 + 4*exp(-((f - fn(k)).^2)/(2*(bw/2.355)^2)));
    Xw = col.*Wsq; yw = y.*Wsq;
    if Areal
        num = real(Xw)'*real(yw) + imag(Xw)'*imag(yw);
        den = real(Xw)'*real(Xw) + imag(Xw)'*imag(Xw) + eps;
        A(k) = num/den;
    else
        den = (Xw'*Xw); if abs(den)<eps, A(k)=0; else, A(k)=(Xw'*yw)/den; end
    end
end
end
