# How to find the best parameters to fit damped oscillations curves

18 views (last 30 days)

Show older comments

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)

##### 0 Comments

### Accepted Answer

Star Strider
on 22 Jan 2020

Edited: Star Strider
on 22 Jan 2020

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:

##### 5 Comments

Star Strider
on 23 Jan 2020

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
on 23 Jan 2020

### More Answers (1)

Bjorn Gustavsson
on 22 Jan 2020

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

##### 0 Comments

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!