Storage sequence of rfft result?

[rouson]

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.

[zoziha]

When you use rfft, the result is an array that is usually shorter than the original DFT output and contains only useful information about the spectrum of the real signal. This is very useful for applications such as audio, seismic waves, and other physical real signal data, because it both saves computing resources and retains enough information to reconstruct the original signal (through the inverse real Fourier transform irfft).

I tried to ask chatgpt about the rfft function in fftpack, and combined with rfft in scipy, I can confirm that rfft in fftpack outputs the correct results.
Chatgpt told me that rfft is slightly different for sequences of odd and even lengths, specifically:

1. For a real sequence of length N, the regular DFT will produce a spectrum of size N containing complex results.
2. The rfft function only returns DFT values corresponding to positive and zero frequencies, because the Fourier transform of real signals has conjugate symmetry, which means that the value at negative frequencies can be deduced from the value at positive frequencies without actually calculating.
3. If the input sequence length N is even, rfft returns a spectrum of length N/2 + 1 (including zero frequency components), and if N is odd, (N + 1)/2 spectrum components are returned.

Here is a simple example I wrote:

program example

    use fftpack
    implicit none

    real(rk) :: x(4) = [9, -9, 1, 3]

101 format(a, i2)
100 format(a, *(es9.2, :, ","))
    print 101, "length  : ", size(x)
    print 100, "rfft    : ", rfft(x)
    print 100, "fft     : ", fft(cmplx(x, kind=rk))
    print 100, "fftshift: ", fftshift(fft(cmplx(x, kind=rk)))

    print *, ""
    print 101, "length  : ", size(x) + 1
    print 100, "rfft    : ", rfft(x, 5)
    print 100, "fft     : ", fft(cmplx(x, kind=rk), 5)
    print 100, "fftshift: ", fftshift(fft(cmplx(x, kind=rk), 5))

end program
! length  :  4
! rfft    :  4.00E+00, 8.00E+00, 1.20E+01, 1.60E+01
! fft     :  4.00E+00, 0.00E+00, 8.00E+00, 1.20E+01, 1.60E+01, 0.00E+00, 8.00E+00,-1.20E+01
! fftshift:  1.60E+01, 0.00E+00, 8.00E+00,-1.20E+01, 4.00E+00, 0.00E+00, 8.00E+00, 1.20E+01
!
! length  :  5
! rfft    :  4.00E+00, 2.98E+00, 9.74E+00, 1.75E+01, 3.39E+00
! fft     :  4.00E+00, 0.00E+00, 2.98E+00, 9.74E+00, 1.75E+01, 3.39E+00, 1.75E+01,-3.39E+00, 2.98E+00,-9.74E+00
! fftshift:  1.75E+01,-3.39E+00, 2.98E+00,-9.74E+00, 4.00E+00, 0.00E+00, 2.98E+00, 9.74E+00, 1.75E+01, 3.39E+00

As we can see from the example, I performed rfft and fft on the real number sequences [9,-9,1,3] and [9,-9,1,3,0], and finally performed fftshift on the results of fft. Combining the results of rfft and fftshift, We can see:

1. Even number sequence [9,-9,1,3]: rfft contains the values of non-negative frequencies, not the 0 value. (4,0) is the 0 frequency value, 4 is included, and the value 16 of (16,0) is included as an extra value, giving [4,8,12,16].
2. Odd number sequence [9,-9,1,3,0]: rfft contains the values of non-negative frequencies. (4,0) is still the 0 frequency value, giving [4,2.98,9.74,1.75,3.39].
   Negative frequency values can be reconstructed from non-negative frequencies, rfft thus saving storage space.

[zoziha]

Based on the above logic, I run the even real number sequence you gave me:

program example

    use fftpack
    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(:)

101 format(a, i2)
100 format(a, *(es9.2, :, ","))
    print 101, "length  : ", size(x)
    print 100, "rfft    : ", rfft(x)
    print 100, "fft     : ", fft(cmplx(x, kind=rk))
    print 100, "fftshift: ", fftshift(fft(cmplx(x, kind=rk)))

end program
! length  :  8
! rfft    :  2.40E+01, 4.00E+00,-2.22E-16, 0.00E+00,-0.00E+00, 0.00E+00,-2.22E-16, 0.00E+00
! fft     :  2.40E+01, 0.00E+00, 4.00E+00,-2.22E-16, 0.00E+00, 0.00E+00, 0.00E+00,-4.44E-16, 0.00E+00, 0.00E+00, 0.00E+00, 2.22E-16, 0.00E+00, 0.00E+00, 4.00E+00, 4.44E-16
! fftshift:  0.00E+00, 0.00E+00, 0.00E+00, 2.22E-16, 0.00E+00, 0.00E+00, 4.00E+00, 4.44E-16, 2.40E+01, 0.00E+00, 4.00E+00,-2.22E-16, 0.00E+00, 0.00E+00, 0.00E+00,-4.44E-16

As can be seen from fftshift, (2.4, 0.0) is a 0 frequency value, we get 2.4 and ignore 0.0, we get rfft values from the results of fftshift in the order of 9, 11, 12, 13, 14, 15, 16, and 1 numbers. (The values -2.22E-16 and -4.44E-16 are considered to be close.)

[rouson]

My confusion is twofold:

1. For FFT packages that I've used previously (e.g., the FFT in the Numerical Recipes book series), transforming the discrete samples of a real function provides a complex-valued result, where the result takes advantage of the symmetries to which you're alluding, i.e., only the complex half-plane needs to be stored and the real N/2 mode can be stored in the imaginary part of the 0 mode. So I'm confused by the fact that rfft returns a real result rather than a complex result.
2. The information that I would expect to see in the complex result is contained in the real rfft result, but not in the storage order that I would expect. The example that I posted corresponds to transforming $3 + cos(x)$, which is equivalent to $3 e^0 + \frac{e^{ix} + e^{-ix}}{2}$. This means that for the Fourier coefficients are $f_0=3$ and $f_1=f_{-1}=0.5$, where $f_{-1}$ need not be stored due to symmetry . Thus, with a complex result in the half-plane, I would expect two non-zero complex Fourier coefficients with the first one ($f_0$) having a real part of 3 and imaginary part 0 and the second one ($f_1$) having a real part of 0.5 and an imaginary part of 0.

Instead, the rfft result is not complex and when I try to reinterpret it as complex using the transfer function and the above logic, the real part of $f_0$ is what I would expect it to be, but the imaginary part is not. In fact, the imaginary part of $f_0$ contains what I would expect to be in the real part of $f_1$.

[rouson]

PR #48 demonstrates how to compute the coefficients in comment 44 above.