使用Numpy rfft2手动重构函数时的失真问题

问题根源 #

1. 未考虑实FFT的共轭对称性：rfft2仅存储实输入FFT的非冗余正频率部分，负频率分量需通过共轭对称关系推导，但你的手动重构仅加入了正频率项，漏掉了对应的负频率共轭分量。
2. ky的遍历范围错误：rfft2对y轴（第一个维度）执行的是完整FFT（而非rfft），因此ky的有效范围是0到Ny-1，你的代码仅遍历到truncation-1，漏掉了高ky值的负频率分量。
3. 频谱泄漏的影响：你的函数cos(X)*cos(Y)中，X/Y的范围与FFT假设的周期不匹配，导致频谱泄漏，单一函数分量被分散到多个FFT系数中，仅保留少量系数必然导致失真，但正确加入共轭分量后，截断重构的效果会符合预期。

修正后的手动重构代码 #

import numpy as np
import matplotlib.pyplot as plt

# Parameters
a = 3.66
b = 1.75
Nx = 256  # Number of sample points in x
Ny = 256  # Number of sample points in y
truncation = 32  # Ensure this is less than the Nyquist frequency: Nx//2 + 1

# Create a grid of points in the domain
x = np.linspace(-a, a, Nx)
y = np.linspace(-b, b, Ny)
X, Y = np.meshgrid(x, y)

# Evaluate the function on the grid
F = np.cos(X) * np.cos(Y)

# Apply FFT
F_fft = np.fft.rfft2(F)

# Automatic reconstruction using irfft2 (for comparison)
F_auto_reconstructed = np.fft.irfft2(F_fft, s=(Nx, Ny))

# Manual Reconstruction
F_reconstructed = np.zeros_like(F, dtype=np.complex128)

# 1. Add positive x frequencies (kx from 0 to truncation-1)
for kx in range(truncation):
    # Positive y frequencies
    for ky in range(truncation):
        exp_term = np.exp(2j * np.pi * (kx * np.arange(Nx)[:, None] / Nx + ky * np.arange(Ny)[None, :] / Ny))
        F_reconstructed += F_fft[ky, kx] * exp_term
    
    # Negative y frequencies (ky from Ny - truncation +1 to Ny-1)
    for ky in range(Ny - truncation + 1, Ny):
        ky_pos = Ny - ky
        # 实函数FFT满足 F(kx, ky) = conj(F(kx, ky_pos))
        coeff = F_fft[ky_pos, kx].conj() if kx != 0 else F_fft[ky_pos, kx]
        exp_term = np.exp(2j * np.pi * (kx * np.arange(Nx)[:, None] / Nx + ky * np.arange(Ny)[None, :] / Ny))
        F_reconstructed += coeff * exp_term

# 2. Add negative x frequencies (kx from Nx - truncation +1 to Nx-1)
for kx in range(Nx - truncation + 1, Nx):
    kx_pos = Nx - kx
    if kx_pos >= truncation:
        continue  # 跳过超出截断范围的正频率对应项
    
    # Positive y frequencies
    for ky in range(truncation):
        coeff = F_fft[ky, kx_pos].conj()
        exp_term = np.exp(2j * np.pi * (kx * np.arange(Nx)[:, None] / Nx + ky * np.arange(Ny)[None, :] / Ny))
        F_reconstructed += coeff * exp_term
    
    # Negative y frequencies
    for ky in range(Ny - truncation + 1, Ny):
        ky_pos = Ny - ky
        coeff = F_fft[ky_pos, kx_pos]  # 双重共轭抵消，等于原系数
        exp_term = np.exp(2j * np.pi * (kx * np.arange(Nx)[:, None] / Nx + ky * np.arange(Ny)[None, :] / Ny))
        F_reconstructed += coeff * exp_term

# 3. 处理x方向的Nyquist频率（如果截断范围包含它）
if truncation > Nx // 2:
    kx = Nx // 2
    # Positive y frequencies
    for ky in range(truncation):
        exp_term = np.exp(2j * np.pi * (kx * np.arange(Nx)[:, None] / Nx + ky * np.arange(Ny)[None, :] / Ny))
        F_reconstructed += F_fft[ky, kx] * exp_term
    # Negative y frequencies
    for ky in range(Ny - truncation + 1, Ny):
        ky_pos = Ny - ky
        coeff = F_fft[ky_pos, kx].conj()
        exp_term = np.exp(2j * np.pi * (kx * np.arange(Nx)[:, None] / Nx + ky * np.arange(Ny)[None, :] / Ny))
        F_reconstructed += coeff * exp_term

# 归一化
F_reconstructed /= (Nx * Ny)

# Plotting for comparison
fig, axes = plt.subplots(1, 3, figsize=(18, 6))
axes[0].imshow(F, extent=[-a, a, -b, b], cmap='viridis')
axes[0].set_title('Original F')
axes[1].imshow(F_auto_reconstructed, extent=[-a, a, -b, b], cmap='viridis')
axes[1].set_title('Automatically Reconstructed F')
axes[2].imshow(np.real(F_reconstructed), extent=[-a, a, -b, b], cmap='viridis')
axes[2].set_title('Manually Reconstructed F')

plt.show()

关键说明 #

• 共轭对称处理：实函数的FFT满足F(kx, ky) = conjugate(F(Nx - kx, Ny - ky))，因此负频率分量的系数可通过正频率系数的共轭获取。
• 完整频率范围：手动重构时需覆盖所有符合截断条件的正负频率分量，而不仅仅是正频率的前N项。
• 频谱泄漏：由于你的函数周期与FFT假设的采样周期不匹配，会产生频谱泄漏，因此即使截断到较大的truncation值，重构结果也会与原函数有细微差异，但这是预期的截断效应，而非代码错误。

内容的提问来源于stack exchange，提问作者Theo Lavier