2
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Do you sample above the Nyquist limit ? $\endgroup$
    – user65203
    Sep 15, 2017 at 7:26

1 Answer 1

2
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ Thanks for the very quick response, Jean! Do you think a more sophisticated numerical integration technique could help here? $\endgroup$
    – Tim W
    Sep 15, 2017 at 8:22
  • $\begingroup$ 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. $\endgroup$
    – JJacquelin
    Sep 15, 2017 at 9:02

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .