scipy 实现fftshift和ifftshift的正确顺序(在python中)

pbpqsu0x  于 2023-06-06  发布在  Python
关注(0)|答案(2)|浏览(271)

我想对一个函数psi(x)进行傅里叶变换,将其乘以一个k空间函数exp(-kx^2-ky^2),然后将乘积逆傅里叶变换回x空间。
但是我的x空间和k空间网格是居中的,我知道我需要fftshiftifftshift来正确实现我的k空间乘法。但我不明白它们是如何工作的,所以我不知道以何种顺序来实现它们。有人能告诉我我做得对不对吗?

import scipy.fftpack as spfft
import numpy as np

#Create a centred k-space grid]

kxmax, kymax = 10,10
kxgrid = np.linspace(-kxmax/2, kxmax/2, NX)
kygrid = np.linspace(-kymax/2, kymax/2, NY)
KX, KY = np.meshgrid(kxgrid, kygrid, indexing='xy')

psi = spfft.ifft2(spfft.fftshift(np.exp(-(KX**2 + KY**2)) * spfft.fftshift(spfft.fft2(psi))))
qaxu7uf2

qaxu7uf21#

不,你没有,但没关系,它可以是非常混乱的。
首先fftifft要求原点在向量的开头(或者在2D情况下,在数组的左上角)。输入psi的原点是否像KX一样居中?如果是这样,它的原点必须移动到ifftshift的开头。(如果没有,那就别管它了。
第二:由于KXKY的原点在它们的中心,所以你必须取消它们的移位:您需要spfft.ifftshift(np.exp(-(KX**2 + KY**2))(注意i)。
最后:输出psi因此将在开始时具有其原点。如果你想让它的原点像KX一样居中,fftshift它。
总结如下:

inputOriginStart = # ...
inputOriginStartFFT = spfft.fft2(psiOriginStart)
filterOriginStartFFT = spfft.ifftshift(np.exp(-(KX**2 + KY**2)))
outputOriginStart = spfft.ifft2(filterOriginStartFFT * inputOriginStartFFT)

其中inputOriginStart是输入psi,假设它的原点在开始处,其中outputOriginStart是输出psi-为清楚起见重新命名。(我总是追求清晰。如果它不起作用,你可以更容易地找出它。)

编辑修复了提问者指出的错误-是的,我有一个错误,在开始时留下了psiOriginStart的原点;则ifftshiftKXKY的中心原点函数。(如果你想取消将outputOriginStart的原点移到中心,那么使用fftshift。)
编辑2将过滤器(KXKY的函数)与数据分开,以使正确的括号显而易见。

如何保持这些直线?有几个技巧要记住:

  • fftifft总是需要输入并给予其起源在开始处的输出。这应该很容易从经验中记住。
  • fftshift获取fft需要/生成的起始原点,并将原点移到中心。同样,我倾向于很容易记住这一点,因为输入fftshift(fft(...))一千次的肌肉记忆。
  • 最后,唯一剩下的事情就是推导出ifftshiftfftshift的逆:它采用中心原点向量/数组,并将原点移到起点。
j9per5c4

j9per5c42#

对信号或频谱进行移位的需要由DFT/IDFT施加的约束引起。这些函数在其输入中具有对参考的期望,样本对应于时间t=0(对于DFT)和频率f=0(对于IDFT)。这两个函数也使用相同的约定输出结果。在介绍移位/逆移位函数之前,让我们先介绍一下DFT/IDFT的预期。

DFT约束

  • 周期性约束 *

假设输入序列是信号的完整周期,并且可以并置以创建周期信号,例如100个样本:

该信号的任何100个样本部分可以用于进行DFT。以上面的洋红标记开始的部分为例:

DFT幅度不改变,但在相位中考虑时间差(最后一行)。对于DFT,信号中的延迟导致对应的频谱乘以复指数exp(j.2 π f0.t)。由于幅度为1,因此乘积仅将2 π f0与频谱系数角相加,因此是相位差。

  • 排序约束 *

这意味着DFT评估延迟:它假设第一个样本用于t=0,并计算相对于它的相位。DFT的输出本身是周期性的,周期为N,N是一个周期中的样本数:

(Only这里示出了幅度,但是相位具有相同的周期性。)

IDFT约束

像以前一样,我们可以取前一序列的任何100个样本部分并进行IDFT:

IDFT具有与DFT相同的输入约束:假设输入频谱是周期性的,第一个样本是针对f=0(DC)的。频谱中的任何频移对应于信号乘以幅度为1的复指数。对于DFT,这仅改变输出的相位。

靠近水平刻度

对应于样本k(k在0和N-1之间)的时间是t=k/N.Ts,Ts是采样周期,并且频率样本k对应于频率k/N.Fs,fs是采样频率。
因此,在DFT的输出中不存在负频率,并且在IDFT的输入中不期望负频率。
然而,由于DFT和IDFT的固有周期性,我们被允许将频谱从范围0到fs转换到范围-1/2fs到+1/2fs。现在我们有一个频谱,其中一半的频率是负的,但仍然位于序列中的正频率之后。
函数fftfreq可用于创建具有正频率和负频率的刻度。该函数与DFT/IDFT没有直接联系,它只是一个帮助器,知道如何确定所需的频率值。
由于周期性N. Ts,时间范围可以以相同的方式从0到Ts移位到-1/2Ts到+1/2Ts。

使用fftshiftifftshift

现在让我们来回答你关于移位的问题。
带有半正和半负索引的比例便于绘制对称视图。但是有一个小问题,负指数放在0和正指数之后。
fftshiftifftshift可用于对元素重新排序:fftshift准备用于绘图目的的序列,ifftshift恢复DFT/IDFT使用/预期的原生顺序,并在第一部分中进行了描述。注意,这些函数除了重新排序元素之外不执行其他动作,它们与FT没有直接关系,尽管它们的名称包含'fft',并且前缀'i'与IDFT没有联系。
该重新排序操作可以被视为序列内的旋转,或者被视为沿着周期性序列(信号或频谱)的线性移位。
人们可能会感到惊讶,有两个功能可以做同样的改变。具有偶数长度的序列具有相同数量的正(包括DC/0)和“负”频率。对于奇数长度序列,存在一个额外的负频率。两个“一半”的长度不一样,这一点必须考虑在内。因此,fftshift可以通过以下方式管理大小:

for k in axes:
    n = tmp.shape[k]
    p2 = (n+1)//2
    mylist = concatenate((arange(p2, n), arange(p2)))
    y = take(y, mylist, k)
return

注意负半部分的位置是如何计算的:p2 = (n+1)//2。和ifftshift是相同的,除了它使用以下位置:p2 = n-(n+1)//2
由于只有用户知道序列是否已经被交换,因此他/她负责使用适当的功能。
底线是:通常不需要移动或取消移动序列。DFT和IDFT期望并输出具有索引0的序列作为第一样本。一种可能的情况是当序列被创建为居中,并且必须由DFT/IDFT处理时。在用作输入之前,必须使用ifftshift进行未移位。
为了对称显示的目的,可能需要移位(fftshift)以使样本居中为0,DC频率或时间t=0。如果移位后的序列必须在DFT/IDFT函数中重用,则可以恢复原始顺序(ifftshift)。仔细跟踪序列当前是否被移位,因为这决定了要使用哪个函数。

相关问题