信息处理 - 为什么我的 DFT/FFT 总是 90° 异相？ - 吾爱随笔录

为什么我的 DFT/FFT 总是 90° 异相？

信息处理 fft 自由度 Python 阶段

2022-02-04 10:25:27

我在一台机器上使用 Python 和 Numpy 进行 FFT，在另一台机器上使用 C#。我正在使用一些虚拟数据来模拟我最终将如何从 C#/UWP 应用程序中的传感器收集数据。这两种方法在虚拟数据上产生一致的结果，这很棒。但是，正弦波频率处的相位始终为 -90°。

这是代码。忽略它本质上是如何做同样的事情 7 次，最终这些阵列将被来自 7 个不同传感器的数据填充。采样频率为 4096 Hz，信号收集时间为 1 秒，时间间隔为 [-32768,32768]。请注意，幅度和相位图上的 x 轴放大到中心频率正负数赫兹。

import numpy as np
import matplotlib.pyplot as plt
import math

centerf = 1350
span = 2*np.pi*centerf*np.linspace(0,4095,4096)/4096
testData = 32768*np.sin(span)
testData = np.array([math.trunc(x) for x in testData])

plt.close("all")
sensors = ["Channel 1", "Channel 2", "Channel 3", "Channel 4", "Channel 5", "Channel 6", "Channel 7"]

data = []
data.append(testData)
data.append(testData)
data.append(testData)
data.append(testData)
data.append(testData)
data.append(testData)
data.append(testData)

theTime = np.linspace(0,len(data[0]),len(data[0]))/len(data[0])

thisFFT = [np.fft.fft(x) for x in data]

f, a = plt.subplots(7, 1,figsize=(10,8),sharex=True)
f2,a2 = plt.subplots(7,1,figsize=(10,8),sharex=True)
f3,a3 = plt.subplots(7,1,figsize=(10,8),sharex=True)

mag = [np.abs(x)/2048 for x in thisFFT]
phase2 = [np.arctan2(x.imag,x.real)*180/np.pi for x in thisFFT]

mag = [x[:2048] for x in mag] # removed mirrored upper half
phase2 = [x[:2048] for x in phase2]

for x in mag:
    x[0] = 0 # remove DC component

[aa.plot(theTime,d,c='r',lw=0.5) for (aa,d) in zip(a,data)]
[aa.plot(d,c='g',lw=0.75) for (aa,d) in zip(a2,mag)]
[aa.plot(d,c='b',lw=0.75) for (aa,d) in zip(a3,phase2)]
for ch,ax in zip(sensors,a.flat):
    ax.set(ylabel=ch)
for ch,ax in zip(sensors,a2.flat):
    ax.set(ylabel=ch)
for ch,ax in zip(sensors,a3.flat):
    ax.set(ylabel=ch)
a.flat[0].set(title="Amplitude over 1s")
a2.flat[0].set(title="FFT Magnitude")
a3.flat[0].set(title="FFT Phase")
a.flat[5].set(xlabel='Time')

[aa.axvline(x=centerf,lw=0.25) for aa in a2]
[aa.axvline(x=centerf,lw=0.25) for aa in a3]
[aa.set_xlim([centerf-20,centerf+20]) for aa in a2]
[aa.set_xlim([centerf-20,centerf+20]) for aa in a3]

f.tight_layout()
f2.tight_layout()
f3.tight_layout()

print("Magnitude at centerf: {c:1.3f}.".format(c=mag[0][centerf]))
print("Phase at centerf: {c:1.3f}°.".format(c=phase2[0][centerf]))

打印语句的输出是

Magnitude at centerf: 32767.397.
Phase at centerf: -90.000°.

幅度很好，图表总是在它应该在的地方显示一个很好的峰值。

我的问题：由于输入数据是始终以 sin(0)=0 开始的正弦波，并且 FFT 算法使用矩形窗口，因此基频处的相位不应该是 0°吗？相位有什么我不理解的吗？

2个回答

傅里叶变换（而 DFT 是具有有限长度序列的离散形式）只是输入信号与每个频率区间的相关性，具有连续形式：

X (ω) = \int_{t = 0}^{T} x (t) e^{- j ω t} d t

$X(\omega) = \int_{t=0}^Tx(t)e^{-j\omega t}dt$

（如果我们考虑有限的持续时间，但是我将继续使用傅里叶变换的连续时间版本，因为它将简化解释）。

在这里，我们看到对于任何频率 $\omega$ 我们正在做一个相关性 $x(t)$ 和 $e^{j\omega t}$ （相关是一个复共轭乘法，然后是积分，或者离散时的求和）。

作为相关性，我们看到当 x(t) 与 $e^{j\omega t}$ 在给定的时间间隔内，当它完全不相关（正交）时，类似地为 0。

现在记下欧拉的恒等式：

e^{- j ω t} = c o s (ω t) - j s i n (ω t)

$e^{-j \omega t} = cos(\omega t) - j sin(\omega t)$

在我们的第一个公式中替换它，我们可以将傅里叶变换扩展为两个相关性：

X (ω) = \int_{t = 0}^{T} x (t) c o s (ω t) d t - j \int_{t = 0}^{T} x (t) s i n (ω t) d t

$X(\omega) = \int_{t=0}^T x(t)cos(\omega t)dt - j\int_{t=0}^T x(t)sin(\omega t)dt$

由此我们看到，当 $x(t) =sin(\omega t)$ , （并且持续时间 T 涵盖一个完整周期或整数个周期），第二项将是完全相关的，而第一项将为零，反之亦然 $x(t) =cos(\omega t)$ . 如果x(t) 为实数，则与第二项的任何相关性都会导致虚构结果。x(t) 上的任何相移，如 $x(t) = sin(\omega t + \phi)$ 在哪里 $\phi$ 除了 $\pm N\pi/2$ 对于 N 任何整数，将导致一个复杂的结果 $X(\omega)$ 有实部和虚部。

请注意，它的离散形式 DFT 如下所示，其中可以解释相同的结果，但对于一些不太熟悉它的人来说可能不太容易看到。

X [k] = Σ_{n = 0}^{N - 1} x [n] e^{- j n ω_{o} \frac{k}{N}}

$X[k] = \Sigma_{n=0}^{N-1}x[n]e^{-jn\omega_o \frac{k}{N}}$

这是一个常见的误解，即典型的 FFT 实现显示相对于 sin() 函数的相位，该函数从 FFT 孔径的“左”边缘零开始。

但是 sin() 函数对应于复指数的虚部。典型的默认设置是严格实数 FFT 结果的相位为零。

cos() 函数也产生一个正弦波，但它是一个正弦波，因此它被定位（或偏移）为围绕 FFT 孔径的中心对称（例如，一个偶数函数而不是一个奇数函数）。cos() 函数对应于复指数的实部。

如果您执行 atan2(x.imag,x.real) 实部参数为零且虚部参数非零，您将获得 +-90 度的相位。

其它你可能感兴趣的问题

上一篇减少递归滤波器启动瞬态的技术参考下一篇在哪里可以了解“模拟原型滤波器”？