How to find the best parameters to fit damped oscillations curves

Hello,

My goal is to fit the best these filtered data (attached) with the equation :

x(1)*exp(-x(2)*t).*sin(2*pi*x(3)*t+x(4)).

Here is my code (below) but this algorithm doesn't give a good fitting.

I don't know if it is a problem of parameters (for example, how to choose x0, ub and lb. I read some information on Matlab website but it doesn't help me)...

Thanks in advance for your help !

fun = @(x,t) x(1)*exp(-x(2)*t).*sin(2*pi*x(3)*t+x(4));
t=[0:1/1000:(1/1000*(length(data_TRT)-1))];
x0 = [2,1,1,1] 
options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt');
lb = [0,0,0,-1]
ub = [2,100,100,1]
[x,resnorm] = lsqcurvefit(fun,x0,t,data_TRT,lb,ub,options)
figure
plot(t,data_TRT,'k-',t,fun(x,t),'r-')
legend('Experimental Data','Modeled Data')
title('Data and Fitted Curve')
xlabel('Temps (secondes)')
Damping_Nigg =x(2)
Frequence_Nigg =x(3)

Answer 1

[Star Strider]

See: How to filter noise from time-frequency data and find natural frequency of a cantilever?

Specifically:

D = load('Louise data_TRT_HELP.mat');
data = D.data_TRT;
t= linspace(0, numel(data), numel(data));
y=data;
Fs=500;
Fn = Fs/2;                                                                      % Nyquist Frequency
L = size(t,2);
y = detrend(y);                                                                 % Remove Linear Trend

fty = fft(y)/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                                             % Frequency Vector
Iv = 1:length(Fv);                                                              % Index Vector

figure(1)
plot(Fv, abs(fty(Iv)))
grid

Wp =  30/Fn;                                                                    % Passband (Normalised)
Ws =  40/Fn;                                                                    % Stopband (Normalised)
Rp =  1;                                                                        % Passband Ripple
Rs = 50;                                                                        % Stopband Ripple
[n,Wp]  = ellipord(Wp,Ws,Rp,Rs);                                                % Filter Order
[z,p,k] = ellip(n,Rp,Rs,Wp);                                                    % Transfer Function Coefficients
[sos,g] = zp2sos(z,p,k);                                                        % Second-Order-Section For Stability

figure(2)
freqz(sos, 2^14, Fs)

yf = filtfilt(sos,g,y);

yu = max(yf);
yl = min(yf);
yr = (yu-yl);                                                                   % Range of ‘y’
yz = yf-yu+(yr/2);
zt = t(yz(:) .* circshift(yz(:),[1 0]) <= 0);                                   % Find zero-crossings
per = 2*mean(diff(zt));                                                         % Estimate period
ym = mean(y);                                                                   % Estimate offset

fit = @(b,x)  b(1) .* exp(b(2).*x) .* (sin(2*pi*x./(b(3)+b(4).*x) + 2*pi/b(5))) + b(6);   % Objective Function to fit
fcn = @(b) sum((fit(b,t) - yf).^2);                                             % Least-Squares cost function
s = fminsearch(fcn, [yr; -1;  per; 0.01;  -1;  ym])                             % Minimise Least-Squares

xp = linspace(min(t),max(t));

figure(3)
plot(t,y,'b')
hold on
plot(t,yf,'g', 'LineWidth',2)
plot(xp,fit(s,xp), 'r', 'LineWidth',1.5)
hold off
grid
title('Sample Data 2')
xlabel('Time')
ylabel('Amplitude')
legend('Original Data', 'Lowpass-Filtered Data', 'Fitted Curve')

producing:

[Louise]

Thanks a lot for your answer ! How do you determine the parameters Wp, Ws, Rp and Rs ?

[Star Strider]

As always, my pleasure!

I decide the passband and stopband frequencies by viewing the Fourier transform of the original data (in figure(1) in this code). I doubt there is any way to define them otherwise. The others are just based on experience, and the amount of stopband attenuation (‘Rs’) I want.

[Louise]

Thanks a lot ! So you have chosen 30 and 40 because it is the interval where there is the largest intensity on the spectra ?

[Star Strider]

As always, my pleasure!

I chose those values to pass the frequencies that appeared to have the most energy, and filter out the high-frequency noise that could interfere with the parameter estimation.

Probably the best way to understand how that works is to experiment with them and see what the results are. Remember that for a lowpass filter, ‘Ws’ must always be greater than ‘Wp’, and with an elliptic filter, by at least 5 Hz (with this sampling frequency) for best results.

[Bjorn Gustavsson]

One thing to consider further is that after the low-pass filtering there might be an "over-filtering" of the initial variations (time < 20). If that's the case then it might be interesting to run the fitting to the original data using the best parameters from the fiting to the filtered data. The general idea is that the filtering reduces the risk of the fiting ending up in some local minima giving you a fit to high-frequency noise, but after the first fitting you would be close enough to the "proper" best-parameters that you now can fit to the original data.

Answer 2

[Bjorn Gustavsson]

If you plot your data you should be able to see that the zero-crossings are increasingly further apart.

That means that your mode-function will not be able to capture the decay nicely. You might have to introduce another parameter in your function to model this chirp.

HTH