44 views (last 30 days)

Show older comments

Lu Gao on 8 Oct 2019

Commented: Pepijn Ekelmans on 28 Jul 2021

Open in MATLAB Online

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

##### 1 Comment Show -1 older commentsHide -1 older comments

Show -1 older commentsHide -1 older comments

Pepijn Ekelmans on 28 Jul 2021

#### Direct link to this comment

https://support.mathworks.com/matlabcentral/answers/484265-incorrect-result-of-fft-phase-calculation#comment_1661462

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 :)

Sign in to comment.

Sign in to answer this question.

### Answers (3)

Daniel M on 9 Oct 2019

Edited: Daniel M on 9 Oct 2019

- fftphase.png

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.

##### 0 Comments Show -2 older commentsHide -2 older comments

Show -2 older commentsHide -2 older comments

Sign in to comment.

Lu Gao on 14 Oct 2019

Open in MATLAB Online

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

##### 0 Comments Show -2 older commentsHide -2 older comments

Show -2 older commentsHide -2 older comments

Sign in to comment.

Daniel M on 15 Oct 2019

- phase.png

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.

##### 0 Comments Show -2 older commentsHide -2 older comments

Show -2 older commentsHide -2 older comments

Sign in to comment.

Sign in to answer this question.

### See Also

### Categories

Physical ModelingSimscape ElectricalElectrical Block LibrariesSwitches and Breakers

Find more on **Switches and Breakers** in Help Center and File Exchange

### Community Treasure Hunt

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

Start Hunting!

**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list

Americas

- América Latina (Español)
- Canada (English)
- United States (English)

Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- Deutsch
- English
- Français

- United Kingdom(English)

Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)

Contact your local office