Scipy FFT - 如何获得相位角

Question

我在使用 python 中的 scipy fft 模块获取简单正弦曲线的相位时遇到了麻烦。我密切关注this 教程并将matlab 代码转换为python。但是，无论我将哪个阶段用于输入，图表始终显示 3。我错过了什么？

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack
import cmath
A=10
fc = 10
phase=60
fs=32#Sampling frequency with oversampling factor of 32

t = np.arange(0,2,1/fs)

#Convert the phase shift to radians from degrees.
phi = phase*np.pi/180

x=A*np.cos(2*np.pi*fc*t+phi)

N=256
X = scipy.fftpack.fftshift(scipy.fftpack.fft(x,N))/N

df=fs/N #Frequency resolution.
sampleindex = np.arange(-N/2,N/2,1) #Ordered index for FFT plot.
f = sampleindex*df #x-axis index continued to ordered frequencies

raw_phases = np.angle(X)

X2=np.copy(X)#Store the FFT results in another array.
#Detect very small numbers and ignore them.
tau = max(abs(X))/10
X2[abs(X)<tau]=0

phase=[cmath.phase(i) for i in X2]
plt.plot(f,phase)
plt.show()

编辑：这是一些更简单的代码。似乎仍然无法获得相位。

y = 28*np.sin(2*np.pi/7*x)
yf = scipy.fftpack.fft(y)
xf = np.linspace(0.0, 1.0/(2.0*T), N/2)
phase = np.angle(yf)
yf = np.abs(yf[:N//2])
phase = phase[:N//2]
xf = xf[1:]
yf = yf[1:]
phase = phase[1:]
yf = yf-np.mean(yf)
#The frequencies seem to always be scaled by 0.1433, not sure why.
c = 2*np.pi/7/0.143301
freqs = xf[yf>np.std(yf)]*c
phases = phase[yf>np.std(yf)]

我得到的频率集中在 2*np.pi/7 附近。但我得到的阶段是：

array([-0.217795  , -0.22007488, -0.22226087,  2.91723935,  2.91524011,
    2.91333367])

虽然根本不应该有阶段。

Answer 1

FFT 测量循环相位，参考输入数据窗口的开头和结尾。如果您的输入正弦波在 FFT 孔径中不是完全整数周期，那么窗口开始和结束的相位之间将存在不连续性，因此 FFT 相位测量将不是您所期望的。

相反，我建议在数据窗口中间（N/2）的所需相位开始（或参考）输入正弦波，而不是开始。然后在 FFT 之前做一个循环 fftshift。然后，生成的 FFT 相位测量将表示相对于原始数据窗口中间的相位，其中没有不连续性；和 atan2() 将表示您的连续输入波形的偶数和奇数分量之间的比率，正如通常预期的那样。

Answer 2

这是显示如何获取角度的最简单的代码。

请注意，我创建的信号 y 中包含整数个句点（正如我建议的 in a comment 和 @hotpaw2 建议的 in their answer）。

这是 OP 代码的问题。

我使用linspace 创建时间轴，使用endpoint=False 选项。这很重要，如果t 包含数字 10，那么我们将不再有精确的整数个句点。使用离散傅立叶变换时，考虑重复信号时会发生什么会有所帮助。只需获取y，并连接其自身的副本：np.hstack((y,y))。这个新信号仍然是单个正弦波的采样，还是我们创造了更复杂的东西？两个副本相遇时会发生什么？

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack

phase = np.pi / 4
t = np.linspace(0, 10, num=200, endpoint=False)
y = np.cos(2 * np.pi * t + phase)
Y = scipy.fftpack.fftshift(scipy.fftpack.fft(y))
f = scipy.fftpack.fftshift(scipy.fftpack.fftfreq(len(t)))

p = np.angle(Y)
p[np.abs(Y) < 1] = 0
plt.plot(f, p)
plt.show()