Storage sequence of rfft result? #44
Description
Could someone help me understand the storage sequence of the rfft function result and how to get the actual complex result from the forward FFT of a real data? I wrote the test program below to transform a 3 + cos(x). I'm unclear why the rfft result is real so I tried using transfer to get the corresponding complex values. The result is close to what I'd expect but not quite.
program main use fftpack, only: rfft, irfft implicit none integer m integer, parameter :: N = 8 double precision, parameter :: two_pi = 2.D0*acos(-1.D0), tolerance = 1.0D-06, x_avg = 3.D0 double precision, parameter :: x(*) = x_avg + [(cos(two_pi*dble(m)/dble(N)), m=0,N-1)] double precision, allocatable, dimension(:) :: rfft_x, x_round_trip double complex, allocatable :: x_hat(:) rfft_x = rfft(x)/dble(N) if (size(rfft_x) /= size(x)) error stop "rfft_x/x size mismatch" x_round_trip = real(irfft(rfft_x)) if (size(x_round_trip) /= size(x)) error stop "x_round_trip/x: size mismatch" if (any(abs(x_round_trip - x) > tolerance)) error stop "round trip out of tolerance" print '(3(10x,a,10x))',"x", "x_round_trip", "rfft_x" do m = 1, size(x) print *, x(m), x_round_trip(m), rfft_x(m) end do print * x_hat = transfer(rfft_x, mold=x_hat, size=size(x)/2) print *, "x_hat = ", x_hat end program which gives the output
x x_round_trip rfft_x 4.0000000000000000 4.0000000000000000 3.0000000000000000 3.7071067811865475 3.7071068286895752 0.50000000000000000 3.0000000000000000 3.0000000000000000 -2.7755575615628914E-017 2.2928932188134525 2.2928931713104248 0.0000000000000000 2.0000000000000000 2.0000000000000000 -0.0000000000000000 2.2928932188134521 2.2928931713104248 0.0000000000000000 3.0000000000000000 3.0000000000000000 -2.7755575615628914E-017 3.7071067811865475 3.7071068286895752 0.0000000000000000 x_hat = (3.0000000000000000,0.50000000000000000) (-2.77555756156289135E-017,0.0000000000000000) (-0.0000000000000000,0.0000000000000000) (-2.77555756156289135E-017,0.0000000000000000) where I would have expected like
x_hat = (3., 0.) (0.5, 0.) (0.,0.) (0.,0.) and a corresponding reordering of rfft_x.