How to calculate average for a multiple of time period

Question

Abhik Saha 2022-8-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1785340-how-to-calculate-average-for-a-multiple-of-time-period

评论： Mathieu NOE 2022-8-25

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mathieu NOE 2022-8-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1785340-how-to-calculate-average-for-a-multiple-of-time-period#answer_1032500

在 MATLAB Online 中打开

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