Python中的快速傅里叶变换（FFT）优化：在信号处理与序列建模中的应用

好的，我们开始今天的讲座，主题是Python中的快速傅里叶变换（FFT）优化及其在信号处理与序列建模中的应用。

引言：FFT的重要性

快速傅里叶变换 (FFT) 是一种高效的算法，用于计算离散傅里叶变换 (DFT)。DFT 将时域信号转换到频域，揭示信号的频率成分。FFT 的重要性在于它极大地降低了计算 DFT 的复杂度，从 O(N^2) 降低到 O(N log N)，其中 N 是信号的长度。这种效率的提升使得 FFT 在信号处理、图像处理、音频分析、通信系统以及各种科学和工程领域中得到广泛应用。在序列建模中，例如时间序列分析和自然语言处理，FFT 也常用于特征提取和模式识别。

DFT 与 FFT 的数学基础

首先，我们回顾一下 DFT 的定义。对于长度为 N 的离散信号 x[n]，其 DFT X[k] 定义为：

X[k] = ∑_{n=0}^{N-1} x[n] exp(-j 2 pi k * n / N) , k = 0, 1, …, N-1

其中：

• x[n] 是时域信号的第 n 个样本。
• X[k] 是频域信号的第 k 个频率分量。
• j 是虚数单位 (√-1)。
• exp() 是复指数函数。

直接计算 DFT 的复杂度为 O(N^2)，因为对于每个频率分量 X[k]，都需要进行 N 次乘法和加法运算。

FFT 通过巧妙地分解 DFT 计算来降低复杂度。最常用的 FFT 算法是 Cooley-Tukey 算法，它是一种分治算法，将 DFT 分解为更小的 DFT，递归地进行计算。 Cooley-Tukey 算法最有效的形式是将长度为 N 的 DFT 分解为两个长度为 N/2 的 DFT（假设 N 是偶数）。这种分解基于以下恒等式：

exp(-j 2 pi k n / N) = cos(-2 pi k n / N) + j sin(-2 pi k * n / N)

通过将 DFT 分解为偶数和奇数索引项，可以将计算量显著减少。 这种递归分解一直进行到达到基本情况，例如长度为 1 或 2 的 DFT。

Python 中的 FFT 实现：numpy.fft

Python 的 numpy 库提供了高效的 FFT 实现，位于 numpy.fft 模块中。该模块包含以下主要函数：

• fft(a, n=None, axis=-1, norm=None): 计算一维 DFT。
• ifft(a, n=None, axis=-1, norm=None): 计算一维 DFT 的逆变换 (IDFT)。
• fft2(a, s=None, axes=(-2, -1), norm=None): 计算二维 DFT。
• ifft2(a, s=None, axes=(-2, -1), norm=None): 计算二维 DFT 的逆变换 (IDFT)。
• fftn(a, s=None, axes=None, norm=None): 计算 N 维 DFT。
• ifftn(a, s=None, axes=None, norm=None): 计算 N 维 DFT 的逆变换 (IDFT)。
• rfft(a, n=None, axis=-1, norm=None): 计算实数信号的一维 DFT。由于实数信号的 DFT 具有共轭对称性，rfft 只返回正频率部分，从而节省计算量。
• irfft(a, n=None, axis=-1, norm=None): 计算实数信号的一维 DFT 的逆变换。
• hfft(a, n=None, axis=-1, norm=None): 计算 Hermitian 信号的一维 DFT。
• ihfft(a, n=None, axis=-1, norm=None): 计算 Hermitian 信号的一维 DFT 的逆变换。
• fftfreq(n, d=1.0): 生成 FFT 频率的数组。
• fftshift(x, axes=None): 将零频率分量移到频谱的中心。
• ifftshift(x, axes=None): 将零频率分量从频谱的中心移回原始位置。

代码示例：一维 FFT

import numpy as np
import matplotlib.pyplot as plt

