How to fit damped oscillation curves

Hey folks,

I'm currently trying to get some oscillation characteristics (first/second/third max/min, logarithmic decrement, damping ratio) out of some data I've measured.

My attempts haven't led to useful results yet, most probably due to the non-sinusoidality, noise and/or parasitic sub-peaks of the graphs.

Is there a good way to handle this in matlab? I'm quite a newbie on this.

I would appreciate any help on this topic! Thank you very much in advance!

Answer 1

[Star Strider]

I’m not certain how closely you want to fit it. One option for fitting a sine curve is Curve fitting to a sinusoidal function.

In the context of my Answer there, the ‘fit’ function and fminsearch call for your data become:

fit = @(b,x)  b(1).*exp(b(2).*x).*(sin(2*pi*x./b(3) + 2*pi/b(4))) + b(5);   % Function to fit

s = fminsearch(fcn, [yr;  decay;  per;  -1;  ym])                           % Minimise Least-Squares

with the addition of ‘decay’, the exponential decay constant. An appropriate initial estimate for it is -1. My function is a ‘smoothed’ sinusoidal function, so the fit will be approximate.

I would do a fft of your data to see what other frequency components are in it. You could then create a regression function that would approximate it more closely.

[MrBlub]

Thank you very much for your code. This is what I get when I try to run it on my example:

Error using decay (line 35)
usage: [rate,P] = decay(pdsys,options)

Error in Testfile (line 25)
s = fminsearch(fcn, [yr;  decay;  per;  -1;  ym])                           %
Minimise Least-Squares

This is how I import my data:

SourceFile = 'example.txt';
fid = fopen(SourceFile);
splitlines = textscan(fid, '%f %f');
x = [splitlines{1}];
y = [splitlines{2}];

Any idea?

[Star Strider]

I don’t have the Robust Control Toolbox, so I wasn’t aware I was shadowing a function. Rename the ‘decay’ variable to something else, perhaps simply ‘DK’, or some variable or function that is not already in your code or toolboxes. Or simply replace it in the initial estimates vector with -1. MATLAB doesn’t care, so long as it has an (appropriate) initial estimate for every parameter you want it to estimate.

Since you’re also naming your data the same as in my earlier code, no other changes should be necessary.

I would check the last value estimated in the ‘zx’ (‘zero-crossings’) assignment, to be sure it corresponds to an actual zero-crossing and is not a spurious ‘end-effect’ value. I can’t tell from your plotted data. If it is a spurious value, simply insert:

zx = zx(1:end-1);

just after the initial ‘zx’ assignment. The rest of the code ideally will work without problems. (I don’t have your data, so I can’t experiment with it.)

[MrBlub]

Alright. Fixed it. This is what it tells me now, even if I put in this 'zx' thing:

Exiting: Maximum number of function evaluations has been exceeded
         - increase MaxFunEvals option.
         Current function value: NaN 

    s =

        3.2655
       -1.0000
           NaN
       -1.0000
        0.2446

The red graph is not plotted. Please find attached an example file of my data for your own experiments. Thank you for your time and help!

[Star Strider]

The data I wrote that for were row vectors, so the ‘zx’ function did not calculate the zero-crossings correctly. I didn’t have your data originally, so I couldn’t write the code specifically for it.

This works:

fidi = fopen('MrBlub test.txt','rt');
D = textscan(fidi, '%f%f', 'Delimiter','\t');
x = D{1};
y = D{2};

yu = max(y);
yl = min(y);
yr = (yu-yl);                                                       % Range of ‘y’
yz = y-yu+(yr/2);
zx = x(yz .* circshift(yz,[1 0]) <= 0);                             % Find zero-crossings
per = 2*mean(diff(zx));                                             % 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);   % Function to fit
fcn = @(b) sum((fit(b,x) - y).^2);                                  % Least-Squares cost function
s = fminsearch(fcn, [yr;  -1;  per;  -1;  ym])                      % Minimise Least-Squares

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

