How can I bulid the exponential envelope of a damped second-order system oscillator?

Hello everyone,

I have a signal made up of points that show the classic second-order system damped oscillations. I would like to build the exponential function that rules the amplitude decay and I did it roughly by interpolating the peaks with an ownmade function that works preatty well. Now the idea is to find out also the coefficients of the exponential decay, but I don't now how to apporach this problem.

I have the complete signal as well as the exponential envelope created by interpolation.

Thank you for the help.

Answer 1

[Star Strider]

One option is presented in Damped harmonic motion curve fit where the exponential envelope function is one of the parameters the code estimates. If you want to do more than that, see Curve Fitting of large Data Measurement?.

[Leonardo Molino]

Thank you very much for the help. May I ask you where can I find some theory behind these very cool algorithms?

[Star Strider]

My pleasure!

I am not certain what you are referring to. The objective function is simply a straightforward expression for a damped sinusoid with a phase term and y-offset. The assignments that set up the initial extimates for the parameters are derived from the data.

The MATLAB functions have their own documentation. I use fminsearch because it does not use the Jacobian to search for a global minimum and so is less likely to be trapped in a local minimum, and because it is part of base MATLAB so no toolboxes are required in order to use it.

[Leonardo Molino]

Your code were so inspiring, thank you!

I would like to show you the resoluts! :)

The data envelope is the one I have obtained with my own function. This is the code. And thanks also for the function fminsearch background :)

D = load('velocity_signal_A.mat');                                           % File Attached
x = D.x;
y = smoothdata(D.y, 'sgolay', 250);    

[ehi, elo] = envelopeFinder(x, y, [], []);
yu = max(y);
yl = min(y);
yr = (yu-yl);                                                                   % Range of ‘y’
yz = y-yu+(yr/2);
zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);                            % Returns Approximate Zero-Crossing Indices Of Argument Vector
zt = x(zci(y));
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) + 2*pi/b(4))) + b(5); 
expenv = @(b,x) [-1 1]*b(1) .* exp(b(2).*x) + b(5);
fcn = @(b) norm(fit(b,x) - y);                                                  % Least-Squares cost function
[s,nmrs] = fminsearch(fcn, [yr; -10;  per;  -1;  ym])                           % Minimise Least-Squares
xp = linspace(min(x),max(x), 500);

figure();
h1 = plot(x,y,'b', 'LineWidth',1.5);
hold on
h2 = plot(xp,fit(s,xp), '--r');
h3= plot(x, expenv(s,x), '-', 'Color',[1 0.5 0], 'LineWidth',2, 'DisplayName', ...
    'Exponential Regression Envelope');
h4 = plot(x, [ehi -elo], '-g', 'LineWidth',1.5, 'DisplayName','Data Envelope');
hold off
grid
xlabel('Time')
ylabel('Amplitude')
legend('Original Data',  'Fitted Curve')
text(0.3*max(xlim),0.7*min(ylim), ...
    sprintf(['$y = ' ...
    ' %.3f\\cdot e^{%.0f\\cdot t}\\cdot sin(2\\pi\\cdot t\\cdot %.0f%.3f)$'],...
    [s(1:2); 1./s(3:4)]), 'Interpreter','latex')
legend([h1; h2; h3(1); h4(1)],'Location','best')

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.

[Leonardo Molino]

Hi!

I have tried to fit this other dynamics but the resoluts are not consistent.

Blue curve are data to fit, while the red one is the fitting.

I have tried with a detrend but the resoluts are the same.

What should I try?

[Star Strider]

My code assumes a constant baseline and a first-order exponential decrease in the envelope. The data in your latest Comment indicate that is not appropriate for your signal, so the model itself will need to change. I am not certain what process produced that result, so the best option would be to use the model of that process for the objective function. The baselline may be a sine curve (or some other nonlinear function such as a polynomial of unknown degree) with a very low frequency, so adding that as a separate additive term would likely be appropriate.

Perhaps using the detrend function first, and observing the results, would help, then fit the detrended result using my original code.

[Leonardo Molino]

Thank you, I will try!

For the records, the equation to be solved is this one here, and rules the dynamics of an aircraft elevator.

The forcing terms are linear functions of the state vector.

I have sucessfully solved this problem with ode45, but as you have correctly said, the analytic expression might be a little bit different from the one of a "regular" pendulum.

[Star Strider]

I appreciate the follow-up.

The model must of course be appropriate for the data.