# Phase retrieval of a periodic signal, not working well

3 views (last 30 days)
Tual monfort on 14 Jun 2024
Answered: Mathieu NOE on 14 Jun 2024
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 ?

Mathieu NOE on 14 Jun 2024
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

### Categories

Find more on Parametric Spectral Estimation in Help Center and File Exchange

R2024a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!