figure(1)
plot(x,y,'b',  xp,fit(s,xp), 'r')
grid

The fft reveals two principal frequencies, at 0.08 and 0.27 Hz with amplitudes of 0.36 and 0.07 respectively. You could write a ‘fit’ anonymous function to include both, but I doubt including the second would significantly improve the fit.

The best objective function would be a model of the process that created your data. The parameters would then have direct physical significance.

[MrBlub]

Okay I understand - thank you again for clarifying this. Your code works fine now, but unfortunately it's not really the result I expected. I'm now realizing this code is approximating a sinusoidal function to my data, but what I need is something different. I attached a little sketch: I want Matlab to find the envelope function or at least the values of the first three amplitudes to determine oscillator characteristics like logarithmic decrement, damping ratio...My attempts (manually, origin, matlab) haven't led to useful results yet, I guess there is some kind of additional noise reduction / smoothing of the data necessary.

[Star Strider]

The ‘envelope’ function are the parts of my anonymous function with the first two parameters:

env_fcn = @(b,x) b(1).*exp(b(2).*x);

so the exponential (or logarithmic) decrement is ‘b(2)’.

My ‘fit’ function is the solution of the differential equation of a second-order system. Combinations of the appropriate parameters should provide all the values mentioned in the Wikipedia article on Damping ratio.

I did an fft on your data that revealed one primary peak and one small secondary peak. The fit you got with my code is likely as good as you can get if you are modeling a single second-order system to your data.

[MrBlub]

Hmm ok. I'm not quite sure if we're talking about the same thing. I played around a little bit with Origin and finally got what I need (for this example). Please have a look on the screenshot attached. This is the procedure:

1. Set baseline to y = 0
2. Find global peaks/amplitudes (5 in this example)
3. Find global zero-crossings (for peak widths)
4. Calculate log. decrement = ln (peak1 / peak3)
5. Perform non-linear fitting (exponential decay function) on positive amplitudes (peak 1, 3 and 5) for illustration

So it's not really about fitting the whole graph to one and the same function, which is, I guess, what your code is trying to do, right?

[Star Strider]

It is necessary to fit all the data to the function in order to get the correct parameter estimates. My code does everything you want. In its parameter estimates, s(s2)=-0.084 is the exponential decay, corresponding to a log decrement of 1.088.

This slightly revised code (adding explicitly the log decrement calculation and plotting the envelope) depicts everything:

fidi = fopen('MrBlub test.txt','rt');
D = textscan(fidi, '%f%f', 'Delimiter','\t');
x = D{1};
y = D{2};

yu = max(y);
yl = min(y);
yr = (yu-yl);                                                       % Range of ‘y’
yz = y-yu+(yr/2);
zx = x(yz .* circshift(yz,[1 0]) <= 0);                             % Find zero-crossings
per = 2*mean(diff(zx));                                             % 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);   % Function to fit
fcn = @(b) sum((fit(b,x) - y).^2);                                  % Least-Squares cost function
s = fminsearch(fcn, [yr;  -1;  per;  -1;  ym])                      % Minimise Least-Squares
envlp = @(b,x) b(1) .* exp(b(2) .* x) + b(5);
logdec = sprintf('Log decrement = %.3f', exp(-s(2)));

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

figure(1)
plot(x,y,'b',  xp,fit(s,xp), 'r')
hold on
plot(x, envlp(s,x), 'g', 'LineWidth',1)
hold off
grid
legend('Data', 'Fitted Function', 'Envelope', 'Location','NE')
text(11, 1.1, logdec)

[MrBlub]

You are absolutey right if the damping (e.g. log. decrement) is a constant value. The bad fitting result indicates that this is not the case in this example. With respect to the physical background of this data, I am now expecting the damping to be a time-dependent value. In this case, the best evaluation of this data is to approximate time-dependent damping parameters (e.g. log. decrement) for every existing peak pairs (n, n+1) using the procedure mentioned above. This would enable me to discuss the time-dependent behavior. Hope this additional explanation helps. Thanks again for spending so much time on this.

