高效计算傅里叶逆变换

Question

在这个问题的非常好的答案here.中建议了以下代码

我的问题是：如果我只是通过评估 np.dot( B, yl ) 来计算傅里叶空间中的导数，就像这段代码中所做的那样。如何通过应用傅里叶逆变换恢复实际空间中的实际导数？

import numpy as np

## non-normalized gaussian with sigma=1
def gauss( x ):
    return np.exp( -x**2 / 2 )

## interval on which the gaussian is evaluated
L = 10
## number of sampling points
N = 21
## sample rate
dl = L / N
## highest frequency detectable
kmax= 1 / ( 2 * dl )

## array of x values
xl = np.linspace( -L/2, L/2, N )
## array of k values
kl = np.linspace( -kmax, kmax, N )

## matrix of exponents
## the Fourier transform is defined via sum f * exp( -2 pi j k x)
## i.e. the 2 pi is in the exponent
## normalization is sqrt(N) where n is the number of sampling points
## this definition makes it forward-backward symmetric
## "outer" also exists in Matlab and basically does the same
exponent = np.outer( -1j * 2 * np.pi * kl, xl ) 
## linear operator for the standard Fourier transformation
A = np.exp( exponent ) / np.sqrt( N )

## nth derivative is given via partial integration as  ( 2 pi j k)^n f(k)
## every row needs to be multiplied by the according k
B = np.array( [ 1j * 2 * np.pi * kk * An for kk, An in zip( kl, A ) ] )

## for the part with the linear term, every column needs to be multiplied
## by the according x or--as here---every row is multiplied element 
## wise with the x-vector
C = np.array( [ xl * An for An in  A ] )

## thats the according linear operator
D = B + C

## the gaussian
yl = gauss( xl )

## the transformation with the linear operator
print(  np.dot( D, yl ).round( decimals=9 ) ) 
## ...results in a zero-vector, as expected

Answer 1

这里只需要定义逆变换的线性算子即可。 按照发布代码的结构，它将是

expinv = np.outer( 1j * 2 * np.pi * xl, kl ) 
AInv = np.exp( expinv ) / np.sqrt( N )

导数的傅里叶变换为

dfF = np.dot( B, yl )

使得实空间中的导数为

dfR = np.dot( AInv, dfF ) 

最终，这意味着“导数”运算符是

np.dot( AInv, B )

这是一个小的N（除了边缘）一个带有条目(-1,0,1)的三对角矩阵，即一个经典的对称数值导数。 增加N，它首先变为1,-2,0,2,1，即导数的更高阶近似。

最终，得到(..., d,-c, b,-a,0,a,-b, c,-d, ...) 类型的交替加权导数