# 生成一个示例信号
N = 512  # 信号长度
fs = 1000 # 采样率
t = np.arange(0, N/fs, 1/fs) # 时间向量
f1 = 50  # 信号频率 1
f2 = 120 # 信号频率 2
x = np.sin(2 * np.pi * f1 * t) + 0.5 * np.sin(2 * np.pi * f2 * t)  # 信号

# 计算 FFT
X = np.fft.fft(x)

# 计算频率轴
freq = np.fft.fftfreq(N, d=1/fs)

# 计算幅度谱
amplitude = np.abs(X)

# 绘制原始信号和幅度谱
plt.figure(figsize=(12, 6))

plt.subplot(2, 1, 1)
plt.plot(t, x)
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.title("Original Signal")

plt.subplot(2, 1, 2)
plt.plot(freq[:N//2], amplitude[:N//2]) # 只显示正频率部分
plt.xlabel("Frequency (Hz)")
plt.ylabel("Amplitude")
plt.title("Amplitude Spectrum")

plt.tight_layout()
plt.show()

这段代码首先生成一个包含两个正弦波的信号，然后使用 np.fft.fft() 计算其 FFT。 接着，使用 np.fft.fftfreq() 生成频率轴，并计算幅度谱。 最后，绘制原始信号和幅度谱。

实数 FFT (rfft) 的优势

如果输入信号是实数，可以使用 np.fft.rfft() 来提高效率。由于实数信号的 DFT 具有共轭对称性，只需要计算正频率部分。 rfft 返回的数组长度为 N//2 + 1，其中 N 是原始信号的长度。

import numpy as np
import matplotlib.pyplot as plt

# 生成一个示例实数信号
N = 512  # 信号长度
fs = 1000 # 采样率
t = np.arange(0, N/fs, 1/fs) # 时间向量
f1 = 50  # 信号频率
x = np.sin(2 * np.pi * f1 * t)  # 实数信号

# 计算 rfft
X = np.fft.rfft(x)

# 计算频率轴
freq = np.fft.fftfreq(N, d=1/fs)
freq = freq[:N//2 + 1]  # 调整频率轴的长度

# 计算幅度谱
amplitude = np.abs(X)

# 绘制幅度谱
plt.figure(figsize=(8, 4))
plt.plot(freq, amplitude)
plt.xlabel("Frequency (Hz)")
plt.ylabel("Amplitude")
plt.title("Amplitude Spectrum (rfft)")
plt.tight_layout()
plt.show()

FFT 优化策略

虽然 numpy.fft 已经经过高度优化，但在某些情况下，仍然可以采取一些策略来进一步提高 FFT 的性能：

1. 输入数据类型： numpy.fft 对 float64 和 complex128 数据类型的优化程度最高。 如果输入数据是其他类型，例如 float32 或 int16，先将其转换为 float64 或 complex128 可以提高性能。

2. 信号长度： FFT 的效率在信号长度为 2 的幂时最高。 如果信号长度不是 2 的幂，可以对其进行补零 (zero-padding) 到下一个 2 的幂。 补零不会改变信号的频率成分，但可以提高 FFT 的计算速度。

3. 避免不必要的拷贝： 在进行 FFT 之前，尽量避免对输入数据进行不必要的拷贝操作。 例如，可以使用 in-place 操作来修改数据，或者使用 view 来创建数据的视图，而不是拷贝。

4. 并行计算： 对于非常大的信号，可以使用并行计算来加速 FFT。 numpy.fft 本身不支持并行计算，但可以使用其他库，例如 dask.fft 或 cupy (如果使用 GPU)，来实现并行 FFT。

5. 选择合适的 FFT 算法： numpy.fft 默认使用 Cooley-Tukey 算法。 但在某些情况下，其他 FFT 算法可能更适合。 例如，如果信号长度是素数，可以使用 Bluestein 算法。 scipy.fft 模块提供了一些其他的 FFT 算法实现。

代码示例：补零 (Zero-Padding)

import numpy as np
import matplotlib.pyplot as plt

# 生成一个示例信号
N = 500  # 信号长度 (不是 2 的幂)
fs = 1000 # 采样率
t = np.arange(0, N/fs, 1/fs) # 时间向量
f1 = 50  # 信号频率
x = np.sin(2 * np.pi * f1 * t)  # 信号

# 补零到下一个 2 的幂
N_padded = 2**np.ceil(np.log2(N)).astype(int)
x_padded = np.pad(x, (0, N_padded - N), 'constant')

# 计算 FFT
X = np.fft.fft(x_padded)

# 计算频率轴
freq = np.fft.fftfreq(N_padded, d=1/fs)

# 计算幅度谱
amplitude = np.abs(X)

# 绘制幅度谱
plt.figure(figsize=(8, 4))
plt.plot(freq[:N_padded//2], amplitude[:N_padded//2])
plt.xlabel("Frequency (Hz)")
plt.ylabel("Amplitude")
plt.title("Amplitude Spectrum (Zero-Padded)")
plt.tight_layout()
plt.show()

代码示例：使用 scipy.fft

import numpy as np
from scipy.fft import fft, ifft, fftfreq
import matplotlib.pyplot as plt

# Generate a sample signal
N = 128
fs = 1000
t = np.arange(0, N/fs, 1/fs)
f1 = 50
x = np.sin(2 * np.pi * f1 * t)

# Compute FFT using scipy.fft
X = fft(x)

# Compute frequency axis
freq = fftfreq(N, 1/fs)

# Plot the amplitude spectrum
plt.figure(figsize=(8, 4))
plt.plot(freq[:N//2], np.abs(X[:N//2]))
plt.xlabel("Frequency (Hz)")
plt.ylabel("Amplitude")
plt.title("Amplitude Spectrum (scipy.fft)")
plt.tight_layout()
plt.show()

scipy.fft 提供了更多高级的 FFT 功能，例如不同的窗口函数和算法选择。

FFT 在信号处理中的应用

FFT 在信号处理中有广泛的应用，包括：

• 频谱分析： FFT 可以用于分析信号的频率成分，识别信号中的主要频率，以及检测信号中的噪声和干扰。
• 滤波： 可以通过在频域对信号进行滤波，然后使用 IDFT 将其转换回时域。 例如，可以使用 FFT 来实现低通滤波器、高通滤波器和带通滤波器。
• 卷积： 时域中的卷积等价于频域中的乘法。 可以使用 FFT 将两个信号转换到频域，然后将它们相乘，再使用 IDFT 将结果转换回时域，从而实现快速卷积。 这在信号处理和图像处理中非常有用。
• 相关性分析： 可以使用 FFT 来计算两个信号之间的互相关性。 互相关性可以用于检测信号中的相似性，以及确定信号之间的时间延迟。

代码示例：使用 FFT 实现卷积

import numpy as np
import matplotlib.pyplot as plt

# 生成两个示例信号
N = 512
x = np.random.randn(N)  # 信号 1
h = np.random.randn(N//4) # 信号 2 (滤波器)

# 使用 FFT 计算卷积
X = np.fft.fft(x, n=N) # 补零到相同长度
H = np.fft.fft(h, n=N)
Y = X * H
y = np.fft.ifft(Y).real # 取实部

# 使用 numpy.convolve 计算卷积 (作为参考)
y_conv = np.convolve(x, h, mode='same')

# 绘制结果
plt.figure(figsize=(12, 6))

plt.subplot(2, 1, 1)
plt.plot(y)
plt.xlabel("Sample")
plt.ylabel("Amplitude")
plt.title("Convolution using FFT")

plt.subplot(2, 1, 2)
plt.plot(y_conv)
plt.xlabel("Sample")
plt.ylabel("Amplitude")
plt.title("Convolution using numpy.convolve")

plt.tight_layout()
plt.show()

FFT 在序列建模中的应用

在序列建模中，FFT 可以用于：

• 特征提取： 将时间序列数据转换到频域，可以提取出一些有用的特征，例如频率、幅度和相位。 这些特征可以用于训练机器学习模型，例如分类器或回归器。
• 模式识别： 通过分析时间序列的频谱，可以识别出一些特定的模式。 例如，可以检测音频信号中的特定音调，或者识别心电图信号中的异常模式。
• 时间序列分解： 可以将时间序列分解为不同的频率成分，从而更好地理解时间序列的结构。
• 数据增强： 通过在频域对数据进行操作，再转换回时域，可以生成新的训练数据。

示例：使用 FFT 进行时间序列分类

假设我们有一组时间序列数据，每个时间序列属于不同的类别。我们可以使用 FFT 来提取时间序列的特征，然后使用分类器来预测时间序列的类别。

import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score

# 生成一些示例时间序列数据
def generate_time_series(n_samples, length, n_classes):
    X = np.zeros((n_samples, length))
    y = np.random.randint(0, n_classes, n_samples)
    for i in range(n_samples):
        f = np.random.uniform(1, 10)
        phase = np.random.uniform(0, 2*np.pi)
        amplitude = np.random.uniform(0.5, 1.5)
        X[i, :] = amplitude * np.sin(2 * np.pi * f * np.arange(length)/length + phase)
    return X, y

# 参数
n_samples = 100
length = 128
n_classes = 3

# 生成数据
X, y = generate_time_series(n_samples, length, n_classes)

# 提取 FFT 特征
X_fft = np.abs(np.fft.fft(X))[:, :length//2]  # 取正频率部分的幅度谱

# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X_fft, y, test_size=0.2, random_state=42)

# 训练随机森林分类器
model = RandomForestClassifier(n_estimators=100, random_state=42)
model.fit(X_train, y_train)

# 预测
y_pred = model.predict(X_test)

# 评估
accuracy = accuracy_score(y_test, y_pred)
print(f"Accuracy: {accuracy}")

表格：FFT 相关函数总结

函数	描述
fft	计算一维 DFT
ifft	计算一维 IDFT
fft2	计算二维 DFT
ifft2	计算二维 IDFT
fftn	计算 N 维 DFT
ifftn	计算 N 维 IDFT
rfft	计算实数信号的一维 DFT (返回正频率部分)
irfft	计算实数信号的一维 IDFT
fftfreq	生成 FFT 频率的数组
fftshift	将零频率分量移到频谱的中心
ifftshift	将零频率分量从频谱的中心移回原始位置
numpy.convolve	直接计算时域卷积 (可用于验证 FFT 卷积结果)
scipy.fft	提供了更多高级的 FFT 功能，例如不同的窗口函数和算法选择. 注意，函数名与 numpy.fft 类似，但属于不同的模块。

代码示例： 使用Dask 加速FFT

对于大规模数据，可以使用Dask进行并行FFT计算。

import numpy as np
import dask.array as da
import time

# 定义数组大小
array_size = 2**14  # 16384
chunk_size = 2**10   # 1024

# 创建一个大的NumPy数组
np_array = np.random.rand(array_size)

# 创建一个Dask数组
dask_array = da.from_array(np_array, chunks=(chunk_size,))

# 使用NumPy计算FFT并计时
start_time_np = time.time()
fft_np = np.fft.fft(np_array)
end_time_np = time.time()

# 使用Dask计算FFT并计时
start_time_da = time.time()
fft_da = da.fft.fft(dask_array).compute()
end_time_da = time.time()

# 输出时间
print(f"NumPy FFT Time: {end_time_np - start_time_np:.4f} seconds")
print(f"Dask FFT Time: {end_time_da - start_time_da:.4f} seconds")

# 验证结果（可选）
print(f"Max diff: {np.max(np.abs(fft_np - fft_da)):.4f}")

总结一下今天的学习内容

今天，我们深入探讨了 Python 中 FFT 的优化及其在信号处理和序列建模中的应用。我们学习了 DFT 和 FFT 的数学基础、numpy.fft 模块的使用方法、FFT 的优化策略，以及 FFT 在信号处理和序列建模中的各种应用。 通过学习，你掌握了使用 FFT 进行信号分析和特征提取的基本技能。

更多IT精英技术系列讲座，到智猿学院