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)
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()
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)
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)
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)
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)?
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!)