python - 在python中用FFT重建原始信号

Question

我有一个分析生成的频谱，其中 x 轴表示角频率，y 表示强度。频谱以某个频率值为中心，这通常称为信号的中心频率(图片上的蓝色图表)。
我想对数据执行 IFFT 到时域，用高斯曲线切割它的有用部分，然后 FFT 回到原始域。我的问题是，在 IFFT(FFT(signal)) 之后，中心频率丢失了:我按形状取回频谱，但它始终以 0 为中心(橙色图)。

目前我对此的解决方案非常糟糕:我缓存了原始的 x 轴，并在 FFT 调用时恢复它。这显然有很多缺点，我想改进它。下面我包含了一个演示问题的小演示。我的问题是:这可以以更优雅的方式解决吗？有没有办法在这个过程中不丢失中心频率？

import numpy as np
from scipy.fftpack import fft, ifft, fftshift, fftfreq
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt

C_LIGHT = 299.793

def generate_data(start, stop, center, delay, GD=0, resolution=0.1):
    window = 8 * np.log(2) / 50
    lamend = 2 * np.pi * C_LIGHT / start
    lamstart = 2 * np.pi * C_LIGHT/stop
    lam = np.arange(lamstart, lamend + resolution, resolution) 
    omega = 2 * np.pi * C_LIGHT / lam 
    relom = omega - center
    _i = np.exp(-(relom) ** 2 / window)
    i = 2 * _i + 2 * np.cos(relom * GD + (omega * delay)) * np.sqrt(_i * _i)
    return omega, i


if __name__ == '__main__':

    # Generate data
    x, y = generate_data(1, 3, 2, 800, GD=0)

    # Linearly interpolate to be evenly spaced
    xs = np.linspace(x[0], x[-1], len(x))
    intp = interp1d(x, y, kind='linear')
    ys = intp(xs)
    x, y = xs, ys
    plt.plot(x, y, label='original')

    # IFFT 
    xt = fftfreq(len(x), d=(x[0]-x[1])/(2*np.pi))
    yt = ifft(y)
    # plt.plot(xt, np.abs(yt))

    # FFT back
    xf = fftshift(fftfreq(len(xt), d=(xt[0]-xt[1])/(2*np.pi)))
    yf = fft(yt)
    plt.plot(xf, np.abs(yf), label='after transforms')
    plt.legend()
    plt.grid()
    plt.show()

Answer 1

我认为fftfreq不会做你认为它会做的事情。 xf为 fft(ifft(y)与 x 相同，您不应该尝试重新计算它。转到另一个域然后再返回时，x 轴不会改变。

另外，请注意 fftfreq返回给定长度和给定样本间距的信号的离散傅立叶变换的频域坐标。它不会做相反的事情，在应用逆离散傅立叶变换后，您不能使用它来确定空间域中的坐标。 (它返回的间距有效，但坐标集无效。)

    plt.plot(x, y, label='original')

    # IFFT 
    yt = ifft(y)
    # plt.plot(np.abs(yt))

    # FFT back
    yf = fft(yt)
    plt.plot(x, np.real(yf), label='after transforms')
    plt.legend()
    plt.grid()
    plt.show()

您的代码的另一个问题是 ifft(y)假定沿 x 轴的一组固定值。您的 x不符合这个。因此，您获得的空间域信号没有意义。

运行你的代码，我看到 x从 3.0 到 1.0，步长为 0.0004777。您必须增加数据，使值从 0.0 运行到 6.0，区域 (3.0, 6.0) 是区域 (0.0, 3.0) 的共轭对称副本。该区域对应于负频率，根据频域的周期性(F[n]==F[n+N]，其中 N 为样本数)。用零填充区域 (0.0, 1.0)。

鉴于频域中的标准化 x 轴，xf = fftfreq(len(xt), d=(xt[1]-xt[0]))应该重建x轴。但是你需要计算xt适本地:xt = np.linspace(0, 1/(x[1]-x[0]), len(x), endpoint=False) (使用 x 标准化 DFT 频率轴)。