4
$\begingroup$

I'm doing an FFT using Python and Numpy on one machine, and C# on another. I'm using some dummy data that mimics how I'll eventually be gathering data from sensors in the C#/UWP application. The two methods yield consistently identical results on the dummy data, which is great. However, the phase at the sine wave's frequency is always -90°.

Here's the code. Ignore how it essentially does the same thing 7 times, eventually those arrays will be filled with data from 7 different sensors. Sampling frequency is 4096 Hz, signals are gathered for 1 second and are on the interval [-32768,32768]. Note that the x-axes on the Magnitude and Phase graphs are zoomed in to the center frequency plus or minus several Hertz.

import numpy as np
import matplotlib.pyplot as plt
import math

centerf = 1350
span = 2*np.pi*centerf*np.linspace(0,4095,4096)/4096
testData = 32768*np.sin(span)
testData = np.array([math.trunc(x) for x in testData])

plt.close("all")
sensors = ["Channel 1", "Channel 2", "Channel 3", "Channel 4", "Channel 5", "Channel 6", "Channel 7"]

data = []
data.append(testData)
data.append(testData)
data.append(testData)
data.append(testData)
data.append(testData)
data.append(testData)
data.append(testData)

theTime = np.linspace(0,len(data[0]),len(data[0]))/len(data[0])

thisFFT = [np.fft.fft(x) for x in data]

f, a = plt.subplots(7, 1,figsize=(10,8),sharex=True)
f2,a2 = plt.subplots(7,1,figsize=(10,8),sharex=True)
f3,a3 = plt.subplots(7,1,figsize=(10,8),sharex=True)

mag = [np.abs(x)/2048 for x in thisFFT]
phase2 = [np.arctan2(x.imag,x.real)*180/np.pi for x in thisFFT]

mag = [x[:2048] for x in mag] # removed mirrored upper half
phase2 = [x[:2048] for x in phase2]

for x in mag:
    x[0] = 0 # remove DC component

[aa.plot(theTime,d,c='r',lw=0.5) for (aa,d) in zip(a,data)]
[aa.plot(d,c='g',lw=0.75) for (aa,d) in zip(a2,mag)]
[aa.plot(d,c='b',lw=0.75) for (aa,d) in zip(a3,phase2)]
for ch,ax in zip(sensors,a.flat):
    ax.set(ylabel=ch)
for ch,ax in zip(sensors,a2.flat):
    ax.set(ylabel=ch)
for ch,ax in zip(sensors,a3.flat):
    ax.set(ylabel=ch)
a.flat[0].set(title="Amplitude over 1s")
a2.flat[0].set(title="FFT Magnitude")
a3.flat[0].set(title="FFT Phase")
a.flat[5].set(xlabel='Time')

[aa.axvline(x=centerf,lw=0.25) for aa in a2]
[aa.axvline(x=centerf,lw=0.25) for aa in a3]
[aa.set_xlim([centerf-20,centerf+20]) for aa in a2]
[aa.set_xlim([centerf-20,centerf+20]) for aa in a3]

f.tight_layout()
f2.tight_layout()
f3.tight_layout()

print("Magnitude at centerf: {c:1.3f}.".format(c=mag[0][centerf]))
print("Phase at centerf: {c:1.3f}°.".format(c=phase2[0][centerf]))

The output from the print statements is

Magnitude at centerf: 32767.397.
Phase at centerf: -90.000°.

The magnitude is fine, the graph always shows a nice peak right where it should be.

My question: since the input data is a sine wave that always starts at sin(0)=0 and the FFT algorithm uses a rectangular window, shouldn't the phase at the fundamental frequency be 0°? Is there something about phase that I'm not understanding?

$\endgroup$
3
  • 4
    $\begingroup$ The phase is relative to a cosinusoid, since the real part of $e^{j\omega t}$ is $\cos(\omega t)$. $\endgroup$
    – Andy Walls
    Commented Oct 12, 2018 at 19:44
  • $\begingroup$ @AndyWalls so an FFT of a "perfect" sine wave signal will always be 90° out of phase? $\endgroup$
    – Ben S.
    Commented Oct 12, 2018 at 19:51
  • $\begingroup$ FFT of a pure sine wave is two imaginary impulses with opposite signs, yes. Note that the phase of a frequency component really reflects the evenness/oddness of that component in the time domain. $\endgroup$
    – Andy Walls
    Commented Oct 12, 2018 at 20:01

2 Answers 2

3
$\begingroup$

The Fourier Transform (and the DFT is the discrete form of this with a finite length sequence) is simply a correlation of the input signal with each frequency bin, of the continuous form:

$$X(\omega) = \int_{t=0}^Tx(t)e^{-j\omega t}dt$$

(If we consider a finite time duration, however I will stay with the continuous time version of the Fourier Transform as it will simplify the explanation).

Here we see that for any freqency $\omega$ we are doing a correlation of $x(t)$ with $e^{j\omega t}$ (correlation is a complex conjugate multiplication followed by integration, or summation when discrete).

As a correlation we see that the result will be maximum when x(t) is fully correlated to $e^{j\omega t}$ over the given time interval, and similarly 0 when it is fully uncorrelated (orthogonal).

Now make note of Euler's identity given as:

$$e^{-j \omega t} = cos(\omega t) - j sin(\omega t)$$

And replacing that in our first formula we can expand the Fourier Transform given into two correlations:

$$X(\omega) = \int_{t=0}^T x(t)cos(\omega t)dt - j\int_{t=0}^T x(t)sin(\omega t)dt$$

From this we see that when $x(t) =sin(\omega t)$, (and duration T covers one full cycle or an integer number of cycles), the second term will be fully correlated, while the first term will be zero, and vice versa when $x(t) =cos(\omega t)$. Any correlation to the second term results in an imaginary result IF x(t) is real. Any phase shift on x(t), as in $x(t) = sin(\omega t + \phi)$ where $\phi$ is anything except $\pm N\pi/2$ for N any integer, will result in a complex result with $X(\omega)$ having real and imaginary components.

Note the discrete form of this, the DFT, is as given below where the same result can be explained, but may not be as easy to see for some that are less familiar with it.

$$X[k] = \Sigma_{n=0}^{N-1}x[n]e^{-jn\omega_o \frac{k}{N}}$$

$\endgroup$
0
$\begingroup$

This is a common misconception, that a typical FFT implementation shows phase relative to a sin() function that starts at zero the "left" edge of the FFT aperture.

But the sin() function corresponds to the imaginary part of a complex exponential. And the typical default is for a strictly real FFT result to have a phase of zero.

The cos() function also produces a sinewave, but a sinewave such that it is positioned (or offset) to be symmetric around the FFT aperture's center (e.g. an even rather than a odd function). The cos() function corresponds to the real part of a complex exponential.

If you do an atan2(x.imag,x.real) with the real parameter zero and the imaginary parameter non-zero, you get a phase of +-90 degrees.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.