FFTW: Проблемы с преобразованием реального в сложное и сложное в реальное 2D

Как указано в заголовке, я использую FFTW (версия 3.2.2) с Fortran 90/95 для выполнения 2D FFT реальных данных (фактически поля случайных чисел). Я думаю, что шаг вперед работает (по крайней мере, я получаю некоторый результат). Однако я хотел проверить все, выполнив 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 . Мне интересно, связана ли проблема с индексацией, по крайней мере, в этом направлении я склоняюсь. В любом случае, если бы кто-нибудь мог предложить какую-то помощь, это было бы замечательно!

Спасибо!


person JRC    schedule 22.09.2011    source источник
comment
Следует помнить, что когда вы выполняете БПФ, за которым следует ОБПФ, результат будет масштабироваться с постоянным коэффициентом относительно исходного ввода. Значение коэффициента масштабирования зависит от того, как определены/реализованы БПФ и ОБПФ, но обычно это N, где N — количество точек в БПФ. Конечно, у вас могут быть и другие, не связанные с этим проблемы.   -  person Paul R    schedule 23.09.2011
comment
Я понял проблему - она ​​совершенно не связана с вызовами БПФ. В моем коде (не в приведенном выше примере) я сделал k одномерным массивом, но обращался к нему как к скалярному значению индекса. В приведенном выше коде не должно быть ничего плохого. Упс! Просто чтобы это все еще было полезно для других людей, я нашел в Интернете код, который содержал много очень полезных примеров вызова подпрограмм FFTW с Fortran: people.sc.fsu.edu/~jburkardt/f_src/fftw3/fftw3_prb.f90   -  person JRC    schedule 23.09.2011
comment
@Jen В приведенном выше коде вы теряете точность, смешивая реальные * 8 и сложные * 8 переменные: комплекс * 8 соответствует реальному * 4, а не реальному * 8. Вместо этого используйте комплекс*16.   -  person ev-br    schedule 24.09.2011
comment
Джен, я даже не смог заставить этот пример кода работать. Вы включили что-нибудь в начале, чтобы оно работало с gfortran? Он также вызывает файл, который не существует в текущей версии, поскольку существует только оболочка .f, а не .f90.   -  person    schedule 09.12.2011
comment
@ Томас Джеймс, я не пробовал пример с gfortran - вместо этого я использую xlf. Вы уверены, что у вас есть правильные команды/синтаксис для ваших шагов компиляции? Обертка .f должна быть в порядке (замените .f90). Если хотите, я могу выслать вам мой Makefile, который я использовал для компиляции кода - просто дайте мне знать. Я пытался вставить это здесь, но это не подходит в качестве комментария.   -  person JRC    schedule 10.12.2011


Ответы (1)


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

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

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

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

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

person BenBoulderite    schedule 11.04.2018
comment
Вы правы, это проблема. Но вы не правы, рекомендуя, что kind=8 — это хорошо и что это всегда означает 8-байтовые компоненты. См. stackoverflow.com/questions/838310/fortran-90-kind-parameter Некоторые компиляторы отказываются компилировать kind=8. - person Vladimir F; 12.04.2018
comment
Спасибо. Мой ответ действительно был полон неточностей. Надеюсь, это лучше. - person BenBoulderite; 12.04.2018
comment
Это намного лучше. - person Vladimir F; 12.04.2018