here a better solution without the need for fft
you should not use the first peak in your code , as it does not correspond to the max amplitude
I prefer to make a "zero crossing" detection a mid y amplitude for both period and phase computation
notice now the period is more accurate
the phase accuracy can be improved if you choose a finer sampling rate
x=0:400;
T=39.5;
PhiO=0.43;
y=5+cos((2*pi/T).*x+PhiO);
%period
[ZxP,ZxN] = find_zc(x,y,mean(y));
Px_av=mean(diff(ZxP))
% phase
yc = interp1(x,y,ZxP);
theta = 2*pi*(1-ZxP(1)/Px_av) - pi/2
% theta = 2*pi*(ZxP(1)/Px_av) + pi/2
plot(x,y,ZxP,yc,'dr');
function [ZxP,ZxN] = find_zc(x,y,threshold)
% put data in rows
x = x(:);
y = y(:);
% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0); %define function: returns indices of +ZCs
ix=zci(y); %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing
ZxP = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
% negative slope "zero" crossing detection, using linear interpolation
zci = @(data) find(diff(sign(data))<0); %define function: returns indices of +ZCs
ix=zci(y); %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing
ZxN = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
end
