NumPy FFT 教程：快速入门指南 – wiki基地

作者： /

NumPy FFT 教程：快速入门指南

摘要

快速傅里叶变换（Fast Fourier Transform, FFT）是一种计算离散傅里叶变换（Discrete Fourier Transform, DFT）及其逆变换的高效算法。DFT 是数字信号处理和许多其他科学计算领域的基础工具，它能够将一个信号（通常是时间序列）分解为其组成频率的频谱。NumPy 库提供了一个强大的 numpy.fft 模块，使得在 Python 中执行 FFT 计算变得简单而高效。本教程将详细介绍 NumPy FFT 的基础知识、核心函数、实际应用示例以及结果解读，旨在帮助您快速掌握使用 NumPy 进行频域分析的基本技能。

目录

  1. 引言：傅里叶变换与 FFT 的重要性
    • 什么是傅里叶变换？（直观理解）
    • 什么是离散傅里叶变换（DFT）？
    • 什么是快速傅里叶变换（FFT）？
    • 为什么需要 FFT？（应用场景）
  2. NumPy FFT 模块 (numpy.fft) 概览
    • 安装与导入 NumPy
    • 核心函数介绍
  3. 核心函数详解与基础示例
    • np.fft.fft()：计算一维 DFT
      • 示例 1：简单信号的 FFT
      • 理解 FFT 输出：复数、对称性
    • np.fft.fftfreq()：计算频率轴
      • 示例 2：为 FFT 结果生成频率坐标
    • np.fft.fftshift() 与 np.fft.ifftshift()：频率排序调整
      • 为什么需要 fftshift？
      • 示例 3：可视化频谱
    • np.fft.ifft()：计算一维逆 DFT
      • 示例 4：从频域恢复时域信号
  4. 频谱分析实战
    • 示例 5：分析包含多种频率成分的信号
    • 示例 6：分析带噪声的信号
    • 解读频谱图：幅度和相位
  5. 进阶话题与注意事项
    • 实数信号的优化：np.fft.rfft() 和 np.fft.irfft()
    • 采样率与奈奎斯特频率
    • 频谱泄漏与窗函数（简介）
    • 幅度归一化
    • 多维 FFT（fft2, ifft2, fftn, ifftn）
  6. 总结与展望

1. 引言：傅里叶变换与 FFT 的重要性

什么是傅里叶变换？（直观理解）

想象一下你听到一段音乐，比如一个和弦。你的耳朵（和大脑）能够分辨出这个和弦是由哪些不同的音符（即不同的音高或频率）组成的。傅里叶变换在数学上做的就是类似的事情：它是一种将一个复杂的信号（如随时间变化的声音波形）分解成一系列简单的正弦波（代表不同频率和幅度）之和的方法。就像棱镜可以将白光分解成彩虹的各种颜色（不同频率的光波）一样，傅里叶变换可以将一个时域信号转换到频域，揭示其内在的频率结构。

什么是离散傅里叶变换（DFT）？

在现实世界中，我们通常处理的是通过采样得到的离散数据点，而不是连续的函数。例如，数字音频就是以固定的时间间隔（采样率）记录声音的幅度。对于这种离散的、有限长度的信号，我们使用离散傅里叶变换（DFT）。DFT 将 N 个离散的时域数据点转换为 N 个离散的频域数据点，这些频域数据点描述了信号中特定频率分量的幅度和相位。

什么是快速傅里叶变换（FFT）？

直接计算 DFT 的复杂度大约是 O(N²)，其中 N 是数据点的数量。当 N 很大时，这个计算量会变得非常巨大。快速傅里叶变换（FFT）并非一种新的变换，而是一系列高效计算 DFT 的算法。最著名的 FFT 算法（如 Cooley-Tukey 算法）可以将计算复杂度显著降低到 O(N log N)。这使得对包含数百万甚至更多数据点的大型数据集进行频域分析成为可能。NumPy 的 numpy.fft 模块内部就实现了这些高效的 FFT 算法。

为什么需要 FFT？（应用场景）

