How do I find the phase for a periodic cos wave using fft?

I have a simple cos wave that has a period, in this case 6. I am offsetting the wave by an angular amount with phase. How can I use the fft function to recover back my initial phase offset? I tried converting the fourier transformed value with the highest amplitude using abs (which does correclty identify the period) back to a phase using angle, but this never seems to work. If someone could enlighten me as to how to do this, I would be so grateful! (Edit: I would like to see how to correctly recover the original period as well)

n = 1e3;
a = linspace(0, 2*pi, n);
period = 6;
phase = rand * 2*pi;
f = rescale(cos((a-phase)*pr));

[Mathieu NOE]

this looks quite similar to this :

Phase retrieval of a periodic signal, not working well - MATLAB Answers - MATLAB Central (mathworks.com)

[Noam]

Thanks for the response! I see in the link you posted that you are using methods other than fft, which is interesting. I could devise other methods as well, but I am trying to learn fft and I would like to use it in this particular scenario, and it seems well within the scope of what fft is capable of solving.

[Mathieu NOE]

sure

NB that in the other post, the author also wanted to use fft, but did not really succeed

Answer 1

[David Goodmanson]

Hi Noam,

here is a version using time and frequency in the usual manner. A couple of ponts here. If you construct a grid as you do with linspace, the cosine will end up having value 0 at the first and last array point. But that does not give a satisfactory result. The array needs to fall one point short at the top end, so that if you think of a set of windows side by side, cos=0 falls into the first point of the next window up. Then the wave is periodic in the desired manner.

In the argument of cosine, the phase is a separate entity and stands apart from multiplication by any factor.

The phase here is set into the equivalent range -pi<= phi<= pi for direct comparison to the angle function.

After using fftshift (which shifts the point corresponding to f=0 to the center of the array; you need to construct a frequency array that matches), positive freq +6 corresponds to array index 507 and negative freq -6 corresponds to array index 495. Since

cos(2*pi*f0*t+phi) = (1/2)*e^(i*(2*pi*f0*t+phi)) + (1/2)*e^(-i*(2*pi*f0*t+phi))

the phase of the positive frequency peak is the same as the orignal phase and the the phase of the negative frequency peak comes in with a minus sign.

n = 1e3;
t = ((0:n-1)/n);
f0 = 6;
phi = -.3
y = cos(2*pi*f0*t + phi);
figure(1)
plot(t,y)
grid on
f = (-n/2:n/2-1);
z = fftshift(fft(y)/n);
figure(2)
stem(f,abs(z))
grid on
xlim([-20 20])
ind2 = max(find(abs(z)>.4))  % find the array index of the positive freq peak
phi2 = angle(z(ind2))
t0 = -phi2/(2*pi*f0)    % since cos peak occurs when (2*pi*f0*t0 + phi) = 0 

[Noam]

Hi David, thank you for this answer! This is much closer to what I am trying to understand. I guess what I am struggling with is how to relate this phi ( phi1, phi2) back to my time series. What I want is to use fft to determine how much my cos wave is offset with respect to the t variable at 0, i.e. what 'time' does my cos wave of frequency f0 'start' at? Or in other words, where is the first peak of the wave after t=0 in terms of t. Thanks again!

[David Goodmanson]

Hi Noam,

I modified the answer to show the time t0 at which the peak occurs. The example has phi = -.3.

The peak in the cos function occurs at the time when its argument = 0. So

2*pi*f0*t0 + phi = 0

and you can then determine t0. The result is t0 = .008 and if you zoom in on the plot in figure 1, that's where the peak is.

[Noam]

Hi David! I finally had some time to get back to this, and your method totally works for my data now!!

Thank you so much! I had to do some translation because my original data is not a time series but actual angular firing. But using your methods I am now successfully recovering the phase!!

Thank you again : )