# Damped sinusoidal regression using method of integral equations

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.

• Do you sample above the Nyquist limit ?
– user65203
Sep 15, 2017 at 7:26

• Of course it could improve, but not much. The scatter of data will become more and more influential. Definitively, this is not a method convenient in case of large number of periods. It should be more accurate simply to count the number of periods $(N)$ on a long time $(T_N)$ and use the approximate value of the period $T\simeq T_N/N$. With an a-priori known value of $T$ , the other parameters in the damped function can be more easily computed. In case of periodic function on large number of periods, other methods of FFT kind are probably better to evaluate accurately the value of the period. Sep 15, 2017 at 9:02