|
| 1 | +program forward_transform_of_real_function |
| 2 | + !! This program computes the forward transform of a real function and constructs |
| 3 | + !! the complex result (re)organized to match the array subscripting to the common |
| 4 | + !! analytical form [1]. Which form one uses in practice requires balancing the |
| 5 | + !! need for speed versus clarity. |
| 6 | + !! |
| 7 | + !! [1] https://docs.scipy.org/doc/scipy/reference/generated/scipy.fftpack.rfft.html#scipy.fftpack.rfft |
| 8 | + use fftpack, only: rfft, irfft |
| 9 | + implicit none |
| 10 | + integer j, k |
| 11 | + integer, parameter :: N = 8 |
| 12 | + double precision, parameter :: two_pi = 2.D0*acos(-1.D0), tolerance = 1.0D-06, f_avg = 3.D0, zero=0.D0 |
| 13 | + double precision, parameter :: f(0:N-1) = f_avg + [(cos(two_pi*dble(j)/dble(N)), j=0,N-1)] |
| 14 | + double precision, dimension(0:N-1) :: f_round_trip, rfft_f |
| 15 | + integer, parameter :: rk = kind(two_pi) |
| 16 | + complex(rk) f_hat(0:N/2) |
| 17 | + |
| 18 | + call assert(mod(N,2)==0, "the algorithm below requires even N") |
| 19 | + |
| 20 | + rfft_f(:) = rfft(f)/dble(N) |
| 21 | + f_hat(:) = [ cmplx(rfft_f(0),zero), [( cmplx(k,k+1), k=lbound(rfft_f,1)+1,ubound(rfft_f,1)-1,2)], cmplx(zero,rfft_f(N-1)) ] |
| 22 | + f_round_trip(:) = dble(irfft(rfft_f)) |
| 23 | + !call assert(any(abs(f_round_trip - f) < tolerance), "inverse of forward FFT must yield the original function") |
| 24 | + |
| 25 | + print *, "f = ", f |
| 26 | + print *, "f_hat = ", f_hat |
| 27 | + print *, "f_round_trip = ", f_round_trip |
| 28 | + |
| 29 | + !print '(3(10x,a,10x))',"f", "f_round_trip", "rfft_f" |
| 30 | + !do m = 1, size(f) |
| 31 | + ! print *, f(m), f_round_trip(m), rfft_f(m) |
| 32 | + !end do |
| 33 | + !print * |
| 34 | + |
| 35 | +contains |
| 36 | + pure subroutine assert(assertion, description) |
| 37 | + logical, intent(in) :: assertion |
| 38 | + character(len=*), intent(in) :: description |
| 39 | + if (.not. assertion) error stop description |
| 40 | + end subroutine |
| 41 | +end program |
0 commit comments