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 测量循环相位，参考输入数据的最开始和最结束 window。如果您的输入正弦波在 FFT 孔径中不完全是整数周期，那么 window 开始和结束的相位之间将存在不连续性，因此 FFT 相位测量将不是您可能的结果期待。

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

Answer 2

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

请注意，我创建了信号 y，其中包含整数个周期（正如我所建议的 , and @hotpaw2 suggested ）。

这是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()

Scipy FFT - 如何获得相位角

Scipy FFT - how to get phase angle

python

fft

scipy