1
votes

I'm looking for a clarification of Fourier transform principles.

I try to do something quite simple: Create a signal (sine wave with a given frequency and phase shift) and recreate its params with Fourier transform. Frequency estimations work fine, but when it comes to phase, it looks like I get systematic shift (-pi/2).

import numpy as np
import matplotlib.pyplot as plt

duration = 4.0 # lenght of window (in sec)
ticks_per_sec = 400.0 # sampling interval 
samples = int(ticks_per_sec*duration)

phase_shift = np.pi / 2 # sin wave shift in angle
freq = 1 # sine wave freq in Hz

phase_shift  = round(ticks_per_sec/freq * phase_shift/(2*np.pi)) # angle value translated to no of ticks
t = np.arange(phase_shift, samples+phase_shift) / ticks_per_sec

s =  1 * np.sin(2.0 * np.pi * freq * t)
N = s.size

fig, axs = plt.subplots(1, 3, figsize=(18, 6))
axs[0].grid(True)
axs[0].set_ylabel("Amplitude")
axs[0].set_xlabel("Time [s]")
axs[0].set_title(f"F: {freq}, Phase shift: {phase_shift} ticks.")
axs[0].plot(np.arange(samples)/ticks_per_sec, s) 

f = np.linspace(0, ticks_per_sec, N)
fft = np.fft.fft(s)
peak_pos = np.argmax(np.abs(fft[:N//2])) 

axs[1].set_ylabel("Amplitude")
axs[1].set_xlabel("Frequency [Hz]")
axs[1].set_title(f"Peak bar: {peak_pos}")
barlist = axs[1].bar(f[:N // 2], np.abs(fft)[:N // 2] * (1 / (N//2)), width=1.5)  # 1 / N is a normalization factor
barlist[peak_pos].set_color('r')


axs[2].set_ylabel("Angle")
axs[2].set_xlabel("Frequency [Hz]")
axs[2].set_title(f"Peak angle: {np.angle(fft[peak_pos])}")
barlist = axs[2].bar(f[:N // 2], np.angle(fft)[:N // 2], width=1.5) 
barlist[peak_pos].set_color('r')

fig.show()

Plotted Results of the code above

Please help me if there's a bug in my code that I can't notice, or I misunderstand something.

Thank you in advance.

1

1 Answers

1
votes

Your code is just fine, this is not a programming issue.

Let's recall that a sine wave can be expressed as a cosine wave with a phase shift (or vice versa), now remember that sine function as an inherent phase shift of -pi/2 in real Fourier basis relatively to cosine.

This means that your code should output a pi/2 phase angle when replacing np.sin by np.cos, i.e. returns input phase_shift, or equivalently, returns a phase angle of zero when specifying phase_shift = np.pi / 2, i.e. phase shift and sine phase compensate each other.