跳到主要内容

NumPy 傅里叶变换

傅里叶变换是信号处理、图像分析和科学计算的核心工具。NumPy 的 numpy.fft 模块提供了离散傅里叶变换（DFT）的高效实现，基于快速傅里叶变换（FFT）算法。本章将深入讲解傅里叶变换的原理和在 NumPy 中的应用。

傅里叶变换基础​

什么是傅里叶变换？​

傅里叶变换将信号从时域（时间空间）转换到频域（频率空间）。简单来说，任何周期函数都可以表示为不同频率的正弦波和余弦波的叠加。

对于离散信号，我们使用离散傅里叶变换（DFT）。其数学定义为：

Xk=∑n=0N−1xn⋅e−2πiknNX_k = \sum_{n=0}^{N-1} x_n \cdot e^{-2\pi i \frac{kn}{N}}

其中 xnx_n 是输入信号，XkX_k 是频域表示，NN 是信号长度。

为什么需要傅里叶变换？​

频谱分析：了解信号包含哪些频率成分，这在音频处理、振动分析中非常重要。

滤波：在频域中去除不需要的频率分量，如去除高频噪声。

卷积加速：时域卷积等于频域乘积，利用这一点可以大幅加速卷积运算。

图像处理：图像压缩、边缘检测、去噪等都离不开傅里叶变换。

时域与频域的关系​

import numpy as np

# 创建一个简单信号：两个正弦波的叠加
sampling_rate = 1000  # 采样率：每秒1000个样本
t = np.linspace(0, 1, sampling_rate, endpoint=False)  # 1秒的时间点

# 信号 = 50Hz正弦波 + 120Hz正弦波
signal = np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 120 * t)

print(f"信号长度: {len(signal)}")
print(f"采样间隔: {t[1] - t[0]:.6f} 秒")
print(f"时间范围: {t[0]} 到 {t[-1]} 秒")

一维傅里叶变换​

fft 函数​

np.fft.fft 计算一维离散傅里叶变换：

import numpy as np

# 创建信号
sampling_rate = 1000
t = np.linspace(0, 1, sampling_rate, endpoint=False)
signal = np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 120 * t)

# 计算FFT
fft_result = np.fft.fft(signal)

print(f"FFT结果形状: {fft_result.shape}")
print(f"FFT结果类型: {fft_result.dtype}")  # 复数类型

# FFT结果是复数，包含幅度和相位信息
print(f"前5个FFT值: {fft_result[:5]}")

频率轴​

使用 np.fft.fftfreq 获取对应的频率值：

import numpy as np

sampling_rate = 1000
n = 1000  # 采样点数

# 获取频率轴
frequencies = np.fft.fftfreq(n, d=1/sampling_rate)

print(f"频率范围: {frequencies.min():.1f} Hz 到 {frequencies.max():.1f} Hz")
print(f"频率间隔: {frequencies[1] - frequencies[0]:.4f} Hz")

# 频率轴的特点：
# [0, 1, 2, ..., n/2-1, -n/2, ..., -2, -1] * (采样率/n)
print(f"\n前10个频率: {frequencies[:10]}")
print(f"后10个频率: {frequencies[-10:]}")

幅度谱和相位谱​

FFT 结果是复数，可以分解为幅度和相位：

import numpy as np

sampling_rate = 1000
t = np.linspace(0, 1, sampling_rate, endpoint=False)
signal = np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 120 * t)

fft_result = np.fft.fft(signal)
frequencies = np.fft.fftfreq(len(signal), d=1/sampling_rate)

# 幅度谱（取绝对值）
amplitude = np.abs(fft_result)

# 功率谱（幅度的平方）
power = amplitude ** 2

# 相位谱
phase = np.angle(fft_result)

# 只看正频率部分
positive_mask = frequencies >= 0
positive_freqs = frequencies[positive_mask]
positive_amplitude = amplitude[positive_mask]

print("幅度谱峰值分析:")
# 找到幅度最大的几个频率
top_indices = np.argsort(positive_amplitude)[-5:][::-1]
for idx in top_indices:
    if positive_freqs[idx] > 0:
        print(f"频率 {positive_freqs[idx]:.1f} Hz: 幅度 {positive_amplitude[idx]:.2f}")

ifft 函数​

np.fft.ifft 是逆傅里叶变换，将频域信号转回时域：

import numpy as np

# 原始信号
t = np.linspace(0, 1, 1000, endpoint=False)
original = np.sin(2 * np.pi * 50 * t)

# FFT
fft_result = np.fft.fft(original)

# IFFT
reconstructed = np.fft.ifft(fft_result)

# 验证重建
print(f"原始信号前5个值: {original[:5]}")
print(f"重建信号前5个值: {reconstructed.real[:5]}")  # 取实部
print(f"重建误差: {np.max(np.abs(original - reconstructed.real)):.2e}")

# 注意：IFFT结果可能是复数（由于数值误差），通常只取实部

fftshift 和 ifftshift​

FFT 结果的默认顺序是：零频率在开头，然后是正频率，最后是负频率。fftshift 将零频率移到中心：

import numpy as np

sampling_rate = 1000
n = 1000

frequencies = np.fft.fftfreq(n, d=1/sampling_rate)
print("原始频率顺序:")
print(f"前5个: {frequencies[:5]}")
print(f"后5个: {frequencies[-5:]}")

# 使用 fftshift 重新排列
shifted_freqs = np.fft.fftshift(frequencies)
print("\nshift后频率顺序:")
print(f"前5个: {shifted_freqs[:5]}")
print(f"后5个: {shifted_freqs[-5:]}")
print(f"中心频率: {shifted_freqs[n//2]}")  # 0 Hz

可视化示例​

import numpy as np

sampling_rate = 1000
t = np.linspace(0, 1, sampling_rate, endpoint=False)
signal = np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 120 * t)

# 计算FFT并shift
fft_result = np.fft.fft(signal)
fft_shifted = np.fft.fftshift(fft_result)
frequencies = np.fft.fftshift(np.fft.fftfreq(len(signal), d=1/sampling_rate))

# 幅度谱（归一化）
amplitude = np.abs(fft_shifted) / len(signal)

print("幅度谱（中心化）:")
# 找到主要频率成分
threshold = 0.1
significant = amplitude > threshold
for freq, amp in zip(frequencies[significant], amplitude[significant]):
    if abs(freq) > 0:
        print(f"频率 {freq:6.1f} Hz: 幅度 {amp:.3f}")

实数信号的 FFT​

对于实数信号，FFT 结果具有对称性：负频率分量是正频率分量的共轭。因此可以使用更高效的 rfft：

rfft 函数​

rfft 只计算正频率部分，节省计算和存储：

import numpy as np

sampling_rate = 1000
t = np.linspace(0, 1, sampling_rate, endpoint=False)
signal = np.sin(2 * np.pi * 50 * t)  # 实数信号

# 标准FFT
fft_result = np.fft.fft(signal)
print(f"fft 结果长度: {len(fft_result)}")  # 1000

# 实数FFT（更高效）
rfft_result = np.fft.rfft(signal)
print(f"rfft 结果长度: {len(rfft_result)}")  # 501

# 频率轴
rfft_freqs = np.fft.rfftfreq(len(signal), d=1/sampling_rate)
print(f"rfft 频率范围: {rfft_freqs[0]:.1f} 到 {rfft_freqs[-1]:.1f} Hz")

# 对比计算结果
print(f"\n幅度对比:")
print(f"fft 第一个值: {np.abs(fft_result[0]):.2f}")
print(f"rfft 第一个值: {np.abs(rfft_result[0]):.2f}")

irfft 函数​

irfft 是 rfft 的逆变换：

import numpy as np

# 原始实数信号
t = np.linspace(0, 1, 1000, endpoint=False)
original = np.sin(2 * np.pi * 50 * t)

# rfft -> irfft
rfft_result = np.fft.rfft(original)
reconstructed = np.fft.irfft(rfft_result, n=len(original))

print(f"原始长度: {len(original)}")
print(f"重建长度: {len(reconstructed)}")
print(f"重建误差: {np.max(np.abs(original - reconstructed)):.2e}")

二维 FFT​

二维 FFT 常用于图像处理：

fft2 函数​

import numpy as np

