高效计算逆傅里叶变换

Question

在这个问题的非常好的答案中建议使用以下代码

我的问题是：如果我只是像这段代码中那样通过计算 np.dot( B, yl ) 来计算傅里叶 space 中的导数。我怎样才能通过应用傅立叶逆变换恢复实数 space 中的实际导数？

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 )

使得实数 space 的导数为

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, ...)

的交替加权导数