How do I get damping constant and periods duration of my damped oscillation

Hey guys,

i m new to matlab and i could need your help. I importet data into matlab which includs the data of a damped oszillation.

i tryed to get my local maxima of the oszillation so i can somehow extract my damping constant. i used the function: islocalmax to find my local maxima.

% Parameterbestimmung für Schwingungsmodell

T = load ("Winkelverlauf_Video2.mat");

t_para = T.T.t/8; % time

theta_para = T.T.theta; % angle of oszillation

figure()

plot(t_para,theta_para);

yline(0.44,'-');

% Lokale Maxima und Minima

t_para_max = t_para(1,1:3000);

theta_para_max = theta_para(1,1:3000);

TF = islocalmax(theta_para_max);

figure()

plot(t_para_max,theta_para_max,'g*')

can u guys help me to find the local maxima of just half a periode and maybe a nice way to find my periode duration, which is about 0,66 s

Answer 1

[Star Strider]

The approach in Calculate time period of a graph signal is one option. That code assumes an exponential envelope, so it will be necessary to change:

exp(b(2).*x)

to:

b(2).*x

to approximate a linearly decreasing envelope.

.

[enrico Rothe]

wow ! Thank you !!!

I m trying to understand what you did there.

x = linspace(0, numel(T1{:,1})-1, numel(T1{:,1}))';

[pks,locs] = findpeaks(T1{:,1});

Period = mean(diff(x(locs)))

i do understand the first line here. but i m not sure what the 2 line means. on the right side of the equation you searching for the peaks in T1 but left side confuses me and how do you get the natural angular frequency of the system in line 3 ?

[enrico Rothe]

if i use this part of your code:

env = envelope(T1{:,1}, 5, 'peak');

figure

plot(x, T1{:,1})

hold on

plot(x,env)

hold off

grid

legend('Signal','Upper Envelope', 'Location','best')

i get that kind of envelope

i think the Problem might be my Signal. Because it is from an image recognition programm it is not very smoth and so it has a lot of peaks. what do you think ?

thats my signal if i soom in:

thanks for your time !!

[Star Strider]

My pleasure!

First, with respect to your first Comment, the first line creates a uniformly-spaced independent variable vector, because that was missing in the data. The second line find the peak values (‘pks’) and their index locations in the original vector (‘locs’). The ‘Period’ assignment returns the estimated period of the oscillation, assuming that they are regularly-spaced. (If they are not, the problem becomes more complicated. That is usually not the situation.)

It is a bit difficult for me to understand your second Comment, since I do not have your data to experiment with. The envelope function call may need to be changed, depending on the sampling frequency of your signal in order to average over the appropriate number of elements and derive the signal envelope rather than tracing the signal. So change (increase) the second argument in the envelope call until it produces the correct signal envelope.

In short, I would have to have at least a rperesentative sample of your signal to experiment with in order to write correct code for it.

[enrico Rothe]

Thanks again ! it worked with np = 100 !

I m now trying to get my damping ration and found some other script.

I uploaded my current script and i think my damping ratio is not quiet correct.

I dont realy understand those 2 lines

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

I do understand, that this equation should show my envelope ( is my Amplitude function) so b(1) is my max_amplitude, right ? So b(2) must be the damping ratio i m looking for ?

B = fminsearch(@(b) norm(theta_para_max(locs) - exp_dk(b,theta_para_max(locs))), rand(2,1).*[1;-1]);

in this line i dont rly understand what <ou are doing. i gues your solving the equation over b(2) ?

Thanks for your valuable help again!

[Star Strider]

With respect to this lilne —

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

I cannot find it in any of the code I posted or referred to. I have no idea where it comes from. The same applies to the fminsearch call. Neither are in the code I posted or referred to, and I have no idea where you discovered them. They are not appropriate to this problem.

In any event, my code can fit the envelope as an exponential, however not as a linear trend.

