Damped sinusoidal regression using method of integral equations

Question

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.

Answer 1

This is not an answer, but a comment which is too long to be posted in the comments section.

I didn't check your code, which is probably correct since it works well for some examples with few periods.

In fact, this is a general feature of the method based on numerical integration. The deviation in numerical calculus mainly comes from the not enough accuracy of the numerical integration, as it is pointed out in the referenced paper, and also from the scatter of experimental data.

In case of periodic or semi-periodic functions on several periods, the accuracy of the computed value of the period (or frequency) is critical. This can be achieved only if they are a sufficient number of points all along the range considered.

When the data covers many periods, the number of points has to be very big, not only proportional to the number of periods but much more, even on each period. This is a drawback of this method of fitting.