FFTW: проблема с реальными к сложным и от сложных к реальным 2D трансфорсмам - PullRequest
Новая
       59

FFTW: проблема с реальными к сложным и от сложных к реальным 2D трансфорсмам

2 голосов
/ 22 сентября 2011

Как гласит заголовок, я использую FFTW (версия 3.2.2) с Fortran 90/95 для выполнения 2D БПФ реальных данных (фактически поле случайных чисел).Я думаю , что шаг вперед работает (по крайней мере, я получаю некоторый вывод).Однако я хотел проверить все, выполнив IFFT, чтобы посмотреть, смогу ли я восстановить исходные данные.К сожалению, когда я призываю комплекс к реальной рутине, ничего не происходит, и я не получаю вывод ошибок, поэтому я немного запутался.Вот некоторые фрагменты кода: 

implicit none

include "fftw3.f"

! - im=501, jm=401, and lm=60

real*8    :: u(im,jm,lm),recov(im,jm,lm)
complex*8 :: cu(1+im/2,jm)
integer*8 :: planf,planb    
real*8    :: dv

! - Generate array of random numbers
dv=4.0
call random_number(u)
u=u*dv
recov=0.0

k=30

! - Forward step (FFT)

call dfftw_plan_dft_r2c_2d(planf,im,jm,u(:,:,k),cu,FFTW_ESTIMATE)
call dfftw_execute_dft_r2c(planf,u(:,:,k),cu)
call dfftw_destroy_plan(planf)

! - Backward step (IFFT)

call dfftw_plan_dft_c2r_2d(planb,im,jm,cu,recov(:,:,k),FFTW_ESTIMATE)
call dfftw_execute_dft_c2r(planb,cu,recov(:,:,k))
call dfftw_destroy_plan(planb)

Вышеприведенный шаг вперед, кажется, работает (r2c), но шаг назад, похоже, не работает.Я проверил это, различая массивы u и recov, которые в итоге не были равны нулю.Кроме того, максимальное и минимальное значения массива recov были равны нулю, что, по-видимому, указывает на то, что ничего не изменилось.

Я просмотрел документацию FFTW и основал свою реализацию на следующей странице http://www.fftw.org/fftw3_doc/Fortran-Examples.html#Fortran-ExamplesМне интересно, если проблема связана с индексацией, по крайней мере, это направление, в котором я склоняюсь.В любом случае, если кто-нибудь может предложить какую-то помощь, это было бы замечательно!

Спасибо!

1 Ответ

1 голос
/ 11 апреля 2018

Не уверен, что это корень всех проблем здесь, но способ, которым вы объявляете переменные, может быть виновником.

Для большинства компиляторов (это, очевидно, даже не стандарт), Complex*8 является старым синтаксисом с одинарной точностью: комплексная переменная занимает в общей сложности 8 байтов, разделенных между действительной и мнимой частью (4 + 4 байта) ).

[Изменить 1 после комментария Владимира Ф. к моему ответу, смотрите его ссылку для деталей:] По моему опыту (т.е. систем / компиляторов, которые я когда-либо использовал), Complex(Kind=8) соответствует объявлению комплексного числа двойной точности (a действительная и мнимая части, каждая из которых занимает 8 байт).

В любой системе / компиляторе Complex(Kind=Kind(0.d0)) должен объявлять комплекс двойной точности.

Короче говоря, у вашего сложного массива неправильный размер. Замените вхождения Real*8 и Complex*8 на Real(kind=8) и Complex(Kind=8) (или Complex(Kind=kind(0.d0)) для лучшей переносимости) соответственно.

...