How to calculate average for a multiple of time period

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

Answer 1

[Mathieu NOE]

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

[Abhik Saha]

Thank you for the code and the justification. yes you are right to calculate average we need to divide by time.

Thanks!

[Mathieu NOE]

My pleasure !