FFT 的应用极其广泛，贯穿于科学和工程的各个领域：

  • 信号处理：分析音频信号的频谱、去除噪声、滤波、调制解调。
  • 图像处理：图像压缩（如 JPEG）、图像滤波、模式识别、特征提取（如频率域特征）。
  • 通信系统：OFDM（正交频分复用）等现代通信技术的基础。
  • 数据压缩：通过保留主要的频率分量并丢弃次要分量来压缩数据。
  • 物理学与工程学：振动分析、声学分析、天文学信号分析、求解偏微分方程。
  • 金融：分析时间序列数据中的周期性模式。

基本上，任何需要理解数据中周期性或频率成分的场景，FFT 都可能派上用场。

2. NumPy FFT 模块 (numpy.fft) 概览

安装与导入 NumPy

如果你还没有安装 NumPy，可以使用 pip 进行安装：

bash
pip install numpy

在 Python 脚本或交互式环境中使用 NumPy FFT 功能，首先需要导入 NumPy 库，通常约定俗成地导入为 np：

python
import numpy as np
import matplotlib.pyplot as plt # 导入绘图库，用于可视化

NumPy 的 FFT 功能都包含在 numpy.fft 子模块中。

核心函数介绍

numpy.fft 模块提供了多种函数，以下是最常用的一些：

  • np.fft.fft(a, n=None, axis=-1, norm=None)：计算一维离散傅里叶变换 (DFT)。
  • np.fft.ifft(a, n=None, axis=-1, norm=None)：计算一维逆离散傅里叶变换 (IDFT)。
  • np.fft.fftfreq(n, d=1.0)：返回 DFT 采样频率。给定窗口长度 n 和采样间距 d，它会返回对应的频率轴坐标。
  • np.fft.fftshift(x, axes=None)：将 FFT 输出中零频率分量移动到频谱中心。这对于可视化非常有用。
  • np.fft.ifftshift(x, axes=None)：fftshift 的逆操作。
  • np.fft.rfft(a, n=None, axis=-1, norm=None)：计算实数输入信号的一维 DFT，利用实数信号 FFT 结果的共轭对称性，只返回非负频率部分，计算更快，占用内存更少。
  • np.fft.irfft(a, n=None, axis=-1, norm=None)：rfft 的逆变换。
  • np.fft.rfftfreq(n, d=1.0)：为 rfft 的输出计算频率轴。
  • np.fft.fft2(a, s=None, axes=(-2, -1), norm=None)：计算二维 DFT（常用于图像）。
  • np.fft.ifft2(a, s=None, axes=(-2, -1), norm=None)：计算二维 IDFT。
  • np.fft.fftn(a, s=None, axes=None, norm=None)：计算 N 维 DFT。
  • np.fft.ifftn(a, s=None, axes=None, norm=None)：计算 N 维 IDFT。

本教程将重点关注一维 FFT 相关的核心函数：fft, ifft, fftfreq, fftshift。

3. 核心函数详解与基础示例

np.fft.fft()：计算一维 DFT

np.fft.fft(a) 接受一个数组 a（通常是一维时间序列数据），并返回其 DFT 结果。返回的结果是一个复数数组，长度与输入数组 a 相同。

示例 1：简单信号的 FFT

让我们从一个非常简单的信号开始：一个直流信号（常数）。

“`python

示例 1: 直流信号的 FFT

N = 100 # 信号长度（采样点数）
signal_dc = np.ones(N) # 幅度为 1 的直流信号

fft_dc = np.fft.fft(signal_dc)

print(“FFT of DC signal (first few elements):”)
print(fft_dc[:5])
print(“\nShape of FFT output:”, fft_dc.shape)

计算频率轴（稍后详细解释）

假设采样间隔为 1 (即采样率为 1 Hz)

freq_dc = np.fft.fftfreq(N, d=1.0)

可视化幅度谱

plt.figure(figsize=(10, 4))
plt.plot(freq_dc, np.abs(fft_dc)) # 使用 np.abs() 获取幅度
plt.title(‘FFT Magnitude of DC Signal’)
plt.xlabel(‘Frequency (Hz)’)
plt.ylabel(‘Magnitude’)
plt.grid(True)
plt.show()
“`

