使用Python中的FFT重构原始信号

Question

使用Python中的FFT重构原始信号

4

我有一个通过分析生成的频谱图，其中 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()

- Péter Leéh

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Cris Luengo · Accepted Answer

我认为fftfreq并不是你想象中的那样。对于fft(ifft(y)的xf与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在0.0004777步长下从3.0到1.0变化。你需要扩充你的数据，使其值从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频率轴)。