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

Question

关于这个主题有很多问题，我已经循环浏览了其中很多问题，获取了有关处理频率的概念性指针（此处和此处）、有关 numpy 函数的文档（此处）、如何提取幅度和相位的信息（此处）、并走出网站，例如this或this。

然而，只有痛苦地用简单的例子向自己“证明它”，并检查不同函数的输出与它们的手动实现相比，才给了我一些想法。

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

Answer 1

\n

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

\n

复指数中的频率是多少？

\n

\n
1. DFT 不一定保留以赫兹为单位的采样频率。频率是指数（ k )。
\n
2. 索引k 的范围为0 to N - 1，并且可以被认为具有周期/集合（该集合是N信号的样本s）的单位。我将省略讨论奈奎斯特极限，但对于真实信号，频率在N / 2之后形成镜像，并在该点之后以负递减值给出（不是隐式周期性框架内的问题）。FFT 中使用的频率不仅仅是k，而是k/N，单位为周期/样本。请参阅此参考资料。示例（参考）：如果对信号进行采样，则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]。
\n
3. 要将这些频率转换为有意义的单位（例如赫兹或毫米），需要将上述周期/样本中的值除以采样间隔T（例如样本之间的距离（以秒为单位））。继续上面的示例，有一个内置调用：np.fft.fftfreq(5, d=T)：如果对模拟信号以等距间隔secs进行采样，总样本为，则归一化频率将为，产生或。5T = 1/2NT = 5 x 1/2 secnp.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]
\n
4. 归一化或非归一化频率用于控制角频率( \xcf\x89_m )，表示为\xcf\x89_m = 2\xcf\x80 k/NT。请注意，这NT是对信号进行采样的总持续时间。该索引k确实会导致基频（\xcf\x89-naught）的倍数，对应于在（此处k = 1）完成\n一次振荡的（余）正弦波的频率。NT
\n

\n

\n

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

\n

\n
1. 给定 FFT 的输出S = fft.fft(s)，输出系数的大小（此处）只是输出系数中复数的欧几里得范数，x 2并根据实际信号 ( ) 的对称性和样本数进行调整1/N：magnitudes = 1/N * np.abs(S)
\n
2. 频率与上面解释的呼叫相匹配np.fft.fftfreq(N)，或者更方便地结合实际的模拟频率单位frequencies = np.fft.fftfreq(N, d=T)。
\n
3. 每个系数的相位是复数极坐标形式的角度phase = np.arctan(np.imag(S)/np.real(S))
\n

\n\n

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

\n

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

\n\n

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

\n

s通过正弦和余弦再现（此处第 150 页）：

\n

    Re = np.real(S[index])\n    Im = np.imag(S[index])\n    \n    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) \n

\n

这是一个完整的例子：

\n

import numpy as np\nimport matplotlib.pyplot as plt\n\nN = 10000           # Sample points     \nT = 1/5000          # Spacing\n# Total duration N * T= 2\nt = np.linspace(0.0, N*T, N, endpoint=False) # Time: Vector of 10,000 elements from 0 to N*T=2.\nfrequency = np.fft.fftfreq(t.size, d=T)      # Normalized Fourier frequencies in spectrum.\n\nf0 = 25             # Frequency of the sampled wave\nphi = np.pi/8       # Phase\nA = 50              # Amplitude\n\ns = A * np.cos(2 * np.pi * f0 * t + phi) # Signal\n\nS = np.fft.fft(s)   # Unnormalized FFT\n\nindex = np.argmax(np.abs(S))\nprint(S[index])\nmagnitude = np.abs(S[index]) * 2/N\nfreq_max = frequency[index]\n\nphase = np.arctan(np.imag(S[index])/np.real(S[index]))\nprint(f"magnitude: {magnitude}, freq_max: {freq_max}, phase: {phase}")\nprint(phi)\n\nfig, [ax1,ax2] = plt.subplots(nrows=2, ncols=1, figsize=(10, 5))\nax1.plot(t,s, linewidth=0.5, linestyle=\'-\', color=\'r\', marker=\'o\', markersize=1,markerfacecolor=(1, 0, 0, 0.1))  \nax1.set_xlim([0, .31])\nax1.set_ylim([-51,51])\nax2.plot(frequency[0:N//2], 2/N * np.abs(S[0:N//2]), \'.\', color=\'xkcd:lightish blue\', label=\'amplitude spectrum\')\nplt.xlim([0, 100])\nplt.show()\n\nRe = np.real(S[index])\nIm = np.imag(S[index])\n\ns_recon = Re*2/N * np.cos(-2 * np.pi * freq_max * t) + abs(Im)*2/N * np.sin(-2 * np.pi * freq_max * t)\n\nfig = plt.figure(figsize=(10, 2.5))\n\nplt.xlim(0,0.3)\nplt.ylim(-51,51)\nplt.plot(t,s_recon, linewidth=0.5, linestyle=\'-\', color=\'r\', marker=\'o\', markersize=1,markerfacecolor=(1, 0, 0, 0.1))  \nplt.show()\n\ns.all() == s_recon.all()\n

\n