MATLAB: How to find the best parameters to fit damped oscillations curves

Question

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

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: