有效地计算傅立叶逆变换 [英] Efficiently compute inverse Fourier transform

问题描述

在此问题的很好回答中建议使用以下代码在这里.

The following code was suggested in the very nice answer to this question here.

关于它的我的问题是:如果我只是通过评估np.dot(B，yl)来计算傅里叶空间中的导数，如本代码所示.通过应用傅立叶逆变换，我如何能够在真实空间中恢复实际导数?

My question about it is: If I was just to compute the derivative in Fourier space as done in this code by evaluating np.dot( B, yl ). How would I be able to recover the actual derivative in real space by applying an inverse Fourier transform?

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

推荐答案

这里只需要为逆变换定义线性运算符.按照已发布代码的结构，将是

Here one only needs to define the linear operator for the inverse transformation. Following the structure from the posted code it would be

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

导数的傅立叶变换是

dfF = np.dot( B, yl )

使得真实空间中的导数为

such that the derivative in real space would be

dfR = np.dot( AInv, dfF ) 

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

Eventually, this means that the "derivative"-operator is

np.dot( AInv, B )

对于小 N (除了边缘)而言，是一个三对角矩阵，其条目为(-1,0,1)，即经典的对称数值导数.增加 N 会首先更改为 1，-2,0,2,1 ，即导数的高阶近似.

which is for small N ( except for the edges ) a tridiagonal matrix with entries (-1,0,1), i.e. a classical symmetric numerical derivative. Increasing N it first changes to 1,-2,0,2,1 i.e. a higher order approximation of the derivative.

最终，一个人得到了(...，d，-c，b，-a，0，a，-b，c，-d，...)类型的交替加权导数>

Eventually, one gets an alternating weighted derivative of type (..., d,-c, b,-a,0,a,-b, c,-d, ...)

这篇关于有效地计算傅立叶逆变换的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！