[edit 1] Добавлены рисунки, чтобы показать исходные данные и полученные данные
[править 2] Я нашел свои ошибки, я использовал fftw_measure вместо fftw_estimate в вызове dfftw_plan_many_dft
[edit 3] исправлена опечатка в коде (замените u на u2d в dfftw_execute_dft_r2c)
Я пытаюсь выполнить 2D fft массива, используя несколько 1D fft вместо использования функции 2D fft, уже существующей в библиотеке fftw. И впоследствии мне нужно выполнить обратное 2D БПФ. Причина этого в том, что (в будущем) мой массив будет слишком большим, чтобы его можно было загрузить за один раз для выполнения 2D fft.
1-й черновик моего кода на данный момент выглядит примерно так:
double precision u2d(nx,ny),u2d2(nx,ny)
double complex qhat2d(nx/2+1,ny),qhat2d2(nx/2+1,ny)
integer N(1)
integer howmany, idist, odist, istride, ostride
integer inembed(2), onembed(2)
integer rank
! some function to read the data into u2d
! perform x-fft
N(1) = NX
howmany = NY
inembed(1) = NX
inembed(2) = NY
istride = 1
idist = NX
ostride = 1
odist = (NX/2+1)
onembed(1) = (NX/2+1)
onembed(2) = NY
rank = 1
write(*,*) 'u', u2d(1,1)
CALL dfftw_plan_many_dft_r2c(PLAN,rank,N(1),howmany,
& u2d,inembed,
& istride,idist,
& qhat2d,onembed,
& ostride,odist,FFTW_ESTIMATE) !
CALL dfftw_execute_dft_r2c(PLAN,u2d,qhat2d) ! x-fft
CALL dfftw_destroy_plan(PLAN)
! perform y-fft
N(1) = NY
howmany = (NX/2+1)
inembed(1) = (NX/2+1)
inembed(2) = NY
istride = (NX/2+1)
idist = 1
ostride = (NX/2+1)
odist = 1
onembed(1) = (NX/2+1)
onembed(2) = NY
rank = 1
CALL dfftw_plan_many_dft(PLAN,rank,N(1),howmany,
& qhat2d,inembed,
& istride,idist,
& qhat2d2,onembed,
& ostride,odist,FFTW_FORWARD,
& FFTW_MEASURE) !
CALL dfftw_execute_dft(PLAN,qhat2d,qhat2d2) ! y-fft
CALL dfftw_destroy_plan(PLAN)
! normally here, perform some filtering operation
! but at the moment, I do nothing
! perform inv-y-fft
N(1) = NY
howmany = (NX/2+1)
inembed(1) = (NX/2+1)
inembed(2) = NY
istride = (NX/2+1)
idist = 1
ostride = (NX/2+1)
odist = 1
onembed(1) = (NX/2+1)
onembed(2) = NY
rank = 1
CALL dfftw_plan_many_dft(PLAN,rank,N(1),howmany,
& qhat2d2,inembed,
& istride,idist,
& qhat2d,onembed,
& ostride,odist,FFTW_BACKWARD,
& FFTW_MEASURE) !
CALL dfftw_execute_dft(PLAN,qhat2d2,qhat2d) ! inv-y-fft
CALL dfftw_destroy_plan(PLAN)
! perform inv-x-fft
N(1) = NX ! I'm not too sure about this value here
howmany = NY
inembed(1) = (NX/2+1)
inembed(2) = NY
istride = 1
idist = (NX/2+1)
ostride = 1
odist = NX
onembed(1) = NX
onembed(2) = NY
rank = 1
CALL dfftw_plan_many_dft_c2r(PLAN,rank,N(1),howmany,
& qhat2d,inembed,
& istride,idist,
& u2d2,onembed,
& ostride,odist,FFTW_ESTIMATE) !
CALL dfftw_execute_dft_c2r(PLAN,qhat2d,u2d2) ! x-fft
CALL dfftw_destroy_plan(PLAN)
write(*,*) 'u2d2', u2d2(1,1)
do i=1,nx
do j=1,ny
u2d2(i,j) = u2d2(i,j) / (nx*ny)
enddo
enddo
write(*,*) 'u2d2', u2d2(1,1) ! here the values u2d2(1,1) is different from u2d(1,1)
! some action to write u2d2 to file
end
Я ожидал, что u2d и u2d2 будут одинаковыми, но я получаю относительно разные значения. Я где-то ошибся?
Оригинал и результаты показаны ниже. Форма выглядит похожей, но значения относительно разные (например, минимум и максимум).
_r2c
) для сложных данных? - person Alexander Vogt   schedule 23.05.2016