Performance evaluation of the inverse real-valued fast Fourier transform on field programmable gate array platforms using open computing language

Regular computational pattern

When the inputs are real, the method proposed in Garrido, Parhi & Grajal (2009) first eliminates redundant operations in the flow graph of the inverse CFFT (highlighted in the blue regions of Fig. 1). After relocating the twiddle factors between stages, a 32-point flow graph of the RIFFT obtained using this method is shown in Fig. 2. This flow graph consists entirely of real datapaths, and the input samples are annotated with a letter (r or i) to indicate whether the value corresponds to the real or imaginary part of the output.

Based on Fig. 2, we first remove the $W^0$ blocks (marked by red boxes) and replace them with the “X” symbol to achieve a more regular distribution of twiddle factors. These $W^0$ blocks, highlighted in the red boxes of Fig. 1, are not considered twiddle factors in the conventional CFFT signal flow graph; instead, they merely indicate that input data should be stored at different output locations. This modification results in the regular butterfly architectures between stage 1 and $lo{g}_{2}N-3$ shown in Fig. 3, where all datapaths are real. Second, we introduce a butterfly counter in each stage. The modified flow graph is presented in Fig. 3, with different colors used to highlight the butterfly structures involved.

In this article, to enable throughput comparisons with off-the-shelf FFT designs (Intel, 2016c), we employ butterfly units capable of processing eight real data points at a time. In Fig. 4, between stage 1 and $lo{g}_{2}N-2$, we observe that processing eight data points simultaneously requires two types of butterfly structures. The first type, shown in Fig. 4A, is used when the real transforms end and the complex sub-transforms begin. The block marked “X” indicates that no multiplication is required, and the data points are instead stored directly at different output locations. The second type, shown in Fig. 4B, is used within the complex sub-transform. It is worth noting that the regular computation pattern is not unique. For example, butterfly units that process four real data points at a time could also be used to construct such patterns.

We develop a scheduling algorithm for these two types of butterfly units. The algorithm determines which type of butterfly unit should be used in each stage. The equation for calculating the butterfly unit type is given as follows:

(3) $$B/\left(1\ll P\right)$$where $N$ denotes the FFT size, $B$ denotes the butterfly number and $P$denotes the stage number. If the above equation evaluates to 0, a type I unit is used; otherwise, a type II unit is selected. For example, in Fig. 3, all butterflies in stage 2 yield an equation value of 0; therefore, only type I units are used in this stage. In stage 1, butterflies 0 and 1 also yield 0, and thus use type I units, while butterflies 2 and 3 yield 1, and therefore use type II units.

