incorrect result of fft phase calculation

MATLAB users,

I have a strange issue in calculating phase from fft data. On a high level:

• I am using 100kHz sampling frequency to sample 2 sinusoidal signals with the frequency of 100Hz.
• Between the signals, there is a phase shift of pi/4.
• I use single side fft as shown on this page, and angle function to calculate the phase.

The time domain signal and amplitude of fft is correct. However, the phase is incorrect. What I am expecting is a flat line for the first signal (as there is no phase shift) and a pike of pi/4 at 100Hz for the second signal (as the phase shift is pi/4).

I am wondering if I have missed something fundamental here.

Thanks for your comment~!

Lu

clear; close all; clc

% sampling frequency and period
fs = 100E3;
Ts = 1.00 / fs; 

% signal frequncy
fsig = 100; 

% 5 peroids of signal
Tsig = 32 / fsig; 

% time 
T0 = (0.00 : Ts : Tsig - Ts)';

% phase shift
Phi = [0, 1/4]' * pi; 

% time shift calculated out of the phase shift
Tau = Phi / (2.00 * pi * fsig);

YfCollection = [];
YtCollection = [];
PhaseCollection = [];

for it = 1 : length(Tau)
    % yt = sin(2.00 * pi * fsig * (T0 + Tau(it))); 
    yt = sin(2.00 * pi * fsig * T0 + Phi(it)); 
    [f0, Yf, Phase] = calcFFT(yt, fs);
    YtCollection = [YtCollection yt];     
    YfCollection = [YfCollection Yf];
    PhaseCollection = [PhaseCollection Phase];
end

function [f, Y, Phase] = calcFFT(Yin, fs)
    Nfft = 2^(nextpow2(length(Yin))-1); 
    Ycmplx = fft(Yin, Nfft); 
    Y2 = abs(Ycmplx)/Nfft;
    Y = Y2(1 : Nfft/2+1);
    Y(2 : end-1) = 2.00 * Y(2 : end-1);  
    Phase2 = unwrap(angle(Ycmplx));
    f = (0 : Nfft/2)' * fs / Nfft; 
    Phase = Phase2(1 : Nfft/2+1);
end

figure(1)
subplot(3, 1, 1)
for i = 1 : size(YfCollection, 2)
    h(i) = semilogx(f0, YfCollection(:, i), ...
        'DisplayName', sprintf('%s = %4.2f%s', '\phi', Phi(i)/pi, '\pi'));
    hold on
end
xlabel('f (in Hz)')
ylabel('|Y(f)|')
legend(h)
ax = gca;
ax.FontSize = 16;
grid minor

subplot(3, 1, 2)
for i = 1 : size(YfCollection, 2)
    h(i) = semilogx(f0, PhaseCollection(:, i)/pi, ...
        'DisplayName', sprintf('%s = %4.2f%s', '\phi', Phi(i)/pi, '\pi'));
    hold on
end
xlabel('f (in Hz)')
ylabel('phase/\pi')
legend(h)
ax = gca;
ax.FontSize = 16;
grid minor

subplot(3, 1, 3)
for i = 1 : size(YtCollection, 2)
    h(i) = plot(1000*T0, YtCollection(:, i), ...
        'DisplayName', sprintf('%s = %4.2f%s', '\phi', Phi(i)/pi, '\pi'));
    hold on
end
xlabel('time (in ms)')
ylabel('y(t)')
legend(h)
ax = gca;
ax.FontSize = 16;
grid minor

[Pepijn Ekelmans]

Hey Lu,

I fully agree with that this is certainly weird. I personally always learned that FFTs and DFTs were based on the phase of a sine wave, but after some careful experimenting around I have found that the DFTs seems to base its phase on cosines. Maybe a bit farfetched, but this might have to do with phase modulation techniques used in telecommunications. There, an absolute phase of 0 degrees is often associated with a cosine. In the end though, it's all relative, so it's more a human choice than a scientific one :)

Answer 1

[Daniel M]

There are two issues at play. The first is a trade-off between temporal resolution and spectral resolution. Your sampling frequency is so high, that your resolution in the frequency domain is only ~6 Hz. You can do several things for this. The simplest is just to decrease your sampling rate. Instead of 100kHz, try 1 kHz. The other thing you could do is increase the number of samples by increasing the number of periods of oscillation.

By simply changing the sampling rate from 1e5 to 1e3, I was able to produce the attached figure. The value of the phase for the orange line is 0.344pi, and for blue it is 0.098pi, for a difference of ~.25pi, which is the correct phase difference.

