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.