[Star Strider]

I guess I’m not understanding what you’re doing. I am modeling your data as a linear time-invariant function. That is going to be difficult to model, unless you already know the time-varying features of the system you are measuring. Since I don’t know what the system is, that is likely not something I can help you with.

You can improve the fit somewhat by adding an additional exponentially decreasing sine curve, but that is about the limit of the ability of any parameter estimation technique. There actually seems to be an initial impulse with a rapid exponential decrease all its own to account for the first peak. Since I have no idea what system you’re modeling, I’ll stop guessing at a solution.

Insert these lines in the appropriate parts of the code to see the result of adding the second function:

fit = @(b,x)  b(1).*exp(b(2).*x).*(sin(2*pi*x./b(3) + 2*pi/b(4))) + b(5) + b(6).*exp(b(7).*x).*(sin(2*pi*x./b(8) + 2*pi/b(9)));   % Function to fit

S0 = [yr;  -1;  per;  -1;  ym;  1;  -1;  per*11;  1];
s = fminsearch(fcn, S0)                      % Minimise Least-Squares

envlp = @(b,x) b(1) .* exp(b(2) .* x) + b(5) + b(6).*exp(b(7).*x);
logdec = sprintf('Log decrement = %.3f', exp(-(s(2)+s(7))));

Add an additional exponentially-decreasing ‘spike’ if you want to experiment with it:

fit = @(b,x)  b(1).*exp(b(2).*x).*(sin(2*pi*x./b(3) + 2*pi/b(4))) + b(5) + b(6).*exp(b(7).*x).*(sin(2*pi*x./b(8) + 2*pi/b(9))) + b(10).*exp(b(11).*x);   % Function to fit

S0 = [yr;  -1;  per;  -1;  ym;  1;  -1;  per*11;  10;  15;  -2];

envlp = @(b,x) b(1) .* exp(b(2) .* x) + b(5) + b(6).*exp(b(7).*x) + b(10).*exp(b(11).*x);
logdec = sprintf('Log decrement = %.3f', exp(-(s(2)+s(7)+s(11))));

Nothing else in the code changes. Except for the anomaly of the initial peak, it’s actually a good fit (adding the ‘spike’ and using the last set of parameters).

I really have nothing more to offer after this, since I’m doing nothing more than guessing at the system you’re measuring, and doing my best to model it.

I leave you to experiment.

[MrBlub]

Ok. I do really appreciate your help and think we can come to a good solution if I explain in more detail where this data comes from. I'll give it a try:

The curves we are talking about are measurements of the electrical voltage of an electric arc (air, atmospheric pressure), that is powered by a pulsed kV capacitor discharge. The circuit can be approximated by an RLC-circuit, where R and C are (as already shown) time-variant values R(t) and C(t) describing the electric arc characteristic. In general, we measure exponentially damped oscillatory trains of waves, provided the (initial) arc resistance is not so great as to cause critical damping or over-damping. What we learn from this example: At the first (positive) peak, the voltage reaches (as expected) the highest voltage level depending on the power supplied by the capacitor discharge. From this point of time to the point where the second (negative) peak arises, the arc seems to have the greatest resistance due to low ionization of the surrounding air. So we do have a huge damping (log. decrement 1-2) between peak 1 and peak 2. While time is moving on, the air becomes more and more conductive due to ionization processes taking place as a result of the ongoing energy dissipation of the arc. From an electrical point of view, the arc resistance is therefore more and more decreasing with time and finally becomes (at some point of time) a constant value (full ionization of the air). Therefore the damping (log. decrement) is also decreasing with time and should finally become a constant value, mathematically spoken log.decrement 1-2 > log. decrement 2-3 > log. decrement 3-4 and so on.

