incorrect result of fft phase calculation (2024)

44 views (last 30 days)

Show older comments

Lu Gao on 8 Oct 2019

  • Link

    Direct link to this question

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

  • Link

    Direct link to this question

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

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

incorrect result of fft phase calculation (2)

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

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

  • Link

    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

  • Link

    Direct link to this answer

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

  • Link

    Direct link to this answer

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

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

Sign in to comment.

Lu Gao on 14 Oct 2019

  • Link

    Direct link to this answer

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

  • Link

    Direct link to this answer

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

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

incorrect result of fft phase calculation (6)

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

Sign in to comment.

Daniel M on 15 Oct 2019

  • Link

    Direct link to this answer

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

  • Link

    Direct link to this answer

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

  • 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

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

Tags

  • dsp
  • fft
  • angle()

Products

  • MATLAB

Release

R2019a

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.


incorrect result of fft phase calculation (8)

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

incorrect result of fft phase calculation (2024)

FAQs

What is the advantage of the FFT calculation compared to a DFT calculation with respect to complex multiplication of 1024 data points? ›

A 1024 point FFT requires about 70 milliseconds to execute, or 70 microseconds per point. This is more than 300 times faster than the DFT calculated by correlation! Not only is NLog2N less than N2, it increases much more slowly as N becomes larger.

What is the phase factor in FFT? ›

The phase factor is exactly what the name says: If I multiply by such a factor I leave. all magnitudes intact but I impart a certain phase in the complex plane to my value. Thinking in the complex plane all I do is: rotate along the unit circle, not stretch or. contract its radius.

What are the disadvantages of FFT? ›

Its limitations include: (1) information is only provided at discrete frequency steps, so further calculation, for example interpolation, is often used to obtain improved estimates of peak frequencies and amplitudes; (2) 'energy' from spectral peaks may 'leak' into adjacent frequencies, potentially causing lower ...

Is FFT less accurate than DFT? ›

The difference in speed can be enormous, especially for long data sets where n may be in the thousands or millions. In the presence of round-off error, many FFT algorithms are much more accurate than evaluating the DFT definition directly or indirectly.

What is the difference between real FFT and complex FFT? ›

Second, the real Fourier transform only deals with positive frequencies. That is, the frequency domain index, k, only runs from 0 to N/2. In comparison, the complex Fourier transform includes both positive and negative frequencies. This means k runs from 0 to N-1.

What is the formula for Fourier phase? ›

X(jω) = XR(jω) + j XI (jω) where XR(jω), XI (jω) ∈ R. X(jω) = |X(jω)| ej∠X(jω) where |X(jω)|, ZX(jω) ∈ R.

How to calculate the phase of a signal? ›

The phase response of a frequency-domain signal X(ω) is defined by ∠X(ω)=tan−1[Im{X(ω)}Re{X(ω)}], where the Re{X(ω)} represents the real part of X(ω) and Im{X(ω)} represents the imaginary part of X(ω).

How do you find the phase value? ›

In Mathematics, face value is the actual value of the digit in a number. For example, if 567 is a number, then the face value of 6 is 6 only, whereas its place value is tens (i.e. 60). Thus, for any number, having a two-digit, three-digit or 'n' number of digits, every digit will have a place value and a face value.

How do you find the phase? ›

Finding the amplitude, period, and phase shift of a function of the form A × sin(Bx - C) + D or A × cos(Bx - C) + D goes as follows: The amplitude is equal to A ; The period is equal to 2π / B ; and. The phase shift is equal to C / B .

References

Top Articles
Latest Posts
Article information

Author: Kelle Weber

Last Updated:

Views: 6131

Rating: 4.2 / 5 (53 voted)

Reviews: 84% of readers found this page helpful

Author information

Name: Kelle Weber

Birthday: 2000-08-05

Address: 6796 Juan Square, Markfort, MN 58988

Phone: +8215934114615

Job: Hospitality Director

Hobby: tabletop games, Foreign language learning, Leather crafting, Horseback riding, Swimming, Knapping, Handball

Introduction: My name is Kelle Weber, I am a magnificent, enchanting, fair, joyous, light, determined, joyous person who loves writing and wants to share my knowledge and understanding with you.