Оценка быстрого преобразования Фурье гауссовской функции в FORTRAN с использованием библиотеки FFTW3

Question

Я пытаюсь написать код ФОРТРАН для оценки быстрого преобразования Фурье гауссовой функции f(r)=exp(-(r^2)) с использованием библиотеки FFTW3.Как все знают, преобразование Фурье гауссовской функции является еще одной гауссовской функцией.

Я рассматриваю оценку интеграла преобразования Фурье от функции Гаусса в сферической координате.

Следовательно, полученный интеграл можно упростить до интеграла [r*exp(-(r^2))*sin(kr)]dr.

Я написал следующий код на языке FORTRAN, чтобы оценить DST дискретного преобразования SINE, которое представляет собой DFT дискретного преобразования Фурье, используя реальный входной массив PURELY.DST выполняется C_FFTW_RODFT00, существующим в FFTW3, с учетом того, что дискретные значения в пространстве позиций являются r = i * delta (i = 1,2, ..., 1024), а входной массив для DST равенфункция r*exp(-(r^2)) НЕ гауссова.Синусоидальная функция в интеграле [r*exp(-(r^2))*sin(kr)]dr является результатом ИНТЕГРАЦИИ по сферическим координатам, и это НЕ мнимая часть exp(ik.r), которая появляется, когда берется аналитическое преобразование Фурье в целом.

Однакорезультат не является гауссовой функцией в импульсном пространстве.

Module FFTW3
 use, intrinsic :: iso_c_binding
include 'fftw3.f03'
end module  

program sine_FFT_transform
use FFTW3
implicit none
integer, parameter :: dp=selected_real_kind(8)

real(kind=dp), parameter :: pi=acos(-1.0_dp)
integer, parameter :: n=1024 
real(kind=dp) :: delta, k
real(kind=dp) :: numerical_F_transform
integer :: i
type(C_PTR) ::  my_plan
real(C_DOUBLE), dimension(1024) :: y
real(C_DOUBLE), dimension(1024) :: yy, yk
integer(C_FFTW_R2R_KIND) :: C_FFTW_RODFT00

my_plan= fftw_plan_r2r_1d(1024,y,yy,FFTW_FORWARD, FFTW_ESTIMATE)

delta=0.0125_dp
do i=1, n        !inserting the input one-dimension position function
y(i)= 2*(delta)*(i-1)*exp(-((i-1)*delta)**2) 
! I multiplied by 2 due to the definition of C_FFTW_RODFT00 in FFTW3
end do

call fftw_execute_r2r(my_plan, y,yy)   
do i=2, n
k = (i-1)*pi/n/delta 
yk(i) = 4*pi*delta*yy(i)/2  !I divide by 2 due to the definition of 
                            !C_FFTW_RODFT00
numerical_F_transform=yk(i)/k
write(11,*) i,k,numerical_F_transform
end do
call fftw_destroy_plan(my_plan)

end program 

Выполнение предыдущего кода дает следующий график, который не предназначен для функции Гаусса. Может кто-нибудь помочь мне понять, в чем проблема?Я думаю, что проблема в основном из-за FFTW3.Возможно я не использовал это должным образом особенно относительно граничных условий.

Answer 1

Просмотр соответствующих страниц на сайте FFTW ( Преобразование из реального в реальное * , Виды преобразования , Реальное-нечетное DFT (DST) ) изаголовочный файл для Фортрана, похоже, что FFTW ожидает FFTW_RODFT00 и т. д., а не FFTW_FORWARD для указания типа преобразования из реального в реальное.Например,

! my_plan= fftw_plan_r2r_1d( n, y, yy, FFTW_FORWARD, FFTW_ESTIMATE )
my_plan= fftw_plan_r2r_1d( n, y, yy, FFTW_RODFT00, FFTW_ESTIMATE )

выполняет дискретное синусоидальное преобразование типа I (DST-I), показанное на приведенной выше странице.Эта модификация, похоже, решает проблему (т. Е. Делает преобразование Фурье гауссианом с положительными значениями).

---

Ниже приведена слегка измененная версия кода OP для эксперимента с вышеуказанной модификацией:

! ... only the modified part is shown...
real(dp) :: delta, k, r, fftw, num, ana
integer :: i, j, n
type(C_PTR) ::  my_plan
real(C_DOUBLE), allocatable :: y(:), yy(:)

delta = 0.0125_dp ; n = 1024   ! rmax = 12.8
! delta = 0.1_dp    ; n = 128    ! rmax = 12.8
! delta = 0.2_dp    ; n = 64    ! rmax = 12.8
! delta = 0.4_dp    ; n = 32    ! rmax = 12.8