So what I’m trying to do is to get a “time-variant log. decrement function” out of this data to mathematically describe this special arc behavior. I hope this is understandable.

[Star Strider]

It is understandable, but it is essentially impossible to model a time-varying process from the output from it, even assuming an impulse input. Without a specific model of the physics of the process (that is beyond my knowledge of physics), there is no way to determine how the parameters vary in time.

I suggest you model the physics of the process, and estimate the parameters from that model. Modeling them as a linear time-invariant system, as I have attempted, won’t work.

I can help you with the MATLAB coding and some of the modeling and parameter estimation, but not the physics. That’s beyond the scope of MATLAB Answers in any event.

Since I can’t help you, I’ll delete this Answer in a few hours, since it’s obviously not relevant to your problem.

[MrBlub]

I see your point and know it will only be an approximate solution as I will not be able to perform such physical modeling and simulation. But I would already be happy if I had a code that numerically determines a, let's say "time-averaged log. decrement", between every single pair of peaks using the formula you can find here

- average damping / log. decrement between peak 1 (pos) and peak 2 (neg)

- average damping / log. decrement between peak 2 (neg) and peak 3 (pos)

- average damping / log. decrement between peak 3 (pos) and peak 4 (neg).

- average damping / log. decrement between peak 4 (neg) and peak 5 (pos)

Don't you think this would be a good first approximation?

[Star Strider]

I don’t understand what you mean by ‘time-averaged log decrement’. This is the best I can do (requires the Signal Processing Toolbox):

fidi = fopen('MrBlub test.txt','rt');
D = textscan(fidi, '%f%f', 'Delimiter','\t');
x = D{1};
y = D{2};

ya = abs(y);

[pks,locs] = findpeaks(ya, 'MinPeakDistance', 200);
logdec = log(pks(1:end-1) ./ pks(2:end));               % Log Decrement
zeta = 1./(sqrt(1 + (2*pi./logdec).^2));                % Damping Ratio

for k1 = 1:length(logdec)
    fprintf(1, '\n\tLog decrement between peaks %2d and %2d = %.4f, Damping ratio = %.4f\n', k1, k1+1, logdec(k1), zeta(k1))
end

figure(3)
plot(x, ya)
hold on
plot(x(locs), ya(locs), 'bp')
plot(x(locs), ya(locs), '-r')
hold off
grid

  Log decrement between peaks  1 and  2 = 1.0986, Damping ratio = 0.1722

  Log decrement between peaks  2 and  3 = 0.0690, Damping ratio = 0.0110

  Log decrement between peaks  3 and  4 = 0.1542, Damping ratio = 0.0245

  Log decrement between peaks  4 and  5 = 0.0870, Damping ratio = 0.0138

The ‘logdec’ assignment is the log-decrement between the absolute values of the adjacent peaks. The ‘zeta’ assignment is the damping ratio for each. Note that the relevant variables are between ‘pks(n)’ and ‘pks(n+1)’, so there are one fewer than the number of peaks. With a nonlinear system where the ratios of the peaks is not close to being constant, taking the time average of data calculated from them would not seem to be appropriate.

[MrBlub]

From a physical point of view, it is in fact a non-linear system due to the complex arc discharge characteristics. If the resistance in this circuit was not an electric arc but a normal constant ohmic resistor, the damping (ratio of the peaks, log. decrement) would be constant and we would have a wonderful damped sinusoid oscillation curve. Unfortunately, the resistance here is an electric arc having a time-dependent resistance. As a consequence, the damping is not constant and we do also have other "instabilities" leading to shown non-normal noise and discontinuities.

Why do you think a "time averaged damping" would not be a good first approximation to describe at least the damping characteristic? If we go back to the example saying that the resistance was a conventional ohmic resistor, we can use the peaks 1 and 2 and calculate the log. decrement using the formula mentioned before . Alternatively we can choose peak 2 and 3 (or 3 and 4, and so on) and we will get the same log. decrement because the damping is a constant value.