# 创建一个简单的二维图案
x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
X, Y = np.meshgrid(x, y)

# 正弦波图案
image = np.sin(2 * np.pi * X / 3) + np.sin(2 * np.pi * Y / 5)

print(f"图像形状: {image.shape}")

# 二维FFT
fft2d = np.fft.fft2(image)
fft2d_shifted = np.fft.fftshift(fft2d)

# 幅度谱
amplitude_2d = np.abs(fft2d_shifted)

print(f"FFT结果形状: {fft2d.shape}")
print(f"幅度谱最大值: {amplitude_2d.max():.2f}")

# 逆变换
reconstructed = np.fft.ifft2(fft2d)
print(f"重建误差: {np.max(np.abs(image - reconstructed.real)):.2e}")

图像滤波示例​

import numpy as np

# 创建带噪声的图像
np.random.seed(42)
x = np.linspace(-5, 5, 128)
y = np.linspace(-5, 5, 128)
X, Y = np.meshgrid(x, y)
clean_image = np.sin(2 * np.pi * X / 3) + np.sin(2 * np.pi * Y / 5)

# 添加高频噪声
noisy_image = clean_image + 0.5 * np.random.randn(*clean_image.shape)

# 频域滤波：去除高频噪声
fft_noisy = np.fft.fft2(noisy_image)
fft_shifted = np.fft.fftshift(fft_noisy)

# 创建低通滤波器
rows, cols = noisy_image.shape
crow, ccol = rows // 2, cols // 2
radius = 30  # 截止频率

Y_grid, X_grid = np.ogrid[:rows, :cols]
mask = (X_grid - ccol)**2 + (Y_grid - crow)**2 <= radius**2

# 应用滤波器
fft_filtered = fft_shifted * mask

# 逆变换
filtered_image = np.fft.ifft2(np.fft.ifftshift(fft_filtered)).real

print("图像滤波结果:")
print(f"原始图像噪声方差: {np.var(noisy_image - clean_image):.4f}")
print(f"滤波后噪声方差: {np.var(filtered_image - clean_image):.4f}")

FFT 的实际应用​

应用1：频谱分析​

分析信号的频率成分：

import numpy as np

def analyze_spectrum(signal, sampling_rate):
    """分析信号的频谱"""
    n = len(signal)
    
    # 计算FFT
    fft_result = np.fft.rfft(signal)
    frequencies = np.fft.rfftfreq(n, d=1/sampling_rate)
    
    # 计算幅度谱（归一化）
    amplitude = np.abs(fft_result) / n * 2
    
    return frequencies, amplitude

# 示例：分析音乐音符
sampling_rate = 44100  # CD音质
duration = 1.0  # 1秒
t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)

# A4音符 (440 Hz) + E5音符 (659 Hz)
note_signal = np.sin(2 * np.pi * 440 * t) + 0.7 * np.sin(2 * np.pi * 659 * t)

freqs, amps = analyze_spectrum(note_signal, sampling_rate)

# 找到主要频率
threshold = 0.1
significant = amps > threshold
print("检测到的频率成分:")
for f, a in zip(freqs[significant], amps[significant]):
    if f > 0 and f < 2000:  # 只显示 2kHz 以下的
        print(f"  {f:.1f} Hz: 幅度 {a:.3f}")

应用2：快速卷积​

利用 FFT 加速大数组的卷积运算：

import numpy as np
import time

def convolve_fft(a, b):
    """使用FFT进行卷积"""
    # 计算输出长度
    n = len(a) + len(b) - 1
    
    # 找到合适的FFT长度（2的幂次效率最高）
    n_fft = 2 ** int(np.ceil(np.log2(n)))
    
    # FFT
    fft_a = np.fft.fft(a, n_fft)
    fft_b = np.fft.fft(b, n_fft)
    
    # 频域相乘
    fft_result = fft_a * fft_b
    
    # IFFT
    result = np.fft.ifft(fft_result)[:n].real
    
    return result

# 性能测试
np.random.seed(42)
a = np.random.rand(10000)
b = np.random.rand(1000)

# 使用FFT卷积
start = time.time()
result_fft = convolve_fft(a, b)
time_fft = time.time() - start