allocate( y( n ), yy( n ) )

! my_plan= fftw_plan_r2r_1d( n, y, yy, FFTW_FORWARD, FFTW_ESTIMATE )
my_plan= fftw_plan_r2r_1d( n, y, yy, FFTW_RODFT00, FFTW_ESTIMATE )

! Loop over r-grid
do i = 1, n
    r = i * delta              ! (2-a)
    y( i )= r * exp( -r**2 )
end do

call fftw_execute_r2r( my_plan, y, yy )

! Loop over k-grid
do i = 1, n

    ! Result of FFTW
    k = i * pi / ((n + 1) * delta)    ! (2-b)
    fftw = 4 * pi * delta * yy( i ) / k / 2   ! the last 2 due to RODFT00

    ! Numerical result via quadrature
    num = 0
    do j = 1, n
        r = j * delta
        num = num + r * exp( -r**2 ) * sin( k * r )
    enddo
    num = num * 4 * pi * delta / k

    ! Analytical result
    ana = sqrt( pi )**3 * exp( -k**2 / 4 )

    ! Output
    write(10,*) k, fftw
    write(20,*) k, num
    write(30,*) k, ana
end do

Компиляция (с gfortran-8.2 + FFTW3.3.8 + OSX10.11):

$ gfortran -fcheck=all -Wall sine.f90 -I/usr/local/Cellar/fftw/3.3.8/include -L/usr/local/Cellar/fftw/3.3.8/lib -lfftw3

Если мы используем FFTW_FORWARD, как в исходном коде, мы получаем

с отрицательным лепестком (где fort.10, fort.20 и fort.30 соответствуют FFTW, квадратуре и аналитическим результатам).Изменение кода для использования FFTW_RODFT00 изменяет результат, как показано ниже, поэтому модификация, кажется, работает (но см. Определение сетки ниже).

---

Дополнительные примечания

• Я немного изменил определение сетки для r и k в моем коде (строки (2-a) и (2-b)), который находится вулучшить точность.Но я все еще не уверен, соответствует ли приведенное выше определение определению, используемому FFTW, поэтому, пожалуйста, прочитайте руководство для деталей ...

• Заголовочный файл fftw3.f03 предоставляет интерфейс для fftw_plan_r2r_1d

  type(C_PTR) function fftw_plan_r2r_1d(n,in,out,kind,flags) bind(C, name='fftw_plan_r2r_1d')
    import
    integer(C_INT), value :: n
    real(C_DOUBLE), dimension(*), intent(out) :: in
    real(C_DOUBLE), dimension(*), intent(out) :: out
    integer(C_FFTW_R2R_KIND), value :: kind
    integer(C_INT), value :: flags
  end function fftw_plan_r2r_1d
  

• (из-за отсутствия поддержки Tex эта часть очень некрасива ...) Интеграл от 4 pi r^2 * exp(-r^2) * sin(kr)/(kr) для r = 0 -> бесконечный равен pi^(3/2) * exp(-k^2 / 4) (полученоиз Wolfram Alpha или отметив, что на самом деле это трехмерное преобразование Фурье для exp (- (x ^ 2 + y ^ 2 + z ^ 2)) на exp (-i * (k1 x +)k2 y + k3 z)) с k = (k1, k2, k3)).Таким образом, хотя и немного нелогично, результат становится положительным гауссовым.

• Я думаю, что r-сетка может быть выбрана гораздо более грубой (например, delta до 0,4), что дает почтис той же точностью, пока она охватывает частотную область преобразованной функции (здесь exp(-r^2)).

Answer 2

Конечно, существуют отрицательные компоненты действительной части БПФ ограниченного гауссовского спектра.Вы просто используете реальную часть преобразования.Так что ваш сюжет абсолютно правильный.

Вы, кажется, ошибаетесь с реальной величиной, которая, конечно, не будет отрицательной.Для этого вам нужно будет fftw_plan_dft_r2c_1d, а затем вычислить абсолютные значения комплексных коэффициентов.Или вы можете ошибочно принять преобразование Фурье с ограниченным ДПФ.

Вы можете проверить здесь, чтобы убедиться в правильности ваших расчетов выше:

http://docs.mantidproject.org/nightly/algorithms/FFT-v1.html

Пожалуйста, имейте в виду, что графики на приведенной выше страницесдвинуты, так что частота 0 находится в середине спектра.

Ссылаясь на себя, числовая интеграция [r*exp(-(r^2))*sin(kr)]dr будет иметь отрицательные компоненты для всех k>1, если нормализуется до 0 для самой высокой частоты.

TLDR: Ваш график является абсолютным состоянием искусства ивстроенный в дискретный и ограниченный функциональный анализ.