手动实现逆FFT时误差沿x轴递增的原因及修正方法咨询

问题描述 #

通过NumPy的FFT组件手动计算逆FFT，将重建信号与原始正弦信号对比时，误差（original - rebuilt）沿x轴逐渐增大；但使用内置函数ifft重建信号时无此现象，请问原因及修正方法？

用户提供的代码如下：

import numpy as np
from scipy import *
from matplotlib import pyplot as plt
tempdata = np.loadtxt("sinx_80.dat")


x_positions = tempdata[:, 0]
delta_x=x_positions[1] - x_positions[0]
roughness_data = tempdata[:, 1]


fft_1d = np.fft.fft(roughness_data)
num_x = len(fft_1d)

freq_x = np.fft.fftfreq(num_x, d=delta_x)
reconstructed_surface=np.zeros(num_x)
x_grid = (x_positions - x_positions[0])

for n in range(num_x):
    amplitude = np.abs(fft_1d[n]) / (num_x)
    phase = np.angle(fft_1d[n])
    reconstructed_surface += amplitude * np.cos(2 * np.pi * (freq_x[n] * x_grid) + phase)

difference = roughness_data - reconstructed_surface


max_error = np.max(np.abs(difference))
rms_error = np.sqrt(np.mean(difference**2))

print(f"Max Error: {max_error}")
print(f"RMS Error: {rms_error}")


plt.figure(figsize=(12, 8))
plt.subplot(3, 1, 1)
plt.plot(x_positions, roughness_data, label="Original Data", color="blue")
plt.title("Original Data")
plt.legend()
plt.grid()


plt.subplot(3, 1, 2)
plt.plot(x_positions, reconstructed_surface, label="Reconstructed Data", color="orange")
plt.title("Reconstructed Data")
plt.legend()
plt.grid()


plt.subplot(3, 1, 3)
plt.plot(x_positions, difference, label="Difference (Original - Reconstructed)", color="red")
plt.title("Difference")
plt.legend()
plt.grid()


plt.show()

误差表现：误差沿x轴逐渐递增（对应图示趋势）

原因分析 #

• 频率相位符号错误：np.fft.fftfreq返回的频率包含正、负两类，逆FFT的余弦合成中，负频率项的相位符号需要反向。你当前对所有频率项统一使用+ phase，导致正负频率的相位无法正确抵消虚部，产生累积的相位偏移，最终表现为误差随x递增。
• 振幅处理逻辑缺失：FFT中正负频率是共轭对，除直流分量（n=0）和Nyquist频率（n=num_x//2）外，每个正频率的振幅需要乘以2才能和对应的负频率共同还原原始信号的幅度，你当前的振幅计算没有做这个处理，会导致信号幅度偏差，但这不是误差递增的核心原因。

修正方法 #

核心修正点 #

1. 区分正、负频率，负频率项使用- phase参与余弦合成；
2. 对正频率（除直流和Nyquist频率）的振幅乘以2，合并共轭频率对的贡献。

修改后的代码 #

import numpy as np
from matplotlib import pyplot as plt  # 移除未使用的scipy导入

tempdata = np.loadtxt("sinx_80.dat")

x_positions = tempdata[:, 0]
delta_x = x_positions[1] - x_positions[0]
roughness_data = tempdata[:, 1]

fft_1d = np.fft.fft(roughness_data)
num_x = len(fft_1d)

freq_x = np.fft.fftfreq(num_x, d=delta_x)
reconstructed_surface = np.zeros(num_x)
x_grid = x_positions - x_positions[0]

for n in range(num_x):
    amplitude = np.abs(fft_1d[n]) / num_x
    phase = np.angle(fft_1d[n])
    # 区分正负频率调整相位符号
    if freq_x[n] < 0:
        reconstructed_surface += amplitude * np.cos(2 * np.pi * freq_x[n] * x_grid - phase)
    else:
        # 正频率中，除直流和Nyquist频率外，振幅乘以2
        if n != 0 and n != num_x // 2:
            amplitude *= 2
        reconstructed_surface += amplitude * np.cos(2 * np.pi * freq_x[n] * x_grid + phase)

difference = roughness_data - reconstructed_surface

max_error = np.max(np.abs(difference))
rms_error = np.sqrt(np.mean(difference**2))

print(f"Max Error: {max_error}")
print(f"RMS Error: {rms_error}")

plt.figure(figsize=(12, 8))
plt.subplot(3, 1, 1)
plt.plot(x_positions, roughness_data, label="Original Data", color="blue")
plt.title("Original Data")
plt.legend()
plt.grid()

plt.subplot(3, 1, 2)
plt.plot(x_positions, reconstructed_surface, label="Reconstructed Data", color="orange")
plt.title("Reconstructed Data")
plt.legend()
plt.grid()

plt.subplot(3, 1, 3)
plt.plot(x_positions, difference, label="Difference (Original - Reconstructed)", color="red")
plt.title("Difference")
plt.legend()
plt.grid()

plt.show()

原理说明 #

FFT结果中，正频率fft_1d[n]和负频率fft_1d[num_x-n]是共轭对称的，手动合成时需要将这一对的贡献合并：A*cos(ωt+φ) + A*cos(-ωt+φ) = 2A*cos(ωt)cos(φ)，才能还原原始实信号。如果不对负频率相位取反，会变成A*cos(ωt+φ) + A*cos(-ωt-φ) = 2A*cos(ωt+φ)，引入额外的相位偏移，随x增大累积形成递增误差。

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