Phase retrieval from FFT

The following code does not work for certain specific angles like 55 deg. But it works for others and I am not sure why. Please advice if someone has run into this issue. I have tried the angle function as well as atan2d function.

Fs1 = 400e3; % Sampling rate

T1 = 1/Fs % Sampling period

N1 = 40000; % Signal length

t1 = (0:N-1)*T % Time vector

s_1 = cos(2*pi*1500*t1 - deg2rad(20)) - sin(2*pi*4000*t1 - deg2rad(55));

%Getting FFT

y1 = fft(s_1);

z1 = fftshift(y1);

l = length(y1);

%f1 = (-ly1/2:ly1/2-1)/ly1*Fs;

f1 = (-l/2:l/2-1)/l*Fs;

stem(f1,abs(z1))

xlabel 'Frequency (Hz)'

ylabel '|y1|'

grid

%Getting phase

tol1 = 1e-6;

z1(abs(z1) < tol1) = 0;

theta1 = rad2deg(angle(z1));

%theta1 = atan2d(imag(z1),real(z1));

stem(f1,theta1)

xlabel 'Frequency (Hz)'

ylabel 'Phase / \pi'

grid

Answer 1

[Sarvesh Kale]

Hi Shrishti Yadav,

As i understand you are failing to understand the output phase associated with sinusoid of frequency 4000 Hz.

Let us look at the following equation

s_1 = cos(2*pi*1500*t1 - deg2rad(20)) - sin(2*pi*4000*t1 - deg2rad(55));

the above can be rewritten as

s_1 = cos(2*pi*1500*t1 - deg2rad(20)) + cos(2*pi*4000*t1 + deg2rad(90) - deg2rad(55));

we have used the property cos(90 + x) = -sin(x) so it simplifies to the following

s_1 = cos(2*pi*1500*t1 - deg2rad(20)) + cos(2*pi*4000*t1 + deg2rad(35));

the number 35 is a result of simple angle addition of 90 and -55.

So in the phase spectra you will see the phase 35 degrees associated with sinusoid of 4000 Hz.

The uploaded code does not run properly unless few modifications are made, here is a revised code snippet

Fs = 400e3;          % Sampling rate 
T = 1/Fs ;           % Sampling period
N = 40000;             % Signal length
t1 = (0:N-1)*T        % Time vector
s_1 = cos(2*pi*1500*t1 - deg2rad(20)) - sin(2*pi*4000*t1 - deg2rad(50));

%Getting FFT
y1 = fft(s_1);
z1 = fftshift(y1);
l = length(y1);

figure ; 
f1 = (-l/2:l/2-1)/l*Fs;
stem(f1,abs(z1))
xlabel('Frequency (Hz)')
ylabel('|y1|')
grid

%Getting phase
tol1 = 1e-6;
z1(abs(z1) < tol1) = 0;
theta1 = rad2deg(angle(z1));

figure;
stem(f1,theta1);
xlabel('Frequency (Hz)')
ylabel('Phase / \pi')
grid

I hope this answers your queries, please accept the query as answered if satisfied !

[Shrishti Yadav]

Hi Sarvesh,

You didn't change anything in my code but the angle to 50 deg. I understand where the number -35 comes from but that is not the issue. The issue is that this phase recovery should work for all angles the same way. the code should be able to extract individual phase differences for each individual wave in the wave form. The corrected code has the same problem- for a 50 deg it gives me 40 and for an input of a 20 deg it outputs 70. I can correct it by adding 90 - (angle) in the beginning but that is not fixing the issue. Also, why is doing that only for the second term of the waveform and not the first? I have not figured out the glitch yet.

Thanks for the comment.