离散傅里叶变换：如何正确使用 fftshift 和 fft

Question

我想对 numpy 数组 Y 进行数值计算。为了进行测试，我使用高斯函数 Y = exp(-x^2)。 （符号）傅立叶变换为 Y' = 常数 * exp(-k^2/4)。

import numpy
X = numpy.arange(-100,100)
Y = numpy.exp(-(X/5.0)**2)

幼稚的方法失败了：

from numpy.fft import *
from matplotlib import pyplot

def plotReIm(x,y):
    f = pyplot.figure()
    ax = f.add_subplot(111)
    ax.plot(x, numpy.real(y), 'b', label='R()')
    ax.plot(x, numpy.imag(y), 'r:', label='I()')
    ax.plot(x, numpy.abs(y), 'k--', label='abs()')
    ax.legend()


Y_k = fftshift(fft(Y))
k = fftshift(fftfreq(len(Y)))
plotReIm(k,Y_k)

real(Y_k)在正值和负值之间跳跃，这对应于跳跃阶段，而符号结果中不存在该跳跃阶段。这当然是不可取的。 （结果在技术上是正确的，abs(Y_k) 给出了预期的幅值 ifft(Y_k) 是 Y。）

这里，函数 fftshift() 渲染数组 k 单调递增并相应地改变 Y_k。对两个向量应用此操作不会改变 zip(k, Y_k) 对。

此更改似乎解决了问题：

Y_k = fftshift(fft(ifftshift(Y)))
k = fftshift(fftfreq(len(Y)))
plotReIm(k,Y_k)

如果需要单调 Y 和 Y_k，这是使用 fft() 函数的正确方法吗？

上述操作的反向操作是：

Yx = fftshift(ifft(ifftshift(Y_k)))
x = fftshift(fftfreq(len(Y_k), k[1] - k[0]))
plotReIm(x,Yx) 

对于这种情况，文档明确指出，Y_k 的排序必须与 fft() 和 fftfreq() 的输出兼容，我们可以通过应用 ifftshift() 来实现。

这些问题已经困扰我很长时间了： fft() 和 ifft() 的输出和输入数组是否始终使得 a[0] 应该包含零频率项 a[1:n/2 +1] 应包含正频率项，a[n/2+1:] 应包含负频率项，按负频率递减的顺序 [numpy 参考]，其中“频率”是自变量？

答案 高斯的傅立叶变换不是高斯并没有回答我的问题。

Answer 1

FFT 可以被认为是产生一组向量，每个向量都有幅度和相位。 fft_shift 操作将零相位角的参考点从 FFT 孔径的边缘更改为原始输入数据矢量的中心。

完成此操作后，结果的相位（以及复矢量的实数分量）有时不会那么“跳跃”，特别是如果某些输入函数被加窗，使其在 FFT 孔径边缘周围不连续。或者，如果输入围绕 FFT 孔径中心对称，则 FFT 结果的相位在 fft_shift 后将始终为零。

fft_shift 可以通过 N/2 的向量旋转来完成，或者通过简单地翻转 FFT 结果中的交替符号位来完成，这可能对 CPU dcache 更友好。

Answer 2

fft（和 ifft）输出的定义如下：http://docs.scipy.org/doc/numpy/reference/routines.fft.html#background-information

这是例程计算的内容，没有了，没有了 较少的。观察到离散傅里叶变换与连续傅里叶变换有很大不同。对于密集采样函数，两者之间存在关系，但除了 fftshift 之外，该关系还涉及相位因子和缩放。这就是您在图中看到的振荡的原因。您可以根据上述 DFT 数学公式自行计算出必要的相位因子。

Answer 3

考虑你的代码，但让我们缩小你的范围

import numpy as np
X = np.array([-1, 0, 1])
Y = np.exp(-(X/5.0)**2)

来实现

ifftshift([-1,0,1]) #[0,1,-1]

你希望你的离散零索引 Y[0] 与 X==0 重合，这可以通过注意 ffthift 和 ifftshift 只是循环移位操作

import numpy as np
from np.fft import fftshift
from np.fft import ifftshift
def _fftshift(x):
    centerElement = len(x)//2
    return np.roll(x, centerElement)

def _ifftshift(x):
    centerElement = (len(x)+1)//2
    return np.roll(x, centerElement)

如果你想摆脱使用 ifftshift 怎么办并且仍然获得与执行 fft(ifftshift(x)) 相同的结果？

#Implementation of fft(ifftshift(x)) without ifftshift

def fft_shift(x):
    N = len(x)
    centerElement = (N+1)//2
    f = np.fft.fft(x)
    discreteFreqs = np.arange(N) #These are frequencies related to untransformed sequence
    angles = - 2 * np.pi * discreteFreqs * centerElement/N
    unitElements  = np.cos(angles)  + np.sin(angles)*1j #Euler's formula
    #in frequency space we just multiply with these roots of unity
    f = np.multiply(f, unitElements)
    return f

所以这不仅仅是一个游戏，因为如果你正确旋转 fft 的结果，你甚至不必使用 ifftshift！当循环移位成本太高或者您希望将连续零的位置移动到 (N+1)//2 以外的位置时，它可能会派上用场。

检查时

ft = np.fft.fft(ifftshift(Y))
_ft = fft_shift(Y)
np.linalg.norm(ft-_ft) < 1e-10

可能会有一些数值错误，但结果实际上是相同的。