如何在 Python 中解释离散傅里叶变换 (FFT) 的结果

Question

关于这个主题有很多问题，我已经循环浏览了很多问题，以获得关于处理频率的概念指针（here 和 here）、关于 numpy 函数的文档（here）、操作方法信息关于提取幅度和相位 (here)，并走出站点，例如 this 或 this。

但是，只有通过简单的示例向自己“证明”这一痛苦并检查不同函数的输出（与其手动实现形成对比），才让我有了一些想法。

答案试图用 Python 记录和分享与 DFT 相关的细节，如果不简单解释，这些细节可能会构成进入障碍。

Answer 1

DFT（FFT 是它的算法计算）是模拟信号 s(t) 的有限离散样本 N 之间的点积（函数时间或空间）和一组复指数的基向量（sin 和 cos 函数）。尽管样本自然是有限的并且可能没有周期性，但它被隐含地认为是周期性重复的离散函数。即使在处理实值信号（通常情况）时，使用复数（欧拉方程）也很方便。在带有np.fft.fft(s) 的信号上实现该函数可能只是为了获得复数的输出系数并陷入其解释中，这可能是令人生畏的。一些步骤是必不可少的：

复指数中的频率是多少？

1. DFT 不一定以赫兹为单位保留采样频率。频率是指数 (k)。
2. 索引 k 的范围从 0 to N - 1 可以被认为具有 cycles / set 的单位（该集是信号 @ 的 N 样本987654333@)。我将省略讨论奈奎斯特极限，但对于真实信号，频率在 N / 2 之后形成镜像，并在该点之后作为负递减值给出（在隐式周期性框架内不是问题） . FFT 中使用的频率不仅仅是k，而是k / N，被认为具有cycles / sample 的单位。见this reference。示例 (reference)：如果对信号进行采样 N = 5 次，则频率为：np.fft.fftfreq(5)，产生[ 0 , 0.2, 0.4, -0.4, -0.2]，即[0/5, 1/5, 2/5, -2/5, -1/5]。
3. 要将这些频率转换为有意义的单位（例如赫兹或毫米），上述周期/样本中的值需要除以采样间隔 T（例如样本之间的距离（以秒为单位））。继续上面的示例，有一个内置调用：np.fft.fftfreq(5, d=T)：如果模拟信号s 以等距间隔T = 1/2 秒采样5 次，总采样为NT = 5 x 1/2 sec，则归一化频率将是np.fft.fftfreq(5, d = 1/2)，产生[0 0.4 0.8 -0.8 -0.4] 或[0/NT, 1/NT, 2/NT, -2/NT, -1/NT]。
4. 归一化或非归一化频率用于控制角频率 (ω_m)，表示为ω_m = 2π k/NT。请注意，NT 是总持续时间 信号被采样的。索引k 确实会产生对应于k = 1 的基频（ω-naught）的倍数 - 完成的（余）正弦波的频率 正好在NT (here) 上振荡一次。

FFT中系数的幅度、频率和相位

1. 给定 FFT S = fft.fft(s) 的输出，输出系数 (here) 的 幅度 只是输出系数中复数的欧几里德范数，已针对实数对称性进行了调整信号（x 2）和样本数量1/N：magnitudes = 1/N * np.abs(S)
2. 频率与上述np.fft.fftfreq(N) 的调用相匹配，或者更方便地结合实际的模拟频率单位frequencies = np.fft.fftfreq(N, d=T)。
3. 每个系数的相位是极坐标形式的复数的角度phase = np.arctan(np.imag(S)/np.real(S))

---

如何在FFT中找到信号s中的主导频率及其系数？

1. 撇开绘图，找到与最高幅度频率对应的索引k可以完成为index = np.argmax(np.abs(S))。例如，要查找幅度最大的4 索引，调用是indices = np.argpartition(S,-4)[-4:]。
2. 并找到实际对应的系数：S[index]，频率为freq_max = np.fft.fftfreq(N, d=T)[index]。

---

得到系数后再现原始信号：

通过正弦和余弦再现s（here 中的第 150 页）：

    Re = np.real(S[index])
    Im = np.imag(S[index])
    
    s_recon = Re * 2/N * np.cos(-2 * np.pi * freq_max * t) + abs(Im) * 2/N * np.sin(-2 * np.pi * freq_max * t) 

这是一个完整的例子：

import numpy as np
import matplotlib.pyplot as plt

N = 10000           # Sample points     
T = 1/5000          # Spacing
# Total duration N * T= 2
t = np.linspace(0.0, N*T, N, endpoint=False) # Time: Vector of 10,000 elements from 0 to N*T=2.
frequency = np.fft.fftfreq(t.size, d=T)      # Normalized Fourier frequencies in spectrum.

f0 = 25             # Frequency of the sampled wave
phi = np.pi/8       # Phase
A = 50              # Amplitude

s = A * np.cos(2 * np.pi * f0 * t + phi) # Signal

S = np.fft.fft(s)   # Unnormalized FFT

index = np.argmax(np.abs(S))
print(S[index])
magnitude = np.abs(S[index]) * 2/N
freq_max = frequency[index]

phase = np.arctan(np.imag(S[index])/np.real(S[index]))
print(f"magnitude: {magnitude}, freq_max: {freq_max}, phase: {phase}")
print(phi)

fig, [ax1,ax2] = plt.subplots(nrows=2, ncols=1, figsize=(10, 5))
ax1.plot(t,s, linewidth=0.5, linestyle='-', color='r', marker='o', markersize=1,markerfacecolor=(1, 0, 0, 0.1))  
ax1.set_xlim([0, .31])
ax1.set_ylim([-51,51])
ax2.plot(frequency[0:N//2], 2/N * np.abs(S[0:N//2]), '.', color='xkcd:lightish blue', label='amplitude spectrum')
plt.xlim([0, 100])
plt.show()

Re = np.real(S[index])
Im = np.imag(S[index])

s_recon = Re*2/N * np.cos(-2 * np.pi * freq_max * t) + abs(Im)*2/N * np.sin(-2 * np.pi * freq_max * t)

fig = plt.figure(figsize=(10, 2.5))

plt.xlim(0,0.3)
plt.ylim(-51,51)
plt.plot(t,s_recon, linewidth=0.5, linestyle='-', color='r', marker='o', markersize=1,markerfacecolor=(1, 0, 0, 0.1))  
plt.show()

s.all() == s_recon.all()