Phase retrieval of a periodic signal, not working well

Hello,

I am trying solve something I thought would be straigfoward, but is not so far ^^.

Let assume we have a perriodic signal as such:

x=0:400;

T=39.5;

PhiO=0.43;

y=5+cos((2*pi/T).*x+PhiO);

And the aim is to get the phase and period of this signal.

Here is for the code for the period retriaval:

%period

[pks,locs] = findpeaks(y,'MinPeakDistance',20) ;

Px_av=mean(diff(locs));

Here is the code for the phase retrieval:

h=fft(y);

plot(abs(h)) % the eleventh element is the frequency of my signal so the eleventh element of the phase "angle(h)" should be the phase of my signal. But if I plot the original signal and the retrieved one, it doesn't match well:

Phi=angle(h);

plot(x,max(y).*cos((2*pi/Px_av).*x+Phi(1,11)));hold on;plot(x,y)

Given that PhiO and Px_av are different, I tried to use interp Phi to match the right frequency, but it doesn't work either.

Do you understand what I do wrong, please ?

Answer 1

[Mathieu NOE]

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))

Px_av = 39.5000

% phase
yc = interp1(x,y,ZxP);
theta = 2*pi*(1-ZxP(1)/Px_av) - pi/2

theta = 0.4203

% 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