I'm trying to perform damped sinusoidal regression using the method outlined in Jean Jacquelin's paper to transform it to a linear regression problem https://fr.scribd.com/doc/14674814/Regressions-et-equations-integrales
This method seems to work great for small numbers of periods - however, as I add additional data (i.e. including additional periods of the sinusoid I am trying to fit to), the solution quickly deteriorates.
This was very counter-intuitive, and I figured I must have made some coding error - my code in MATLAB is below, using a simulated non-damped sinusoid. However, I cannot find an error, and the code does give very good fits for data with a small number of periods.
Has anyone successfully implemented this method? It seems very promising, so any advice you can give would be greatly appreciated!
per = .4;
len = 1; % or poor results with len = 4
x = (0:.001:len)';
y = 2*cos((2*pi)/per*x)+.5*randn(size(x));
n = length(x);
S = zeros(size(x)); % S(1)=0
SS = zeros(size(x)); % SS(1)=0
for k=2:n
S(k) = S(k-1) +(y(k)+y(k-1))*(x(k)-x(k-1))/2;
SS(k) = SS(k-1)+(S(k)+S(k-1))*(x(k)-x(k-1))/2;
end
M = [sum(SS.^2) sum(SS.*S) sum(SS.*x) sum(SS); ...
sum(SS.*S) sum(S.^2) sum(S.*x) sum(S) ; ...
sum(SS.*x) sum(S.*x) sum(x.^2) sum(x) ; ...
sum(SS) sum(S) sum(x) n];
q = [sum(SS.*y);sum(S.*y); sum(x.*y); sum(y)];
v = M\q;
alpha = v(2)/2; omega = sqrt(-v(1)-alpha^2);
beta = sin(omega*x).*exp(alpha*x);
eta = cos(omega*x).*exp(alpha*x);
be = sum(beta.*eta);
M2 = [sum(beta.^2) be;
be sum(eta.^2)];
q2 = [sum(beta.*y); sum(eta.*y)];
v2 = M2\q2;
rho = sqrt(v2(1)^2+v2(2)^2);
phi = atan(v2(2)/v2(1)) + pi*(v2(1)<0);
fit = rho*sin(omega*x+phi).*exp(alpha*x);
figure
hold on
plot(x,y)
plot(x,fit)
edit: thought I might add - this is not due to testing the method on a non-damped sinusoid and causing a near-zero alpha parameter. Swapping the simulated data with
y = (2*cos((2*pi)/per*x)+.5*randn(size(x))).*exp(x*damp);
yields the same issue of increased data leading to a worse fit.