FFTR1: One-Dimensional Real Fourier Transforms

The trick to doing a DFT of real data is to once again use the DL formula. Specifically the routine begins by taking the original array and passing it as is to fftc1. Of course fftc1 treats the data as an array of complex numbers. Note, however, that the real parts of those numbers are the elements of $F^e$ and the imaginary parts are the elements of $F^o$, so it is possible to extract the actual DFT of $f$ from the results of fftc1. Let $H_b$ represent the output of fftc1 and let $F_b$ represent the desired output, i.e. the transform of the real function $f$. The input we gave to fftc1 was the combination $f^e + i f^o$ so by the linearity of Fourier transforms we know that $H = F^e + i F^o$. From the DL formula we also know that $F_b = F^e_b + W^b F^o_b$. How can we extract $F^e$ and $F^o$ from $H$ in order to plug them into the DL formula? The answer lies in the symmetry relation obeyed by the Fourier transforms of real data (which includes $F^e$ and $F^o$), $F_{-b} = F_b^*$. So

$$H_b = F^e_b + i F^o_b$$ (11)

$$H_{-b} = F^{e *}_b + i F^{o *}_b$$ (12)

and thus

$$F^e_b = {1 \over 2} \left(H_b + H_{-b}^*\right)$$ (13)

$$F^o_b = -{i \over 2} \left(H_b - H_{-b}^*\right)$$ (14)

Plugging these equations into the DL formula gives

$$F_b = {1 \over 2} \left[H_b + H_{-b}^* - i W^b \left(H_b - H_{-b}^*\right)\right].$$ (15)

Now all that's left is to figure out where all of these terms are stored in the array. Remember that all the array indices in these formulas are complex. In the case of $F$ finding the elements is trivial because only non-negative values of $k$ are given as output so the elements $F_b$ are stored in order $b=0,1,...,N/2-1$. (The case $b=N/2$ is treated separately below.) Recall, however, that $H$ is the transform of a complex $N/2$ size array, so its elements range from $b=-N/4$ to $N/4$ in the wrap-around storage order described in section 2.3.2. In other words the value of $H_{-b}$ is stored in the location $N/2-b$. So

$$F_b = {1 \over 2} \left[H_b + H_{N/2-b}^* - i W^b \left(H_b - H_{N/2-b}^*\right)\right].$$ (16)

Plugging in $N/2-b$ for $b$ in this equation and noting that $W^{N/2-b} = \exp[2 \pi i (N/2-b)/N] = -(W^b)^*$ gives

$$F_{N/2-b} = {1 \over 2} \left[H_b^* + H_{N/2-b} - i (W^b)^* \left(H_b^* - H_{N/2-b}^*\right)\right].$$ (17)

Using these two equations the program steps over all the elements of the array overwriting $H_b$ and $H_{N/2-b}$ with $F_b$ and $F_{N/2-b}$. The cases $b=0$ and $b=N/2$ are a little different, however, because $F_0$ and $F_{N/2}$ are both real. The function fftr1 places $F_0$ in the first (real) slot of the $b=0$ element of the array and $F_{N/2}$ in the second (imaginary) slot. From the formulas above, noting that $H_{N/2}=H_0$ by periodicity, we see that

$$F_0 = Re(H_0) + Im(H_0)$$ (18)

$$F_{N/2} = Re(H_0) - Im(H_0).$$ (19)

Inverting this process seems like a mess but actually works out to be nearly identical. I won't reproduce all the algebra here but if you invert the formulas above to find $H$ from $F$ you will find that except for a couple of minus signs and a factor of two for the $k=0$ and $k=N/2$ terms the results look identical to the forward case. As a result the inverse transform simply uses the same loop to rearrange the array and then calls fftc1 to do an inverse transform and get the function $f(x)$ back.