如何用Python解释离散傅里叶变换(FFT)的结果
如何用Python解释离散傅里叶变换(FFT)的结果
提问于 2022-01-22 06:55:49
回答 1
Stack Overflow用户
发布于 2022-01-22 06:55:49
DFT (FFT算法计算)是模拟信号N的有限离散样本数s(t) (时间或空间函数)与复指数基向量(sin和cos函数)之间的一个点乘积。虽然样本自然是有限的，并且可能没有周期性，但它被隐式地认为是一个周期性重复的离散函数。即使在处理实值信号(通常情况下)时，也可以方便地处理复数(Euler方程)。用np.fft.fft(s)在信号上实现这个函数，只会得到复数的输出系数，并使其无法解释，这可能是很可怕的。一些步骤是必不可少的：
复指数中的频率是多少？
- DFT不一定保持赫兹的采样频率。频率是指数(，k，)。
- 指数k范围从
0 to N - 1，可以认为是具有循环/集的单位(集合是信号s的N样本)。我将省略讨论Nyquist极限，但对于实际信号，频率在N / 2后形成镜像，并在此之后作为负递减值给出(不是隐式周期性框架内的问题)。在快速傅立叶变换中使用的频率不是简单的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]。 - 为了将这些频率转换成有意义的单位(例如，Hetz或mm)，上述循环/样本中的值需要除以采样间隔T (例如样本之间的距离)。继续上面的例子，有一个内置的调用：
np.fft.fftfreq(5, d=T)：如果模拟信号s被以等距的5时间间隔T = 1/2秒采样一个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]。 - 采用归一化或非归一化频率控制角频率 (ω_m)，表示为ω_m = 2πk/NT。请注意，
NT是采样信号的总持续时间。指数k确实导致基频的倍数(ω-naught)对应于k = 1--即(co-)正弦波的频率，正好在NT(这里)上完成一次振荡。
FFT中系数的大小、频率和相位
- 给定FFT
S = fft.fft(s)的输出，输出系数(这里)的幅值只是输出系数中的复数的欧氏范数，这些复数是根据实信号(x 2)的对称性和样本数(1/N：magnitudes = 1/N * np.abs(S))而调整的。 - 这些频率与上述
np.fft.fftfreq(N)解释的呼叫相匹配，或者更方便地结合实际的模拟频率单元，frequencies = np.fft.fftfreq(N, d=T)。 - 每个系数的相位是极型
phase = np.arctan(np.imag(S)/np.real(S))中复数的角度。
如何找到FFT中信号s中的主频率及其系数？
- 除作图外，找出与频率最高的频率对应的指数
k，可作为index = np.argmax(np.abs(S))完成。例如，要找到规模最大的4索引，调用是indices = np.argpartition(S,-4)[-4:]。 - 并求出实际的对应系数：
S[index]和频率freq_max = np.fft.fftfreq(N, d=T)[index]。
在获得系数后再现原始信号：
通过正弦和余弦复制s (第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()页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/70810666
复制相关文章
相似问题
腾讯云开发者
Copyright © 2013 - 2026 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有
深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 
粤公网安备44030502008569号
腾讯云计算（北京）有限责任公司 京ICP证150476号 | 京ICP备11018762号