It returns the exponential constant (-825.4210e-006) and period of the sinusoid (161.6629e-003).

LD = load('Winkelverlauf_Video2.mat');

T = LD.T;

t = T.t(:);

theta = T.theta(:);

theta = fillmissing(theta,'linear');

x = t;

[pks,locs] = findpeaks(theta);

Period = mean(diff(x(locs)))

env = envelope(theta, 250, 'peak');

figure

plot(x, theta)

hold on

plot(x,env)

hold off

grid

legend('Signal','Upper Envelope', 'Location','best')

yc(:,1) = theta;

[mv,ix] = max(yc);

lcs = find(islocalmax(yc(yc(:,1)>0), 'MinProminence',0.05));

lims = [ix max(lcs)];

yc = yc(lims(1):lims(2));

x = linspace(0, numel(yc)-1, numel(yc))';

for k = 1:size(yc,2)

yk = yc(:,k);

y = yk;

% y = detrend(y);

ym = mean(y); % Estimate offset

y = y - ym;

yu = max(y);

yl = min(y);

yr = (yu-yl); % Range of ‘y’

yz = y-yu+(yr/2);

zci = @(v) find(diff(sign(v))); % Returns Approximate Zero-Crossing Indices Of Argument Vector

zt = x(zci(y));

per = 2*mean(diff(zt)) % Estimate period

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

fcn = @(b) norm(fit(b,x) - yk); % Least-Squares cost function

s0 = [yr; -1E-2; per; -1; ym]; % Initial Parameter Estimates

[s,nmrs] = fminsearch(fcn, s0) % Minimise Least-Squares

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

figure

plot(x,yk,':b', 'LineWidth',1.5)

hold on

plot(xp,fit(s,xp), '--r')

hold off

grid

xlabel('Time (s)')

ylabel('Amplitude')

% title(ttlstr{k})

legend('Original Data', 'Fitted Curve')

text(0.1*max(xlim),0.7*min(ylim), sprintf('$y = %.3f\\cdot e^{%.3f\\cdot x}\\cdot sin(2\\pi\\cdot x\\cdot %.3f%+.3f) %+.4f$', [s(1:2); 1./s(3:4); s(5)]), 'Interpreter','latex')

end

.

[enrico Rothe]

https://de.mathworks.com/matlabcentral/answers/754134-plot-and-calculate-damping-ratio?s_tid=srchtitle

here is the question i thought i could ve used for my problem to find delta, the damping ration.

thanks for your code i ll try my best to understand it !

[enrico Rothe]

Hallo Star strider,

i rly appriciat your help but because i m new to matlab and programming the part after : lims = [ix max(lcs)]; (basically to for loop) is hard to understand for me.

my goal is to find the global damping ration.

i guess your for loop....

for k = 1:size(yc,2)

yk = yc(:,k);

y = yk;

% y = detrend(y);

ym = mean(y); % Estimate offset

y = y - ym;

yu = max(y);

is giving me a function which approximates my given data ??

i'm confused wtih the output of the second figure.

it might be because i forgot to transform the 240fps/30fps data to seconds which i did here

is there an easy why to get my global damping ration ?

best regards

[enrico Rothe]

i forgot the updated code...

[Star Strider]

Referring to my code in my previous Comment, put these assignments after the loop to calculate δ, time constant τ,and ξ —

[fitpks,fitlocs] = findpeaks(fit(s,xp));

delta = log(fitpks(end)/fitpks(1))/(numel(fitpks)-1);

tau = 1/s(2);

xi = 1/sqrt(1 + (2*pi/delta)^2);

They evaluate to:

Those are calculated from the exponential model fit.

I also came up with a way to fit the linear portion of the linear descent as opposed to using the exponential function —

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

producing this plot —

Not perfect, although the revised model a much better fit than the exponential model.

These produce different values for the derived constants, however the exponential-model derived constants are more accurate because the fit is exponential, not linear, as is the revised model.

.