理解 FFT 输出：

  • 复数：fft 函数的输出是一个复数数组。每个复数包含两个信息：
    • 幅度 (Magnitude)：表示该频率分量的强度或能量。通常通过 np.abs() 计算。
    • 相位 (Phase)：表示该频率分量的相对于原点的起始位置（时间偏移）。通常通过 np.angle() 计算。对于许多应用，幅度信息更为重要。
  • 输出顺序与对称性：对于实数输入信号（大部分情况），FFT 的输出具有共轭对称性。fft_output[k] 和 fft_output[N-k] 是共轭复数（实部相同，虚部相反）。这意味着幅度谱 (np.abs(fft_output)) 是对称的。
    • fft_output[0]：对应于零频率（DC 分量）。在上面的例子中，我们看到第一个元素是 (100+0j)，幅度为 100（等于信号值 1 乘以信号长度 N），其余元素接近于零。这是因为纯直流信号只包含零频率成分。
    • fft_output[1] 到 fft_output[N//2]：对应于正频率分量。
    • fft_output[N//2+1] 到 fft_output[N-1]：对应于负频率分量（根据 DFT 定义排列）。

np.fft.fftfreq()：计算频率轴

fft 的输出数组的索引 k 并不直接等于频率。我们需要使用 np.fft.fftfreq(n, d) 来计算每个索引对应的实际频率。

  • n：窗口长度（FFT 输入/输出的点数，通常是信号长度 N）。
  • d：采样间隔（时间单位）。如果已知采样率 fs（样本/秒），则 d = 1 / fs。默认 d=1.0。

该函数返回一个频率数组，其顺序与 fft 输出的顺序相匹配：先是零频率，然后是正频率，最后是负频率（从最高负频率开始）。

示例 2：为 FFT 结果生成频率坐标

让我们分析一个正弦波信号。

“`python

示例 2: 单一频率正弦波的 FFT 及频率轴

fs = 100 # 采样率 (Hz)
T = 2 # 信号持续时间 (seconds)
N = fs * T # 总采样点数

t = np.linspace(0.0, T, N, endpoint=False) # 时间向量 (不包含 T)

frequency1 = 5 # 正弦波频率 (Hz)
amplitude1 = 1.0
signal_sine = amplitude1 * np.sin(2 * np.pi * frequency1 * t)

计算 FFT

fft_sine = np.fft.fft(signal_sine)

计算频率轴

d = 采样间隔 = 1 / fs

freq_sine = np.fft.fftfreq(N, d=1/fs)

print(“Shape of frequency axis:”, freq_sine.shape)
print(“Frequency bins (first few and last few):”)
print(np.round(freq_sine[:5], 2))
print(np.round(freq_sine[-5:], 2))

可视化幅度谱

plt.figure(figsize=(12, 5))
plt.plot(freq_sine, np.abs(fft_sine))
plt.title(‘FFT Magnitude of 5 Hz Sine Wave (Unshifted)’)
plt.xlabel(‘Frequency (Hz)’)
plt.ylabel(‘Magnitude’)
plt.grid(True)
plt.show()
“`

你会看到频谱图在 +5 Hz 和 -5 Hz 处有两个显著的峰值。这正是因为实数信号的 FFT 结果是对称的。峰值的高度与信号幅度和长度有关（稍后讨论归一化）。

np.fft.fftshift() 与 np.fft.ifftshift()：频率排序调整

fftfreq 返回的频率数组顺序是 [0, f_1, ..., f_max, -f_max, ..., -f_1]。这对于计算是自然的，但对于人类观察和绘图来说，我们通常更习惯看到零频率在中间，两侧分别是负频率和正频率，即 [-f_max, ..., -f_1, 0, f_1, ..., f_max] 的顺序。

np.fft.fftshift(x) 函数就是用来实现这种重新排序的。它将数组 x 的前半部分和后半部分进行交换。对于 FFT 结果和 fftfreq 的结果同时使用 fftshift，就可以得到以零频率为中心的频谱图。np.fft.ifftshift(x) 是其逆操作。

示例 3：可视化频谱

“`python

示例 3: 使用 fftshift 可视化频谱

— 复用 示例 2 的 signal_sine, fft_sine, freq_sine —

fs = 100
T = 2
N = fs * T
t = np.linspace(0.0, T, N, endpoint=False)
frequency1 = 5
amplitude1 = 1.0
signal_sine = amplitude1 * np.sin(2 * np.pi * frequency1 * t)
fft_sine = np.fft.fft(signal_sine)
freq_sine = np.fft.fftfreq(N, d=1/fs)

—

对 FFT 结果和频率轴都应用 fftshift

fft_sine_shifted = np.fft.fftshift(fft_sine)
freq_sine_shifted = np.fft.fftshift(freq_sine)

print(“Shifted frequency bins (first few and last few):”)
print(np.round(freq_sine_shifted[:5], 2))
print(np.round(freq_sine_shifted[-5:], 2))

可视化移位后的幅度谱

plt.figure(figsize=(12, 5))
plt.plot(freq_sine_shifted, np.abs(fft_sine_shifted))
plt.title(‘FFT Magnitude of 5 Hz Sine Wave (Shifted)’)
plt.xlabel(‘Frequency (Hz)’)
plt.ylabel(‘Magnitude’)
plt.xlim(-fs/2, fs/2) # 通常将 x 轴限制在奈奎斯特频率范围内
plt.grid(True)
plt.show()
“`

现在，频谱图的中心是 0 Hz，两个峰值清晰地位于 -5 Hz 和 +5 Hz。这更符合我们对频谱的直观理解。

np.fft.ifft()：计算一维逆 DFT

np.fft.ifft(a) 函数可以从频域表示 a（通常是 fft 的结果）恢复出原始的时域信号。这是傅里叶变换可逆性的体现。

示例 4：从频域恢复时域信号

“`python

示例 4: 使用 ifft 恢复信号

— 复用 示例 2 的 signal_sine 和 fft_sine —

fs = 100
T = 2
N = fs * T
t = np.linspace(0.0, T, N, endpoint=False)
frequency1 = 5
amplitude1 = 1.0
signal_sine = amplitude1 * np.sin(2 * np.pi * frequency1 * t)
fft_sine = np.fft.fft(signal_sine)

—

计算逆 FFT

recovered_signal = np.fft.ifft(fft_sine)

检查恢复的信号

由于浮点计算误差，结果可能包含微小的虚部，对于实数输入，我们通常取其实部

recovered_signal_real = recovered_signal.real

比较原始信号和恢复信号

plt.figure(figsize=(12, 6))

plt.subplot(2, 1, 1)
plt.plot(t, signal_sine)
plt.title(‘Original Sine Wave Signal’)
plt.xlabel(‘Time (s)’)
plt.ylabel(‘Amplitude’)
plt.grid(True)

plt.subplot(2, 1, 2)
plt.plot(t, recovered_signal_real)
plt.title(‘Recovered Signal using IFFT’)
plt.xlabel(‘Time (s)’)
plt.ylabel(‘Amplitude’)
plt.grid(True)

plt.tight_layout()
plt.show()

检查误差

error = np.max(np.abs(signal_sine – recovered_signal_real))
print(f”Maximum absolute error between original and recovered signal: {error:.2e}”)
“`

你会看到恢复的信号几乎与原始信号完全一致，误差非常小，证明了 fft 和 ifft 的可逆性。

4. 频谱分析实战

现在我们将这些工具应用于更复杂的信号。

示例 5：分析包含多种频率成分的信号

创建一个包含两个不同频率和幅度的正弦波叠加的信号。

“`python

示例 5: 分析包含两种频率的信号

fs = 200 # 采样率 (Hz)
T = 1.0 # 信号持续时间 (seconds)
N = int(fs * T) # 总采样点数
t = np.linspace(0.0, T, N, endpoint=False)

信号成分

freq1 = 10 # Hz
amp1 = 2.0
freq2 = 30 # Hz
amp2 = 1.5
signal_multi = amp1 * np.sin(2np.pifreq1t) + amp2 * np.sin(2np.pifreq2t)

计算 FFT

fft_multi = np.fft.fft(signal_multi)
freq_multi = np.fft.fftfreq(N, d=1/fs)

Shift

fft_multi_shifted = np.fft.fftshift(fft_multi)
freq_multi_shifted = np.fft.fftshift(freq_multi)

计算幅度 (归一化)

对于正频率部分，幅度大约是 FFT 幅度绝对值 * 2 / N

对于 0 Hz (DC) 和奈奎斯特频率（如果 N 是偶数），幅度是 FFT 幅度绝对值 / N

magnitude_multi = np.abs(fft_multi_shifted) / N

修正非零频率的幅度 (乘以 2)

注意：这里简化处理，没有精确区分 0 和奈奎斯特频率

magnitude_multi[freq_multi_shifted != 0] *= 2

可视化

plt.figure(figsize=(12, 6))

plt.subplot(2, 1, 1)
plt.plot(t[:100], signal_multi[:100]) # 只显示部分时域波形
plt.title(‘Time Domain Signal (First 100 samples)’)
plt.xlabel(‘Time (s)’)
plt.ylabel(‘Amplitude’)
plt.grid(True)

plt.subplot(2, 1, 2)

plt.plot(freq_multi_shifted, np.abs(fft_multi_shifted)) # 未归一化幅度

plt.stem(freq_multi_shifted, magnitude_multi, markerfmt=’.’, basefmt=’k-‘, linefmt=’b-‘) # 使用 stem 图更清晰
plt.title(‘Frequency Domain – Magnitude Spectrum (Shifted & Normalized)’)
plt.xlabel(‘Frequency (Hz)’)
plt.ylabel(‘Approx. Amplitude’)
plt.xlim(-fs/2, fs/2)

plt.ylim(0, N/2 * max(amp1, amp2) * 1.1) # 调整 Y 轴范围观察峰值 (未归一化)

plt.ylim(0, max(amp1, amp2) * 1.1) # 调整 Y 轴范围观察峰值 (归一化后)
plt.grid(True)

plt.tight_layout()
plt.show()

找到幅度最大的几个频率

peak_indices = np.argsort(magnitude_multi)[-4:] # 幅度最大的4个点的索引 (因为对称性)
peak_freqs = freq_multi_shifted[peak_indices]
peak_mags = magnitude_multi[peak_indices]

print(“Detected peak frequencies (Hz):”)
print(np.round(peak_freqs, 1))
print(“Corresponding approximate amplitudes:”)
print(np.round(peak_mags, 1))
“`

频谱图清晰地显示了在 +10 Hz, -10 Hz, +30 Hz, -30 Hz 处的峰值。归一化后的峰值高度近似等于原始信号中对应频率成分的幅度 (2.0 和 1.5)。

示例 6：分析带噪声的信号

真实世界的信号往往夹杂着噪声。FFT 对于从噪声中识别主要频率成分非常有用。

“`python

示例 6: 分析带噪声的信号

— 复用 示例 5 的参数和信号 —

fs = 200
T = 1.0
N = int(fs * T)
t = np.linspace(0.0, T, N, endpoint=False)
freq1 = 10
amp1 = 2.0
freq2 = 30
amp2 = 1.5
signal_clean = amp1 * np.sin(2np.pifreq1t) + amp2 * np.sin(2np.pifreq2t)

—

添加高斯白噪声

noise_amplitude = 0.8
noise = noise_amplitude * np.random.randn(N)
signal_noisy = signal_clean + noise

计算 FFT

fft_noisy = np.fft.fft(signal_noisy)
freq_noisy = np.fft.fftfreq(N, d=1/fs)

Shift

fft_noisy_shifted = np.fft.fftshift(fft_noisy)
freq_noisy_shifted = np.fft.fftshift(freq_noisy)

计算幅度 (归一化)

magnitude_noisy = np.abs(fft_noisy_shifted) / N
magnitude_noisy[freq_noisy_shifted != 0] *= 2

可视化

plt.figure(figsize=(12, 6))

plt.subplot(2, 1, 1)
plt.plot(t[:200], signal_noisy[:200]) # 显示部分带噪时域波形
plt.title(‘Noisy Time Domain Signal (First 200 samples)’)
plt.xlabel(‘Time (s)’)
plt.ylabel(‘Amplitude’)
plt.grid(True)

plt.subplot(2, 1, 2)
plt.plot(freq_noisy_shifted, magnitude_noisy) # 使用 plot 更能显示噪声基底
plt.title(‘Frequency Domain – Magnitude Spectrum (Shifted & Normalized)’)
plt.xlabel(‘Frequency (Hz)’)
plt.ylabel(‘Approx. Amplitude’)
plt.xlim(-fs/2, fs/2)
plt.ylim(0, max(amp1, amp2) * 1.2) # 调整 Y 轴范围
plt.grid(True)

plt.tight_layout()
plt.show()

再次尝试找到峰值

注意：噪声会影响峰值检测，可能需要更复杂的算法

peak_indices_noisy = np.argsort(magnitude_noisy)[-4:]
peak_freqs_noisy = freq_noisy_shifted[peak_indices_noisy]
peak_mags_noisy = magnitude_noisy[peak_indices_noisy]

print(“Detected peak frequencies in noisy signal (Hz):”)
print(np.round(peak_freqs_noisy, 1))
print(“Corresponding approximate amplitudes:”)
print(np.round(peak_mags_noisy, 1))
“`

即使信号被噪声污染，频谱图中 10 Hz 和 30 Hz 的峰值依然清晰可见，远高于背景噪声的“噪声基底”。这展示了 FFT 在噪声背景下提取周期性信号的强大能力。

解读频谱图：幅度和相位

  • 幅度谱 (Magnitude Spectrum)：通常是我们最关心的，通过 np.abs(fft_output) 获得。它显示了信号中各个频率分量的强度。峰值的位置指示了信号中存在的主要频率，峰值的高度（经过适当归一化后）指示了该频率分量的幅度。
  • 相位谱 (Phase Spectrum)：通过 np.angle(fft_output) 获得。它表示每个频率分量的初始相位（相对于时间零点的偏移）。相位的解释通常比幅度更复杂，并且对信号的微小时间位移非常敏感。在许多入门级应用中，相位信息可能不是关注的重点，但在通信、滤波设计等领域至关重要。

“`python

计算并绘制相位谱 (示例 5 的信号)

phase_multi_shifted = np.angle(fft_multi_shifted) # 计算相位 (弧度)

plt.figure(figsize=(12, 5))
plt.plot(freq_multi_shifted, phase_multi_shifted)
plt.title(‘Phase Spectrum (Shifted)’)
plt.xlabel(‘Frequency (Hz)’)
plt.ylabel(‘Phase (radians)’)
plt.xlim(-fs/2, fs/2)
plt.ylim(-np.pi, np.pi)
plt.grid(True)
plt.show()
“`

你会看到相位谱通常看起来比较“杂乱”，尤其是在幅度很小的频率点。主要频率成分处的相位值才具有明确意义。

5. 进阶话题与注意事项

实数信号的优化：np.fft.rfft() 和 np.fft.irfft()

我们之前提到，对于实数输入信号，fft 的输出是共轭对称的。这意味着负频率部分的信息是冗余的。np.fft.rfft(a) 利用了这一点，它只计算并返回非负频率部分（从 0 Hz 到奈奎斯特频率）的 DFT 结果。

  • 优点：计算速度更快，返回的数组大小约为 fft 结果的一半，节省内存。
  • 输出：
    • 如果输入长度 N 是偶数，rfft 返回 N//2 + 1 个复数。
    • 如果输入长度 N 是奇数，rfft 返回 (N+1)//2 个复数。
  • 频率轴：需要使用 np.fft.rfftfreq(n, d) 来获取对应的频率。
  • 逆变换：使用 np.fft.irfft(a) 进行逆变换。irfft 会自动根据输入的复数数组重建共轭对称的完整频域表示，然后执行逆变换，确保得到实数输出。

“`python

使用 rfft 分析示例 5 的信号

signal_multi = amp1 * np.sin(2np.pifreq1t) + amp2 * np.sin(2np.pifreq2t) # N=200 (偶数)

rfft_multi = np.fft.rfft(signal_multi)
rfftfreq_multi = np.fft.rfftfreq(N, d=1/fs)

print(“Shape of rfft output:”, rfft_multi.shape) # 应为 N//2 + 1 = 101
print(“Shape of rfftfreq output:”, rfftfreq_multi.shape) # 应为 101

幅度归一化 (rfft)

magnitude_rfft = np.abs(rfft_multi) / N

对于 rfft, 除了 0 Hz 和奈奎斯特频率（如果存在），其他都需要乘以 2

检查奈奎斯特频率是否存在 (N 为偶数时，最后一个点是奈奎斯特频率)

nyquist_freq_present = N % 2 == 0
if nyquist_freq_present:
# 不乘 2 的索引是 0 和 N//2 (即最后一个索引)
indices_to_double = slice(1, -1)
else:
# 不乘 2 的索引只有 0
indices_to_double = slice(1, None)

magnitude_rfft[indices_to_double] *= 2

plt.figure(figsize=(10, 4))
plt.stem(rfftfreq_multi, magnitude_rfft, markerfmt=’.’, basefmt=’k-‘, linefmt=’b-‘)
plt.title(‘Magnitude Spectrum using rfft (Normalized)’)
plt.xlabel(‘Frequency (Hz)’)
plt.ylabel(‘Approx. Amplitude’)
plt.grid(True)
plt.show()

逆变换

recovered_signal_r = np.fft.irfft(rfft_multi)
print(f”Max error using rfft/irfft: {np.max(np.abs(signal_multi – recovered_signal_r)):.2e}”)
“`

rfft 的频谱图只显示了正频率部分（包括 0 Hz），但包含了所有必要的信息。对于处理实数信号，推荐使用 rfft 和 irfft 以提高效率。

采样率与奈奎斯特频率

  • 采样率 (fs)：每秒采集数据点的数量。它决定了我们能够分析的最高频率。
  • 奈奎斯特-香农采样定理：为了不失真地表示一个包含最高频率 f_max 的信号，采样率 fs 必须大于 2 * f_max (fs > 2 * f_max)。
  • 奈奎斯特频率 (f_nyquist)：等于采样率的一半 (fs / 2)。它是使用给定采样率能够可靠分析的最高频率。如果原始信号中包含高于奈奎斯特频率的成分，这些高频成分会在 FFT 结果中“折叠”或“混叠”到较低的频率上，导致频谱失真。
  • 实践：在进行 FFT 分析前，确保你的采样率足够高，能够覆盖你感兴趣的所有频率。通常选择 fs 至少是你关心的最高频率的 2.5 到 3 倍。在绘图时，通常将频率轴的范围限制在 [-fs/2, fs/2]（对于 fft）或 [0, fs/2]（对于 rfft）。

频谱泄漏与窗函数（简介）

FFT 假设输入信号是无限周期信号的一个完整周期。如果你的信号段不包含整数个周期（即信号的起始和结束点不连续），直接进行 FFT 会导致能量从真实的频率“泄漏”到邻近的频率仓中，在频谱图上表现为峰值变宽、旁边出现旁瓣。

窗函数 (Windowing) 是一种通过在信号应用 FFT 之前乘以一个特定形状的“窗”来减轻频谱泄漏的技术。窗函数在信号段的边缘逐渐减小到零，使得信号看起来更像是周期的。

常用的窗函数有 Hanning (np.hanning), Hamming (np.hamming), Blackman (np.blackman) 等。

“`python

简单演示加窗效果 (使用示例 5 信号)

window = np.hanning(N)
signal_windowed = signal_multi * window

fft_windowed = np.fft.fft(signal_windowed)
fft_windowed_shifted = np.fft.fftshift(fft_windowed)

幅度比较 (未归一化)

plt.figure(figsize=(12, 5))
plt.plot(freq_multi_shifted, np.abs(fft_multi_shifted), label=’Rectangular Window (No Window)’)
plt.plot(freq_multi_shifted, np.abs(fft_windowed_shifted), label=’Hanning Window’, alpha=0.7)
plt.title(‘Effect of Windowing on Magnitude Spectrum’)
plt.xlabel(‘Frequency (Hz)’)
plt.ylabel(‘Magnitude (Unnormalized)’)
plt.xlim(0, 40) # 放大观察峰值区域
plt.legend()
plt.grid(True)
plt.show()
“`

加窗会稍微降低峰值的高度并使其变宽，但能显著抑制旁瓣，提高频率分辨率（区分靠近的频率的能力）。选择哪种窗函数取决于具体的应用需求（例如，主瓣宽度和旁瓣抑制水平的权衡）。

幅度归一化

np.fft.fft 的直接输出幅度与输入信号的长度 N 成正比。为了得到与原始时域信号幅度直接对应的频谱幅度，需要进行归一化。

常用的归一化方法（针对 fft 的输出，并已 fftshift）：

  • 零频率 (DC) 分量：amplitude[0] = abs(fft_output[0]) / N
  • 正频率分量：amplitude[k] = 2 * abs(fft_output[k]) / N (对于 1 <= k < N//2)
  • 负频率分量：amplitude[k] = 2 * abs(fft_output[k]) / N (对于 N//2 < k < N)
  • 奈奎斯特频率（仅当 N 是偶数时存在，位于索引 N//2）：amplitude[N//2] = abs(fft_output[N//2]) / N

对于 rfft 的输出（索引 k 对应 rfftfreq）：

  • 零频率 (DC) 分量：amplitude[0] = abs(rfft_output[0]) / N
  • 其他频率分量：amplitude[k] = 2 * abs(rfft_output[k]) / N (对于 1 <= k < len(rfft_output)-1 或 1 <= k < len(rfft_output) 取决于 N 奇偶)
  • 奈奎斯特频率（仅当 N 是偶数时，位于最后一个索引）：amplitude[-1] = abs(rfft_output[-1]) / N

我们之前的示例代码中已经演示了这种归一化。

多维 FFT（fft2, ifft2, fftn, ifftn）

NumPy 的 FFT 功能不限于一维信号。

  • np.fft.fft2(a)：计算二维数组 a（例如，灰度图像）的二维 FFT。结果是一个复数二维数组，表示图像在水平和垂直方向上的空间频率成分。常用于图像滤波、压缩等。
  • np.fft.ifft2(a)：二维逆 FFT。
  • np.fft.fftn(a) 和 np.fft.ifftn(a)：推广到任意 N 维数组。

“`python

简单二维 FFT 示例

创建一个简单的二维模式 (例如，一个中心亮斑)

img_size = 64
X, Y = np.meshgrid(np.linspace(-1, 1, img_size), np.linspace(-1, 1, img_size))
radius = 0.2
image = (X2 + Y2 < radius**2).astype(float) # 圆形亮斑

fft_image = np.fft.fft2(image)
fft_image_shifted = np.fft.fftshift(fft_image) # 将零频移到中心

magnitude_spectrum_2d = np.log1p(np.abs(fft_image_shifted)) # 取对数使低频细节可见

plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.imshow(image, cmap=’gray’)
plt.title(‘Original Image’)
plt.axis(‘off’)

plt.subplot(1, 2, 2)
plt.imshow(magnitude_spectrum_2d, cmap=’gray’)
plt.title(‘2D FFT Magnitude Spectrum (Log Scale, Shifted)’)
plt.axis(‘off’)

plt.tight_layout()
plt.show()
“`

二维 FFT 的结果显示了图像的空间频率分布。低频成分（靠近中心）对应图像的平滑区域，高频成分（远离中心）对应图像的边缘和细节。

6. 总结与展望

本教程详细介绍了使用 NumPy fft 模块进行快速傅里叶变换的基础知识和实践方法。我们涵盖了：

  • FFT 的基本概念及其在信号分解中的作用。
  • NumPy 中用于一维 FFT 的核心函数：fft, ifft, fftfreq, fftshift。
  • 如何生成信号、计算 FFT、获取频率轴以及可视化频谱。
  • 通过实例演示了如何分析包含单一频率、多种频率以及噪声的信号。
  • 讨论了实数信号优化的 rfft/irfft，以及奈奎斯特频率、频谱泄漏、窗函数和幅度归一化等重要概念。
  • 简要介绍了多维 FFT (fft2)。

掌握 NumPy FFT 是进行数字信号处理、数据分析和科学计算的一项关键技能。通过本教程的学习，您应该能够开始使用 FFT 来探索数据中的频率信息。

下一步可以探索的方向包括：

  • 深入研究不同的窗函数及其特性。
  • 学习更高级的频谱分析技术，如功率谱密度 (PSD) 估计。
  • 将 FFT 应用于特定领域的问题，如音频处理（滤波、均衡器）、图像处理（去噪、特征提取）等。
  • 探索 SciPy 库 (scipy.signal)，它提供了更多高级的信号处理功能，包括滤波、谱估计、窗函数等，通常与 NumPy FFT 结合使用。

希望本指南为您打开了通往频域分析世界的大门！

发表评论

您的邮箱地址不会被公开。 必填项已用 * 标注

滚动至顶部