# 使用np.convolve
start = time.time()
result_np = np.convolve(a, b, mode='full')
time_np = time.time() - start

print(f"FFT卷积时间: {time_fft:.4f}s")
print(f"np.convolve时间: {time_np:.4f}s")
print(f"结果误差: {np.max(np.abs(result_fft - result_np)):.2e}")

应用3：信号去噪​

import numpy as np

def denoise_signal(signal, threshold_ratio=0.1):
    """使用频域滤波去噪"""
    # FFT
    fft_result = np.fft.fft(signal)
    
    # 计算阈值
    amplitude = np.abs(fft_result)
    threshold = threshold_ratio * np.max(amplitude)
    
    # 软阈值：去除小的频率分量
    mask = amplitude > threshold
    fft_denoised = fft_result * mask
    
    # 逆变换
    denoised = np.fft.ifft(fft_denoised).real
    
    return denoised

# 创建带噪声的信号
np.random.seed(42)
t = np.linspace(0, 2, 1000, endpoint=False)
clean = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 10 * t)
noisy = clean + 0.3 * np.random.randn(len(t))

# 去噪
denoised = denoise_signal(noisy, threshold_ratio=0.05)

print("去噪效果:")
print(f"原始噪声方差: {np.var(noisy - clean):.4f}")
print(f"去噪后方差: {np.var(denoised - clean):.4f}")

应用4：相关性计算​

FFT 可用于快速计算相关函数：

import numpy as np

def correlate_fft(a, b):
    """使用FFT计算互相关"""
    n = len(a) + len(b) - 1
    n_fft = 2 ** int(np.ceil(np.log2(n)))
    
    # FFT（注意b需要共轭）
    fft_a = np.fft.fft(a, n_fft)
    fft_b = np.fft.fft(b, n_fft)
    
    # 频域相乘（相关是共轭相乘）
    fft_result = fft_a * np.conj(fft_b)
    
    # IFFT
    result = np.fft.ifft(fft_result)[:n].real
    
    return result

# 示例：信号检测
t = np.linspace(0, 1, 1000, endpoint=False)
template = np.sin(2 * np.pi * 10 * t[:100])  # 模板信号

# 在信号中嵌入模板
signal = np.zeros(1000)
signal[500:600] = template
signal += 0.1 * np.random.randn(1000)  # 添加噪声

# 计算相关
correlation = correlate_fft(signal, template)

# 找到最大相关位置
max_pos = np.argmax(correlation)
print(f"模板检测位置: {max_pos} (实际嵌入位置: 500)")

FFT 参数详解​

n 参数：输出长度​

import numpy as np

signal = np.array([1, 2, 3, 4, 5])

# 默认：输出长度等于输入长度
fft_default = np.fft.fft(signal)
print(f"默认长度: {len(fft_default)}")

# 指定更长的输出（零填充）
fft_padded = np.fft.fft(signal, n=10)
print(f"填充后长度: {len(fft_padded)}")

# 零填充的效果：频谱更密集
freqs_default = np.fft.fftfreq(len(signal))
freqs_padded = np.fft.fftfreq(10)

print(f"\n默认频率间隔: {freqs_default[1] - freqs_default[0]:.4f}")
print(f"填充后频率间隔: {freqs_padded[1] - freqs_padded[0]:.4f}")

norm 参数：归一化方式​

import numpy as np

signal = np.array([1, 2, 3, 4])

# 默认（backward）：正向变换不归一化，逆向变换除以n
fft_backward = np.fft.fft(signal, norm='backward')
ifft_backward = np.fft.ifft(fft_backward, norm='backward')

# 正交（ortho）：正逆变换都除以 sqrt(n)
fft_ortho = np.fft.fft(signal, norm='ortho')
ifft_ortho = np.fft.ifft(fft_ortho, norm='ortho')

# 前向（forward）：正向变换除以n，逆向变换不归一化
fft_forward = np.fft.fft(signal, norm='forward')

print("不同归一化方式:")
print(f"backward FFT[0]: {fft_backward[0]:.2f}")
print(f"ortho FFT[0]: {fft_ortho[0]:.2f}")
print(f"forward FFT[0]: {fft_forward[0]:.2f}")

