Search code examples
pythonscipyfft

How to determine the actual phase angles of component waves of a signal processed with scipy.fft?


I used scipy.fft to transform a signal and find its components. The signal was randomly chosen, and the used code is shown below. When manually constructing the component waves, I am confused about the phase angle I need to add to each component.

tt = 10*np.pi
N = 1000
fs = N/tt
space = np.linspace(0,tt,N)
signal = lambda t: 2*np.sin(t) + np.cos(2*t)

The signal looks like this:
original signal

fft1 = fft.fft(signal(space))  
fft1t = fft1/len(space) # scaling the amplitudes
fft1t[1:] = fft1t[1:]*2 # scaling the amplitudes
freq = fft.fftfreq(space.shape[-1], d=1/fs)

The stemplot shows the correctly identified terms with their amplitudes, and the corresponding frequencies are identified from freq as 0.1592 for the 5th term and 0.3183 for the 10th term.

plt.stem(np.abs(fft1t)[:int(len(fft1t)/2)]) # plotting only up to the Nyquist freq.
plt.xlim(-.5,20)
plt.grid()

stemplot

Now when I plot try to manually reconstruct the signal, and plot is on top of the original, I get a phase mismatch by 2 times the phase angle calculated from the fft results.

uk = space**0 * fft1t[0].real
for n in [5,10]: # 5th and 10th term of the Fourier series
    uk += (fft1t[n].real * np.cos(n * np.pi * space / (N/fs/2) + 0*np.angle(fft1t[n])) 
         + fft1t[n].imag * np.sin(n * np.pi * space / (N/fs/2) + 0*np.angle(fft1t[n])))
sns.lineplot(x=space, y=signal(space))
sns.lineplot(x=space, y=uk) 

phase angle 0

When I add one phase angle I get:

uk = space**0 * fft1t[0].real
for n in [5,10]: # 5th and 10th term of the Fourier series
        uk += (fft1t[n].real * np.cos(n * np.pi * space / (N/fs/2) + 1*np.angle(fft1t[n])) 
             + fft1t[n].imag * np.sin(n * np.pi * space / (N/fs/2) + 1*np.angle(fft1t[n])))
sns.lineplot(x=space, y=signal(space))
sns.lineplot(x=space, y=uk)

phase angle 1

And when I add two times the phase angle, then I get a perfect match:

uk = space**0 * fft1t[0].real
for n in [5,10]: # 5th and 10th term of the Fourier series
        uk += (fft1t[n].real * np.cos(n * np.pi * space / (N/fs/2) + 2*np.angle(fft1t[n])) 
             + fft1t[n].imag * np.sin(n * np.pi * space / (N/fs/2) + 2*np.angle(fft1t[n])))
sns.lineplot(x=space, y=signal(space))
sns.lineplot(x=space, y=uk)

phase angle 2

Please help me understand why am I required to double the calculated phase angle here, and what is the general rule for determining the phase angle to add. In another question (example) I recently posted, I only had to add one times the phase angle. Is this connected to the sampling space, where instead of sampling from [-T,T] (like in the example from the link) we sample in [0,T] (like in this example)?


Solution

  • The term you compute for each frequency component is wrong.

    If you want to use the phase, then you should also use the magnitude:

    p = n * np.pi * space / (N/fs/2)
    np.abs(fft1t[n]) * np.cos(p + np.angle(fft1t[n]))
    

    That is, each frequency component represents a single shifted sine wave. A phase of 0 corresponds to an even function, hence we use the cosine, not the sine, to build that function.

    If you want to use complex multiplication, then you don’t need to compute the phase separately, it is encoded in the complex numbers. The equation now is

    fft1t[n] * (np.cos(p) + 1j * np.sin(p))
    

    (which is equivalent to the multiplication with exp(…) in the other answer). But this form will produce both a real and imaginary component. The imaginary part would be cancelled when considering the corresponding negative frequency component. But you already discarded that component, and doubled the positive one, meaning the positive and negative components are already combined in some way. And thus we need to ignore the imaginary part that comes out:

    fft1t[n] * (np.cos(p) + 1j * np.sin(p)).real
    

    If we want to be a bit more efficient and never compute that imaginary component, then we need to properly expand this complex multiplication:

    fft1t[n].real * np.cos(p) - fft1t[n].imag * np.sin(p)
    

    (note the minus!)