The overall structure of the proposed OpenCL kernel is shown in Fig. 5. Stage 0 and the last two stages have fixed computational patterns and can be hard-coded outside of the loop. The computations between stage 1 and $lo{g}_{2}N-3$ are carried out using the two types of butterfly units shown in Fig. 4. The for-loop consists of three parts: (a) the butterfly scheduling algorithm described above. (b) the FFT data access scheme. (c) the twiddle factor access scheme. If these three components are properly designed, the loop can be fully unrolled (#pragma unroll) by high-level synthesis tools to construct pipelined architectures. In this case, the initiation interval (Intel, 2016a) equals one, indicating that the pipeline has no data or memory dependencies (Intel, 2016b) and is capable of computing one sample per clock cycle.

In the flow graph, we use a decimation-in-time (DIT) radix-8 RFFT to illustrate the idea, but similar approaches can also be applied to the decimation-in-frequency RFFT algorithm. Moreover, once the fixed computational pattern is extracted, other radices of butterflies, such as radix-2, can also be employed for computation, with corresponding data access schemes and twiddle factor access schemes developed accordingly.

FFT data access scheme

We directly adopt the data access scheme described in Intel (2016c). A sliding window and data reordering algorithm are designed to provide the butterfly unit with the correct FFT data in each clock cycle. The sliding window is essentially a buffer of length $N+8\ast (lo{g}_{2}N-2)$ and the data stored in this buffer are shifted by one element at each iteration. The reader is referred to Intel (2016c) for further details.

Twiddle factor access scheme

Twiddle factors are precomputed and stored in a ROM bank with consecutive addresses, as shown in Table 1. For example, $W^N$ is stored at the address N.

Based on the proposed scheduling algorithm and the structure of the butterfly units, we establish mathematical relationships among the stage number of the signal flow graph, the butterfly index in each stage, and the memory addresses of the twiddle factors. First, for the type I butterfly unit, the address of the upper twiddle factor is given by:

(4) $$B\ast \left(1\ll \left(lo{g}_{2}N-P-2\right)\right)$$where B denotes the butterfly number and P denotes the stage number. For example, in Fig. 3 for butterfly 3 in stage 2, according to Eq. (3) it is a type I butterfly unit and the above equation evaluates to 6. This means the upper twiddle factor in the type I unit is $W^6$ and the lower twiddle factor is $W^{12}$. Similarly, for butterfly 1 in stage 1, it is also a type I unit and the above equation evaluates to 4, indicating that the upper twiddle factor is $W^4$ and the lower twiddle factor is $W^8$. For the type II butterfly unit, the address of the upper twiddle factor is given by

(5) $$\left(\mathrm{B}\mathrm{\%}\left(1\ll \mathrm{P}\right)\right)\ast \left(1\ll \left(lo{g}_{2}N-P-1\right)\right).$$

For example, for butterfly 3 in stage 1, it is a type II butterfly unit and the above equation evaluates to 8. This means both the upper and lower twiddle factor in the type II unit is $W^8$.

Experiment setup

The importance of FFT in DSP applications has led to a number of well-known benchmarks.

CPU benchmark: FFTW (Massachusetts Institute of Technology, 2016) is a widely used free-software library that computes the discrete Fourier transform and its various special cases on scalar processors like CPUs. Its performance is competitive with vendor-optimized programs, but unlike these programs, FFTW is not tuned to a fixed machine. FFTW provides separate and specialized routines for calculating both CFFT and RFFT.

GPU benchmark: The CUFFT (Nvidia Developer, 2016) is the NVIDIA CUDA FFT product. It provides a simple interface for computing FFTs on an NVIDIA GPU, which allows users to quickly leverage the floating-point power and parallelism of the GPU. Contrary to FFTW, CUFFT is designed for vector processors like GPUs and hence it is much faster than FFTW. CUFFT offers several computation modes: C2C (complex to complex), C2R (complex to real) and R2C (real to complex). C2C is the regular complex FFT, C2R is the RFFT algorithm and R2C is the inverse of the RFFT, which could be implemented using a similar strategy presented in this article.

FPGA benchmark: Intel has an OpenCL implementation of one-dimensional CFFT (Intel, 2016c). The example processes 4,096 complex single-precision floating-point values. The input data is ordered and the output data is in bit-reversed order. The example contains a single radix-4 butterfly unit capable of processing 4 complex data points (eight real data points) per clock cycle. In order to compare with this benchmark, in our proposed design we have also chosen to use butterfly units capable of calculating eight real data points per clock cycle.

The OpenCL kernel which contains the proposed algorithm is pre-compiled using Intel’s OpenCL tool kits. The OpenCL host program is running on the CPU and the kernel is running on the FPGA. When compiling the kernel, three optimization techniques are used: (1) Loop unrolling (Intel, 2016d): loop unrolling involves replicating a loop body multiple times and decreases the number of iterations at the expense of increased hardware resource consumption. We use the “#pragma unroll” directive to unroll the loop in Fig. 4 and the compiler automatically builds the pipeline. The data or memory dependency (Intel, 2016b) is removed and the initialization interval is equal to one (Intel, 2016a). (2) Optimize floating point operations: with the “–fpc” (Intel, 2016e) compilation flag turned on, the Intel OpenCL compiler tries to remove floating point rounding operations and conversions whenever possible. (3) Non-interleaved memory (Intel, 2016f): in our case, we want to partition the banks manually as two non-interleaved (and contiguous) memory regions, using the “-no-interleaving” compilation flag, to achieve better load balancing. Specifically, in this work the kernel will read input data from one off-chip RAM bank and store the results in the other RAM bank. During computation all the intermediate data points are stored on the FPGA chip and the kernel does not have to interact with the off-chip RAM multiple times. By configuring FPGA memory in this way, the Von Neumann bottleneck is minimized.

Results

In this section, we compare the proposed RFFT with FFTW, CUFFT and Intel’s CFFT.

Performance tests of our proposed RFFT algorithm are carried out on a Nallatech PCIe-385N A7 FPGA board and it is compared to Intel’s CFFT on the same board. We also compare our implementation to FFTW and CUFFT tested on a laptop with an Intel Core i7-4510U CPU running at 2 GHz and a Nvidia GeForce GT 720M GPU running at 1.4 GHz. The kernel execution time of various FFT sizes and platforms is reported in Table 2. There are generally two ways to benchmark FFT. One approach is to use points/sec, defined by N/(time for one FFT in seconds). The other is to use flops/sec, defined by 5Nlog₂N/(time for one FFT in seconds). Our results match closely to Intel’s reported results in Intel (2016c). For a 4,096-point CFFT, Intel has reported a throughput of 81 GFLOPS on the BittWare S5-PCIe-HQ D5/D8 board while our result is 83 GFLOPS on a different board. For a 4,096 input size, in terms of points/sec, the proposed RFFT achieves approximately a 1.39x speedup compared with Intel’s CFFT, a 2.36x speedup compared with the RIFFT routine of CUFFT and a 24.91x speedup compared with the RIFFT routine of FFTW.

Xilinx has a document regarding their implementations of CFFT using HDL (Xilinx, 2011). The implementation details outlined in Xilinx (2011) and Salehi, Amirfattahi & Parhi (2013) differ from ours in several aspects. First, these two documents utilize Virtex family FPGAs while we utilize Stratix-V FPGAs. Second, these documents utilize a 16-bit or 20-bit word length to represent data while we employ a single-precision floating-point data format. Consequently, our calculated results can be effortlessly validated against FFTW or CUFFT. Despite these major differences, a rough comparison is still feasible. In Xilinx (2011), the maximum frequency is around 400 MHz for a 1,024-point CFFT, resulting in a theoretical throughput of 400 MPoints/sec. In contrast, according to Table 2 in Intel’s original implementation, the throughput for a 1,024-point CFFT is 1,259 MPoints/sec, approximately three times faster than XFFT. In Salehi, Amirfattahi & Parhi (2013), the article reports a throughput of 922 MPoints/sec for a 64-point RFFT using the radix-8 algorithm and two levels of parallelism. In our experiments, we employ a single instance of the kernel without replicating the pipeline (Intel, 2016g).

Details regarding the hardware resources available on the target FPGA board, including the count of DSP blocks and memory blocks, are provided in reference (Nallatech, 2016). The resource utilization ratio of the proposed RIFFT and Intel’s CFFT on the same FPGA is reported in Table 3, demonstrating that the RIFFT algorithm consumes fewer hardware resources for each input size. Table 4 shows the average power consumption in different configurations. The proposed method achieves a 2.92× better energy efficiency over CUFFT on the GPU, and a 18.98× better energy efficiency over FFTW on the CPU respectively. Within the specified input range, FPGA exhibits clear advantages over GPU and CPU in terms of both throughput and power efficiency.

But it is also important to realize that while OpenCL or C/C++ enables device portability, achieving true performance portability requires careful tuning and adaptation to the underlying hardware architecture. OpenCL is designed as a cross-platform programming model, allowing the same kernel code to be compiled and executed on CPUs, GPUs, FPGAs, DSPs, and other accelerators from vendors such as Intel, AMD, NVIDIA, and Xilinx. This is achieved by abstracting the underlying hardware through a standard API and runtime. However, while the same code can run across devices, it rarely achieves comparable performance everywhere, since efficiency depends strongly on hardware-specific characteristics. For instance, the proposed kernel structure leverages pipelined parallelism, which is well suited for FPGAs. In contrast, on GPUs the kernel should be restructured to exploit single instruction, multiple data (SIMD) parallelism in order to fully utilize the architecture.