在 Python 中使用 FFT 和 IFFT 在偶数点采样时,离散函数变得复杂,同时自由传播实数函数

问题描述 投票:0回答:1

我的代码涉及使用傅里叶变换和傅里叶逆变换传播实函数。具体来说,该函数演化为

∂ψ(z,t)/∂t - v ∂ψ(z,t)/∂z =0

我通过对上面给出的方程进行傅立叶变换来解决这个问题

∂ψ(k,t)/∂t + ikv ψ(k,t)=0  

给出了解决方案

ψ(k,t)=e^(-ikvt)ψ(k,0)

现在执行傅里叶逆变换将给出传播函数

ψ(z-vt,0)

但是,当我在 Python 中执行此操作时(请参阅下面的代码),最终函数(应用 FFT 和 IFFT 后)似乎有一个微小的复杂部分,特别是当函数在偶数个点采样时。这个错误在我的模拟过程中不断累积,导致答案不太准确。然而,当我以奇数个点对函数进行采样时,这个问题似乎就消失了。有人可以帮忙解决这里发生的事情吗?

下面是一个简单的代码,使用函数的范数说明了我的观点。运行这段代码,我们可以看到,当函数以 10 个点采样时,经过 FFT 和 IFFT 后函数实部的范数是准确的,有时只能精确到小数点后 2 或 3 位 (而整个函数(实数+复数)的范数仍然守恒)。另一方面,在奇数采样的情况下,范数精确到小数点后 16 位。如果有人可以向我解释发生了什么以及如何在偶数点采样时使用此函数获得更好的准确性,那就太好了

import numpy as np

def free_prop(vectr,Nd,vel=0.92,dt=1):
    dz=1/Nd
    psi_k = np.fft.fft(vectr)
    k_vals = 2.0*np.pi*np.fft.fftfreq(Nd, d=dz)
    return np.fft.ifft(np.exp(-1j*dt*k_vals*vel)*psi_k)

vectr_even=np.random.rand(10) ## Even case
vectr_odd=np.random.rand(11) ## Odd case

print('Norm of input array (with even sampling):',np.linalg.norm(vectr_even))
print('Norm of output complex array (with even sampling):',np.linalg.norm(free_prop(vectr_even,10)))
print('Norm of real part of output array (with even sampling):',np.linalg.norm(np.real(free_prop(vectr_even,10))))
print()
print('Norm of input array (with odd sampling):', np.linalg.norm(vectr_odd))
print('Norm of output complex array (with odd sampling):',np.linalg.norm(free_prop(vectr_odd,11)))
print('Norm of real part of complex array (with odd sampling):',np.linalg.norm(np.real(free_prop(vectr_odd,11))))

从这里可以看出,偶数情况下的范数(仅考虑实部时)在小数点后第二位有所不同,如果它是浮点错误,我会感到非常惊讶。

python fft wave propagation norm
1个回答
0
投票

k_vals

Nd=6

[0., 6.28318531, 12.56637061, -18.84955592, -12.56637061, -6.28318531]

Nd=5

[0., 6.28318531, 12.56637061, -12.56637061, -6.28318531]

Nd
为偶数时,存在一个负频率,没有对应的正频率。众所周知,IFFT 的输入必须是共轭对称的,输出才能是实值。对于奇数长度数组,计算出的指数将自动是共轭对称的。对于偶数长度数组,存在一个没有正对应项的频率。我们需要使该频率分量为实值,以便指数是共轭对称的,并且逆 FFT 是实值的。

def free_prop(vectr, Nd, vel=0.92, dt=1):
    dz=1/Nd
    psi_k = np.fft.fft(vectr)
    k_vals = 2 * np.pi * np.fft.fftfreq(Nd, d=dz)
    k_kernel = np.exp(-1j * dt * k_vals * vel)
    if np.mod(Nd, 2) == 0:
       k_kernel[Nd // 2] = k_kernel[Nd // 2].real
    return np.fft.ifft(k_kernel * psi_k)
© www.soinside.com 2019 - 2024. All rights reserved.