How to determine phase of a sine wave using FFT

Dear MATLAB experts,

For a small project that I want to do, I need to compute the phase of a sine wave. However, I was unable to obtain a correct result. I have posted my code with this message (please see attachment) for your review and comments. If anyone could help me with this, I would greatly appreciate it.

Thank you!

Fuh

[Fuh-Cherng]

Here is the code.

Answer 1

[Star Strider]

The code hasn’t posted. Use the angle function on the output of the fft function:

t  = [...];                                                 % Time Vector
s  = [...];                                                 % Signal Vector

Ts = mean(diff(t));                                         % Sampling Time
Fs = 1/Ts;                                                  % Sampling Frequency
Fn = Fs/2;                                                  % Nyquist Frequency
L  = length(s);

fts = fft(s)/L;                                             % Normalised Fourier Transform
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                         % Frequency Vector
Iv = 1:length(Fv);                                          % Index Vector

amp_fts = abs(fts(Iv))*2;                                   % Spectrum Amplitude
phs_fts = angle(fts(Iv));                                   % Spectrum Phase

Note — This is obviously UNTESTED CODE but it should work.

[Star Strider]

Fuh-Cherng’s ‘Answer’ moved here:

Here is the code. I hope it goes through this time. Sorry about the hassle!

I tested my code again (using the formulas you provided), the resultant phase remained the same, but still incorrect. In the input signal of my example code, the phase is set to 45 degrees. The estimated phase via angle() through fft() however showed a number like 36.2858 degrees. If I plug in a different phase value, the estimated phase is still incorrect. What did I do wrong?

% prepare a sine wave
fs = 5000;                  % sampling rate
dur = 1.2;                  % duration in seconds
n = fs * dur;               % number of data points
t = (0 : n-1) / fs;         % time vector
f = 15;                     % frequency in Hz
a = 1;                      % amlplitude
phi = 45;                   % phase in degrees

x = a * sin(2*pi*f*t + phi * (180/pi));

binwidth = 1 / dur;          % frequency binwidth in Hz
f_index = f / binwidth + 1;  % frequence index

y = abs(fft(x)) / (n/2);            % fft amplitude correction factor equals n/2
phase = angle(fft(x)) * 180/pi;

figure; 
subplot(311); plot(x);
subplot(312); plot(y);
subplot(313); plot(phase);
disp(phase(f_index));       % the estimated phase is incorrect. What did I do wrong???

[Star Strider]

One problem I see are these lines:

phi = 45;                   % phase in degrees

x = a * sin(2*pi*f*t + phi * (180/pi));

You entered ‘phi’ in degrees, then converted it to degrees again. The correct conversion from degrees to radians would be:

x = a * sin(2*pi*f*t + phi * (pi/180));

It’s difficult for me to interpret your plots, but the output of the disp call is now:

-45

Problem solved?

[SMAS]

Hello Star Strider, I can see your answer, I am doing something like Fuh , but I don´t understand the minus and then if I use the same code and only change the phi = 90; the disp is now: -4.8192e-13. I believe the correct answer must be 90, sorry I don´t understand this. Can you help me?

[Min]

Because different trig forms can produce same results, phase output depends on how Matlab interprets the signal. For example, sin(x) and -sin(-x) produce the same output. But also -sin(x) = cos(x + pi/2). That is why the phase output in the above example was -45 (translated by -pi/2) when the value used in the original function was 45 (pi/4). Similarly, your phi = 90 degree (pi/2 in radian) was translated by -pi/2, producing 0. Hope this helps.

[H S]

Dear Star Strider,

I think for the zero frequency (DC), your code after "amp_fts = abs(fts(Iv))*2;" should be corrected as follows:

amp_fts(1)=amp_fts(1)/2;

[Star Strider]

Not according to fft (link).

Answer 2

[Fuh-Cherng]

Thank you so much for pointing out the mistake I made! Greatly appreciated!

[Star Strider]

As always, my pleasure.

I would appreciate it then if you would Accept my Answer.

You can do that by deleting your Answer and then Accepting mine. (No one gets any credit for Accepting their own Answer, anyway.)