But the reason the values are not flat (blue) and pi/4 (orange) is where the second issue comes into play. It is related to round off error of the signal, and you can read all about it the following website, where they address this issue by applying a tolerance threshold before calculating phase.

https://www.gaussianwaves.com/2015/11/interpreting-fft-results-obtaining-magnitude-and-phase-information/

Answer 2

[Lu Gao]

Hi Daniel,

Thanks for your quick response. It solved part of my question, but I still do not get correct answer of the "absolute" phase of my signal. My code is attached below, and particularly:

• To increate the resolution in the frequency domain, I simulated 1000 peroids in the time domain
• In the fft, I set the treshold at the half decay of the maximum magnitude of 1-side FFT
• I also unwrapped the phase
• I also tried cos(omega * t) instead of sin(omega * t)

As you can see in the attached figure, I still see a phase spike at 100Hz for phi = 0 signal (i.e. the blue curve), and a spike for phi = pi/e signal. While the difference between these phases is pi/4, but still, the absolute value of the phase is not correct, and my result is not consistent to the ref. as you provided in the previous message.

What have I missed?

Lu

clear; close all; clc

% sampling frequency and period
fs = 1000;
Ts = 1.00 / fs; 

% signal frequncy
fsig = 100; 

% peroids of signal
Tsig = 1000 / fsig; 

% time 
T0 = (0.00 : Ts : Tsig - Ts)';

% phase shift
Phi = [0, 1/4]' * pi; 

% time shift calculated out of the phase shift
Tau = Phi / (2.00 * pi * fsig);

YfCollection = [];
YtCollection = [];
PhaseCollection = [];

for it = 1 : length(Tau)
    % yt = sin(2.00 * pi * fsig * (T0 + Tau(it))); 
    yt = sin(2.00 * pi * fsig * T0 + Phi(it)); 
    [f0, Yf, Phase] = calcFFT(yt, fs);
    YtCollection = [YtCollection yt];     
    YfCollection = [YfCollection Yf];
    PhaseCollection = [PhaseCollection Phase];
end

figure(1)
subplot(2, 1, 1)
for i = 1 : size(YfCollection, 2)
    h(i) = stem(f0, YfCollection(:, i), 'filled', ...
        'LineWidth', 1.2, ...
        'DisplayName', sprintf('%s = %4.2f%s', '\phi', Phi(i)/pi, '\pi'));
    hold on
end
xlabel('f (in Hz)')
ylabel('|Y(f)|')
legend(h)
grid minor
set(gca, 'XScale', 'log', 'FontSize', 20)

subplot(2, 1, 2)
for i = 1 : size(YfCollection, 2)
    h(i) = stem(f0, PhaseCollection(:, i)/pi, 'filled', ...
        'LineWidth', 1.3, ...
        'DisplayName', sprintf('%s = %4.2f%s', '\phi', Phi(i)/pi, '\pi'));
    hold on
end
xlabel('f (in Hz)')
ylabel('phase/\pi')
legend(h)
set(gca, 'XScale', 'log', 'FontSize', 20)
grid minor

function [f, Y, Phase] = calcFFT(Yin, fs)
    Nfft = 2^(nextpow2(length(Yin))-1); 
    f = (0 : Nfft/2)' * fs / Nfft; 
    Ycmplx = fft(Yin, Nfft); 
    % magnitude
    Y2 = abs(Ycmplx)/Nfft;
    Y = Y2(1 : Nfft/2+1);
    Y(2 : end-1) = 2.00 * Y(2 : end-1);    
    % phase
    threshold = max(Y) / 2;
    Ycmplx1 = Ycmplx(1 : Nfft/2+1);     
    for i = 1 : length(Ycmplx1)
        if(Y(i) <= threshold)
            Phase(i, 1) = 0;
        else
            Phase(i, 1) = angle(Ycmplx(i)); 
        end
    end
    Phase = unwrap(Phase);
end

Answer 3

[Daniel M]

I made only a few changes to your code and got satisfactory results.

Nfft = length(Yin); The zero padding was doing odd things to your signal. It's good for performant processing time (not applicable in this case), but the trade off is precision. Since you are doing a precise experiment here, remove the zero padding.

fs = 5000; You were saturating your signals because you had 100 Hz signal but only 1000 Hz sampling rate. So I made it 5000 and it seems to work fine.

The phase now reads -0.25 for orange, and -0.5 for blue. This makes sense given that sin(x) = cos(x - pi/2). As for the reason the phase is "based upon a cosine wave".... well, I think that goes deep into how the DFT is implemented. Probably has to do with cosine being a symmetric function. But you'll have to do some more research on that.

In fact, when you switch your sin waves with cos, you get 0 and 0.5 for the phase.