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
