How to extract the correct phase of a sine wave from data

101 views (last 30 days)
Hi everyone
I'm trying to write a code which allows me to determine the phase of a sine wave from its vector.
given: the sinusoidal signal vector, its frequency
unknown: its phase
I have tried with this code but i don't understand why do i end up getting wrong/ unprecise results
a=0; b=800; %time intervalls
N=20000; %N number of samples
t=linspace(a,b,N); %time vector
ts=(b-a)/N; %sampling time
fs=1/ts; %sampling frequency
fn=fs/2; %Nyquist frequency
fx=1; phix= 180; %Frequency and phase of x(t)
x=sin(2*pi*fx.*t+phix* (pi/180)); %Original signal
binwidth = 1 / (b-a); % frequency binwidth in Hz
f_index = fx / binwidth +1; % frequency index
y = abs(fft(x)) / (N/2);
phase = angle(fft(x))* 180/pi+90;
disp(phase(f_index));
for example here i get 187.2000 instead of 180
I'd be very grateful for your help and clarifications!

Accepted Answer

Freewing
Freewing on 24 Oct 2020
The reason you get imprecise results is because you are using discrete fourier transform. When you get the phase of the frequency bin from fft, you get, roughly speaking, the phase of a frequency bin, but not the phase of a single frequency. You can get much better estimation if you analyse a few neighboring bins, but that’s a complex process.
If you’d like to compute exact phase and you know the signal is exactly sinusoidal, you’d better just use asin(x(1)).
  1 Comment
bob98
bob98 on 25 Oct 2020
thank you for your clarification,but i can't use asin(x(1)) since i don't know the amplitude of the signal.

Sign in to comment.

More Answers (2)

Ameer Hamza
Ameer Hamza on 24 Oct 2020
That is because you are doing a discrete Fourier transform. There can be a bit of spectral leakage to other frequency components if the sampling rate is low. For example, run your code after setting.
N=20000000;
and check the result.

Brenton Barritt
Brenton Barritt on 18 Jan 2022
Your code is fine except for: t=linspace(a,b,N);
It seems right to finish the t vector on exactly 800, but in fact 800 is the first point of the next sequence. As an example, to produce one decade of integers starting from zero, linspace(0, 10, 10) seems right but gives ugly numbers! linspace(0, 9, 10) is actually what we need. Including '10' is wrong, that's the start of the next decade.
So what you really want to generate is: t = linspace(a, b-b/N, N);
Now your code gives 180.0000, with no need to blame spectral leakage!

Products


Release

R2019b

Community Treasure Hunt

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

Start Hunting!