# How to calculate average for a multiple of time period

13 views (last 30 days)
Abhik Saha on 24 Aug 2022
Commented: Mathieu NOE on 25 Aug 2022
Hello, Everyone, I have a second order diffrential equation from which I get position as a function of time and velocity as a function of time. (code below) But after that when I try to find the time period(T)=2*pi/omega from the data then I could not find a correct way to calculate that using the data generated. After that I want to calculate average of position and average of velocity. Now the average of position is define in this way
<x>=, same thing for the velocity where n1 is the integer number then I need to plot average value of x for different n1? I am new to matlab please help regarding this. I have copied the code below.
----------------main program----------------------------
t=linspace(0.2,10,1000);
y0=[1 0];
[tsol, ysol]=ode45(@(t,y0) firstodefun2(t,y0), t, y0);
figure(1);clf;
plot(tsol,ysol(:,1),'b')
ylabel('Position')
figure(2);clf;
plot(tsol,ysol(:,2),'r')
ylabel('velocity')
------------------function------------------------------
function dy=firstodefun2(t,y0)
G=1;gamma=1;omega=100;
dy=zeros(2,1);
dy(1)=y0(2);
dy(2)=G*sin(omega*t)-gamma*y0(2);
end

Mathieu NOE on 24 Aug 2022
here your are , my friend
the frequency is estimated by taking the fft of the velocity signal and getting the dominant peak.
but the average displacement definition is not correct : the integral must be divided by the time of integration to get the average (otherwise you get an integral of the displacement but not it's average value)
hope it helps
%----------------main program----------------------------
dt = 1e-3;
t=(0:dt:10);
samples = numel(t);
y0=[1 0];
[tsol, ysol]=ode45(@(t,y0) firstodefun2(t,y0), t, y0);
pos = ysol(:,1);
velo = ysol(:,2);
figure(1);clf;
plot(tsol,pos,'b')
ylabel('Position')
figure(2);clf;
plot(tsol,velo,'r')
ylabel('velocity')
% get omega back from velocity fft
Fs = 1/dt;
[f,P] = one_sidded_fft(velo,Fs);
% find the peak frequency
[PKS,LOCS] = findpeaks(P,'MinPeakHeight',max(P)/2);
freq = f(LOCS);
omega_fft = 2*pi*freq
figure(3),
plot(f,P)
title('velocity fft')
xlabel('Frequency (Hz)');
ylabel('Amplitude')
% average Position (integral over n1 samples from t0)
n1 = 100; % must be less or equal to samples
t1 = t(1:n1);
pos1 = pos(1:n1);
av_Position = trapz(t1,pos1)/t1(end)
%------------------function------------------------------
function dy=firstodefun2(t,y0)
G=1;gamma=1;omega=100;
dy=zeros(2,1);
dy(1)=y0(2);
dy(2)=G*sin(omega*t)-gamma*y0(2);
end
function [f,P] = one_sidded_fft(X,Fs)
samples = length(X);
% regular (non windowed) fft%
P = abs(fft(X)/samples);
% %% hanning window fft%
% P = abs(fft(X.*hanning(samples))*4/samples);
% one sidded fft spectrum % Select first half
if rem(samples,2) % nfft odd
select = (1:(samples+1)/2)';
else
select = (1:samples/2+1)';
end
P = P(select);
P(2:end-1) = 2*P(2:end-1);
f = (select - 1)*Fs/samples;
end
##### 2 CommentsShowHide 1 older comment
Mathieu NOE on 25 Aug 2022
My pleasure !