# 验证逆变换
print(f"\nbackward 逆变换验证: {np.allclose(ifft_backward, signal)}")
print(f"ortho 逆变换验证: {np.allclose(ifft_ortho, signal)}")

常见问题与技巧​

频率分辨率​

频率分辨率由信号长度和采样率决定：

import numpy as np

sampling_rate = 1000  # Hz
duration = 1.0        # 秒
n = int(sampling_rate * duration)

# 频率分辨率
freq_resolution = sampling_rate / n
print(f"频率分辨率: {freq_resolution} Hz")

# 要分辨更接近的频率，需要更长的信号
duration_long = 2.0
n_long = int(sampling_rate * duration_long)
freq_resolution_long = sampling_rate / n_long
print(f"更长时间分辨率: {freq_resolution_long} Hz")

# 或者增加采样率（但不能提高分辨率，只有更长的观测时间才能）

频谱泄漏​

当信号不是完整周期时会发生频谱泄漏：

import numpy as np

sampling_rate = 1000
t = np.linspace(0, 1, sampling_rate, endpoint=False)

# 整周期信号（50 Hz 在 1 秒内有 50 个完整周期）
signal_exact = np.sin(2 * np.pi * 50 * t)

# 非整周期信号（50.5 Hz）
signal_leak = np.sin(2 * np.pi * 50.5 * t)

# 对比频谱
fft_exact = np.abs(np.fft.rfft(signal_exact))
fft_leak = np.abs(np.fft.rfft(signal_leak))

freqs = np.fft.rfftfreq(sampling_rate, d=1/sampling_rate)

print("频谱泄漏对比:")
exact_peak_idx = np.argmax(fft_exact)
leak_peak_idx = np.argmax(fft_leak)
print(f"整周期峰值频率: {freqs[exact_peak_idx]:.1f} Hz")
print(f"非整周期峰值频率: {freqs[leak_peak_idx]:.1f} Hz")

# 使用窗函数可以减少泄漏
window = np.hanning(sampling_rate)
signal_windowed = signal_leak * window
fft_windowed = np.abs(np.fft.rfft(signal_windowed))
print(f"加窗后峰值频率: {freqs[np.argmax(fft_windowed)]:.1f} Hz")

窗函数​

import numpy as np

# 常用窗函数
n = 100

# 矩形窗（无窗）
rectangular = np.ones(n)

# 汉宁窗
hanning = np.hanning(n)

# 汉明窗
hamming = np.hamming(n)

# 布莱克曼窗
blackman = np.blackman(n)

print("窗函数特性:")
print(f"汉宁窗最大值: {hanning.max():.4f}")
print(f"汉明窗最大值: {hamming.max():.4f}")
print(f"布莱克曼窗最大值: {blackman.max():.4f}")

FFT 函数参考​

函数说明
fft(a)一维FFT
ifft(a)一维逆FFT
fft2(a)二维FFT
ifft2(a)二维逆FFT
fftn(a)N维FFT
ifftn(a)N维逆FFT
rfft(a)实数FFT（只返回正频率）
irfft(a)实数逆FFT
fftfreq(n)返回FFT频率轴
rfftfreq(n)返回rfft频率轴
fftshift(a)将零频率移到中心
ifftshift(a)fftshift的逆操作

小结​

本章介绍了 NumPy 中傅里叶变换的核心功能：

  1. 基础概念：理解时域与频域的转换关系
  2. 一维FFT：fft、ifft、fftfreq 的使用
  3. 实数FFT：rfft 和 irfft 的高效用法
  4. 二维FFT：图像处理中的 fft2 和 ifft2
  5. 实际应用：频谱分析、快速卷积、信号去噪
  6. 技巧：频率分辨率、频谱泄漏、窗函数

傅里叶变换是信号处理和科学计算的基础工具，理解其原理和正确使用方法对于数据分析至关重要。

练习​

  1. 分析一段音频信号，找出其中的主频率成分
  2. 使用 FFT 实现图像的低通滤波和高通滤波
  3. 比较不同窗函数对频谱泄漏的影响
  4. 使用 FFT 实现两个大数组的快速卷积
  5. 实现一个简单的频谱分析仪，显示信号的时域波形和频域频谱