Everything changes now when replacing the normal ohmic resistor (R = const) with an ignited electric arc (R = R(t)). We are expecting a huge resistance in the first pulse (due to low ionization / conductivity) corresponding to a huge damping. While time moves on the resistance will more and more decrease, so does the damping. If you want to discuss this behavior using the available data (and this is all I have), how would you do this?

My suggestion was to calculate a "time averaged damping": We know the damping is not constant. So starting from the first peak, we know the damping will decrease in time (in an unknown manner) from each single measured voltage point until the next (neg) peak is reached. So what I suggest is to calculate the log. decrement as usual and say: We dont actually know how in detail the damping changes in the time period from peak 1 to peak 2, but we can say that, as a time-average, the damping is x.

Does this make sense?!

[Star Strider]

You will need to have some function describing ‘R(t)’ to effectively model your system.

I suggested that taking the mean of the damping might not be appropriate because of the time-varying nature of the system.

There are a few ways to configure an RLC circuit. What is your particular model?

[MrBlub]

I know there is some freeware available to simulate RLC circuits. The problem is to model a function describing the electric arc resistance R(t). This results from the time-dependent ionization state of the air which depends on a huge number of complex factors:

- non-linear temperature distribution of the air in the spark gap

- non-linear pressure distribution of the air in the spark gap

- non-linear temperature distribution of the electrodes seperating the spark gap

- non-linear surface condition of the electrodes due to heating and ongoing micro-erosion processes

- stochastic formation of additional discharge filaments (not only one single "arc channel" occuring)

- non-linear arc length(s) due to stochastic motion of the burning arc filament(s)

...

That is why I mentioned before we will not be able to set up an appropriate model to descripte R(t) analytically.

[Star Strider]

Simulating RLC circuits is not the problem, but modeling a specific RLC circuit requires that we know how the components are configured. Is it series, parallel, or some combination? That has a decisive effect on its characteristics and frequency-dependent impedance.

If we can’t at least approximate R(t) with some sort of function, then any attempt at modeling its behaviour is impossible. Spark gaps are not new, and formed the basis of wireless transmitters up to about a century ago. It amazes me that there is no mathematical model of the electrical characteristics of one.

In any event, without knowing the configuration of the RLC circuit model of your system, and having a realistic model of R(t). I can go no further.

Answer 2

[the cyclist]

If you have the Statistics and Machine Learning Toolbox, I would suggest that you look at using the nlinfit function to do a fit like this. You can define your own specific functional form, and fit to that. There are examples on that page that can help you with syntax.

Answer 3

[John D'Errico]

Is there a good way to handle a hard problem where you have no model for the process, where there is nasty, non-normal noise, random bumps, discontinuities? Hmm. not really.

Personally, I'd start by building a shape that does reflect this process. A sine wave, because it is symmetrical, does not do so, even though it is periodic. The ramp up and ramp down behaviors are clearly not the same here. So, were it for me to do (I have seen no actual data, so I cannot even play around with what you have) I would break that curve you have shown into 5 distinct segments. Essentially, I would

1. Locate where the curve crosses zero between each segment.

2. Then shift each piece on top of each other, translating them in time so they now overlap.

3. Now, create a model of this "process". I'd probably play with a principle components solution, that would allow me to model these 5 curves, and find one characteristic curve that represents the overall shape. The PCA analysis would find a basic curve shape, then scaled by a different amount for each segment.

4. I would then build a spline model from that estimated curve. My SLM toolbox would be useful to do this, since I would constrain the curve to pass through zero at each end, use periodicity constraints, etc.

5. Finally, I would now model the entire curve as a piecewise amalgamation of those segments.

I can't suggest much more than that without any actual data.

[Star Strider]

@John — The data are in ‘test.txt’, attached to this comment by MrBlub on 21 Jul 2015 at 17:40.

Answer 4

[Justine Nyakundi]

did you get a solution to this?