How would one handle phase shift fitting an otherwise linear model?

This is a question I was suprised to not find on SE Mathematics or here and it may just be that my theory is weak and this isn't possible.

I have time series data for multiple signals. They are the same piece of physics, so I expect they are the same and would like to describe a function of few parameters that well describes the function. However, the experiment was not ran on highly accurate time and they all are similar waveforms but with different phases. That is, the signals are of the form,

each test case with a different phase, δ, and noise, ϵ, for a few basis function that I've chosen.

My first thought was to use that Fourier transform reduces phase shift to a phase factor, which gives me a linear model - I could just fit the model in fourier space by taking the FT of the basis functions

I went about this by doing the usual least squares form with the complex transpose instead of real transpose

Where y is the observed data, and are the model parameters, , and then I can IFT and the phase would be account for but this doesn't appear to work as there is no phase shift. I'm not sure if I've merely chosen poor basis functions, but even if they aren't great, I would expect some sort of phase shift, for better of for worse.

Where is my issue?

Because it may prove to be a programming mistake, I've included a code sample implementing this

% load signal from experiment, expSig, ftexpSig
dataNum = 7;
fname = sprintf('..\\MatFiles\\data_%02d.mat',dataNum);
load(fname)

% generate domain for model signal
Npts = 2^7;
dt = t(2)-t(1);
tmodel = (0:Npts-1)'*dt-.5;

% create model signal
sqPulse = heaviside(tmodel)-heaviside(tmodel-1);
model = (1-cos(2*pi*tmodel)).*sqPulse;
ftmodel = fft(model);
ftmodel = ftmodel/norm(ftmodel); % not necessary, but eases my conscious

% solve system 
basisFuncs = ftmodel;
coef = (basisFuncs'*basisFuncs)\(basisFuncs'*(ftexpSig));

% apply definition X b
fit = ifft(ftmodel*coef);

figure(1) % spectra of exp and fit
omega = (-Npts/2:Npts/2-1)';
semilogy(omega, fftshift(abs(ftmodel*coef)),...
    omega, fftshift(abs(ftexpSig)))

figure(2) % time series of exp and fit
plot(t, expSig, ...
    tmodel, -real(ifft(ftmodel*coef)),...
    tmodel, sqPulse)

Not sure if I'd rather it be the theory or the code. Thanks.