High performance fourier transforms on GPUs with GLSL

Universidade do Minho Escola de Engenharia Departamento de Informática Jorge Francisco Teixeira Bastos da Mota High Performance Fourier Transforms on GPUs with GLSL February 2023 Universidade do Minho Escola de Engenharia Departamento de Informática Jorge Francisco Teixeira Bastos da Mota High Performance Fourier Transforms on GPUs with GLSL Master dissertation Integrated Master’s in Informatics Engineering Dissertation supervised by António José Borba Ramires Fernandes February 2023 COPYRIGHT AND TERMS OF USE FOR THIRD PARTY WORK This dissertation reports on academic work that can be used by third parties as long as the internationally accepted standards and good practices are respected concerning copyright and related rights. This work can thereafter be used under the terms established in the license below. Readers needing authorization conditions not provided for in the indicated licensing should contact the author through the RepositóriUM of the University of Minho. LICENSE GRANTED TO USERS OF THIS WORK: CC BY-SA https://creativecommons.org/licenses/by-sa/4.0/ a ACKNOWLEDGEMENTS First of all, I would like to thank my supervisor, Professor António Ramires, who accepted to guide me in this project, presenting on several occasions fundamental suggestions for the improvement of this dissertation, making the realization of this dissertation a true learning process. Furthermore, I am extremely grateful to my family for their unconditional support, affection, and interest in pushing me in the right direction throughout my life. I would also like to thank all my friends, who accompanied me along this path, who listened to my concerns and advised me when necessary. b STATEMENT OF INTEGRITY I hereby declare having conducted this academic work with integrity. I confirm that I have not used plagiarism or any form of undue use of information or falsification of results along the process leading to its elaboration. I further declare that I have fully acknowledged the Code of Ethical Conduct of the University of Minho. c ABSTRACT The Fast Fourier Transform is a family of algorithms indispensable for the computation of the Discrete Fourier Transform. As a result, these transforms are the core of many applications in several areas and are required to be computed efficiently in many scenarios. The continuous evolution of GPUs has increased the popularity of parallelizable algorithm implementations on this type of hardware. Traditionally GPUs were associated to graphics background, however, with the popularization of the compute functionality of this hardware, most modern GPUs now have this capability, hence, algorithms now are more likely to be implemented in the general-purpose compute pipeline of GPUs. As a result, many applications take advantage of compute programming in GPGPU-capable frameworks such as GLSL, a high-level shading language frequently used in the context of computer graphics. In this dissertation we provide, refine and compare GPU-driven implementations of the family of FFT algorithms in GLSL, with the goal to provide programmers with efficient and simplified compute kernels for this transform, from the classic Cooley-Tukey algorithm to more suitable algorithms for the GPU such as the Stockham algorithm with higher radix. Accordingly, we also use the cuFFT NVIDIA framework for reference in the comparisons of the GLSL algorithms implementations with the goal to analyse their significance on the tradeoff of using specialized implementations of the FFT algorithms or integrating dedicated software tools for any case of application. Finally, we demonstrate how all improvements discussed in this dissertation culminate in performance improvement in a real-time rendering technique that heavily depends on multiple of these transforms in the Nau3D engine as a case of study. K E Y W O R D S FFT, GPGPU, GLSL, cuFFT, performance, compute, Cooley-Tukey, Stockham. d RESUMO A Transformada Rápida de Fourier é um algoritmo ou uma família de algoritmos indispensáveis para o cálculo da Transformada Discreta de Fourier. Assim, essas transformadas são o núcleo de muitas aplicações em diversas áreas e precisam ser calculadas de forma eficiente em muitos cenários. A evolução contínua dos GPUs aumentou a popularidade das implementações de algoritmos paralelizáveis neste tipo de hardware. Tradicionalmente, os GPUs eram associadas ao fundo gráfico, no entanto, com a popularização da funcionalidade de compute desse hardware, os GPUs mais modernos agora têm essa capacidade, portanto, os algoritmos agora são mais propensos a serem implementados na compute pipeline de propósito geral dos GPUs. Como resultado, muitas aplicações aproveitam a programação em compute em frameworks compatíveis com GPGPU como GLSL, uma linguagem de shading de alto nível usada recorrentemente no contexto de computação gráfica. Nesta dissertação fornecemos, refinamos e analisamos implementações em GPU da família de algoritmos FFT em GLSL, com o objetivo de fornecer aos programadores compute kernels eficientes e simplificados para esta transformada, desde o clássico algoritmo de Cooley-Tukey até algoritmos mais adequados para o GPU. Da mesma forma, também usamos a framework cuFFT NVIDIA como referência nas comparações das implementações dos algoritmos em GLSL com o objetivo de analisar a sua importância no tradeoff entre usar implementações especializadas dos algoritmos FFT ou integrar ferramentas de software dedicadas para qualquer caso de aplicação. Por fim, demonstramos como todas as melhorias discutidas nesta dissertação culminam na melhoria de desempenho numa técnica de renderização em tempo real que depende de FFTs na engine Nau3D como caso de estudo. PA L AV R A S -C H AV E FFT, GPGPU, GLSL, cuFFT, performance, compute, Cooley-Tukey, Stockham. e CONTENTS 1INTRODUCTION 5 1.1 Contextualization 5 1.2 Motivation 5 1.3 Aim of the work 6 1.4 Document structure 6 2THE FOURIER TRANSFORM 7 2.1 Continuous Fourier Transform 8 2.2 Discrete Fourier Transform 9 2.2.1 Matrix multiplication 10 2.3 Fast Fourier Transform 12 2.3.1 Radix-2 Decimation-in-Time FFT 13 2.3.2 Radix-2 Decimation-in-Frequency FFT 17 2.4 Stockham algorithm 20 2.5 Radix-4 instead of Radix-2 22 2.6 Two real inputs within one complex 26 2.7 2D Fourier Transform 27 3IMPLEMENTATION IN GLSL 28 3.1 Cooley-Tukey 28 3.1.1 All stages in one pass 32 3.2 Radix-2 Stockham 33 3.3 Radix-4 Stockham 34 4ANALYSIS AND COMPARISON 38 4.1 cuFFT 38 4.2 GLSL implementation results 39 4.3 Case of study 43 4.3.1 Tensendorf waves 43 4.3.2 Results 44 5CONCLUSIONS AND FUTURE WORK 46 iii CONTENTS iv 5.1 Results 46 5.2 Future work 47 Bibliography 48 I APPENDICES A G L S L F F T 52 B C U F F T 67 2 THE FOURIER TRANSFORM It is noticeable the presence of Fourier Transforms in a great variety of apparent unrelated fields of application, even the FFT is often called ubiquitous due to its effective nature of solving a great hand of problems for the most intended time complexity. Some applications include polynomial multiplication Jia (2020), numerical integration, time-domain interpolation, and x-ray diffraction. Moreover, it is present in several fields of study such as Applied Mechanics, Signal Processing, Sonics and Acoustics, Biomedical Engineering, Instrumentation, Radar, Numerical Methods, Electromagnetics and more (Shakshi (2016), Lee et al. (2017), Brigham (1988)). As already mentioned before, this dissertation focuses more on the application of FFT in the context of Computer Graphics, where this transform is applied in several cases such as image processing, image filtering, and ocean waves simulation. In Signal Analysis when representing a signal with amplitude as a function of time, it can be translated to the frequency domain, a domain that consists of signals of sines and cosines waves of varied frequencies, as illustrated in Figure 1, but to calculate the coefficients of those waves we use the Fourier Transform. Figure 1: Time to frequency signal decomposition Source: NTiAudio Sine and cosine waves are in simple waveforms, they can then be manipulated with relative ease. This process is constantly present in communications due to the transmission of data over wires and radio circuits through signals and most devices nowadays perform it frequently. In this introductory chapter, we present a preface on the Fourier Transform in Section 2.1 and describe the discrete version of the Fourier Transform in Section 2.2, which is the focus of this dissertation. Furthermore, we present the state of the art of the most popular algorithms in Section 2.3 and progressively demonstrate how 7 2.1. Continuous Fourier Transform 8 to improve them in the following sections Section 2.4 and Section 2.5. Regarding optimization of the FFT, we then proceed to detail how a complex-to-complex FFT computation can simultaneously transform two distinct real input sequences in Section 2.6. At the end of this chapter, in Section 2.7 we finalize with an explanation of how one-dimensional algorithms can be applied in the context of 2D input sequences such as an image 2.1 CONTINUOUS FOURIER TRANSFORM The Fourier Transform is a mathematical method to decompose a function into frequency components. Intuitively, the Inverse Fourier Transform is the corresponding method to reverse that process and reconstruct the original function from the one in frequency domain representation. Although there are many forms, the Fourier Transform key definition (Adams (2020)) can be described as Equation 1. X(f) = Z+∞ −∞x(t)e−i f tdt x(t) = 1 2πZ+∞ −∞X(f)e−i f td f (1) where •X(f),∀f∈R→ function in frequency domain representation, also called the Fourier Transform of x(t); •x(t),∀t∈R→function in time or space domain representation; •i→imaginary unit i=√−1. This formulation shows the usage of a complex-valued domain, making the Fourier Transform range from real to complex values, one complex coefficient per frequency X:R→C If we take into account Euler’s formula (Equation 2), we can rewrite the Fourier Transform as represented in Equation 3. eix =cos x+isin x(2) X(f) = Z+∞ −∞x(t)(cos(−f t) + isin(−f t))dt (3) Hence, we can break the Fourier Transform apart into two formulas that give each coefficient of the sine and cosine components as functions without dealing with complex numbers. 2.2. Discrete Fourier Transform 9 X(f) = Xa(f) + iXb(f) Xa(f) = Z+∞ −∞x(t)cos(f t)dt Xb(f) = Z+∞ −∞x(t)sin(f t)dt (4) This model of the Fourier Transform applied to infinite domain functions is called Continuous Fourier Transform. 2.2 DISCRETE FOURIER TRANSFORM The Fourier Transform of a finite sequence of equally-spaced samples of a function is called the Discrete Fourier Transform (DFT). It converts a finite set of values in time/space domain to frequency domain representation. It is an usefull version of the Fourier Transform since it deals with a discrete amount of data. This transform is described in Equation 5. Xk= N−1 ∑ n=0 xn·e−i2π Nkn xn=1 N N−1 ∑ k=0 Xk·ei2π Nkn (5) Notably, the discrete version of the Fourier Transform has some obvious differences since it deals with a discrete time sequence. The first difference is that the sum covers all elements of the input values instead of integrating the infinite domain of the function, but we can also notice that the exponential, similar to the aforesaid, divides the values by N ( N being the total number of elements in the sequence) due to the inability to look at frequency and time f t continuously we instead take the k’th frequency over n. We can expand this formula as: Xk=x0+x1ei2π Nk+... +xN−1ei2π Nk(N−1) Having this sum simplified we then only need to resolve the complex exponential, and we can do that by replacing the ei2π Nkn by the Euler formula as mentioned before to reduce the maths to a simpler sum of real and imaginary numbers. Xk=x0+x1(cos b1+isin b1) + ... +xN−1(cos bN−1+isin bN−1)(6) where bn=2π Nkn Finally, the result will be a complex number. 2.2. Discrete Fourier Transform 10 EXAMPLE Let us now follow an example of calculation of the DFT for a sequence x with N number of elements. x=h1 0.707 0 −0.707 −1−0.707 0 0.707i N=8 With this sequence, we now want to transform it into the frequency domain, therefore we apply the Discrete Fourier Transform to each element xn→Xk , thus, for each k ’th element of X we apply the DFT for every element of x. X0=1·e−i2π 8·0·0+0.707 ·e−i2π 8·0·1+... +0.707 ·e−i2π 8·0·7 = (0+0i) X1=1·e−i2π 8·1·0+0.707 ·e−i2π 8·1·1+... +0.707 ·e−i2π 8·1·7 = (4+0i) ... X7=1·e−i2π 8·7·0+0.707 ·e−i2π 8·7·1+... +0.707 ·e−i2π 8·7·7 = (4+0i) And that will produce our complex-valued output in the frequency domain. X=h0i4+0i0i0i0i0i0i4+0ii 2.2.1 Matrix multiplication The example shown above is done sequentially as if each frequency component is computed individually, but there is a way to express the DFT using matrix multiplication (Rao and Yip (2018)). Since the operations are done equally without any extra step we can group all analyzing function sinusoids ( e−i2π Nkn ), also referred to as twiddle factors. W=        ω0·0 Nω1·0 N. . . ω(N−1)·0 N ω0·1 Nω1·1 N. . . ω(N−1)·1 N . . .. . ..... . . ω0·(N−1) Nω1·(N−1) N. . . ω(N−1)·(N−1) N        =       1 1 . . . 1 1ω. . . ω(N−1) . . .. . ..... . . 1ω(N−1). . . ω(N−1)·(N−1)       2.2. Discrete Fourier Transform 11 where ωN=e−i2π N The substitution variable ωallows us to avoid writing extensive exponents. The symbol W represents the transformation matrix of the Discrete Fourier Transform (Rao and Yip (2018)), also called the DFT matrix, and its inverse can be defined as. W−1=1 N·        1 1 . . . 1 1ωN. . . ω(N−1) N . . .. . ..... . . 1ω(N−1) N. . . ω(N−1)·(N−1) N        where ωN=ei2π N We separate the twiddle factors in a matrix and express the DFT as a matrix multiplication between the twiddle factors and the input sequence. When the size is known before the transformation, it is possible to pre-calculate and reuse the matrix of the factors X=W·x x=W−1·X It is also worth noting that normalizing the DFT and IDFT matrix be by √N instead of just normalizing the IDFT by N , will make W a unitary matrix Horn and Johnson (2012). However, this normalization by √N is not common in FFT implementations since it is more simple and efficient to normalize the inverse by N , than the forward and the inverse by √N. EXAMPLE Continuing the previous example, we can adapt the application of the DFT to the matrix multiplication form. W=       1 1 . . . 1 1ω8. . . ω7 8 . . .. . ..... . . 1ω7 8. . . ω49 8       where ω8=ei2π 8 X=W·x=W·       1 0.707 . . . 0.707       =       0 4+0i . . . 4+0i       2.3. Fast Fourier Transform 12 It is conspicuous that the complexity time for each complex multiplication of every singular term of the sequence with the complex exponential value is O(N2) , hence, the computation of the Discrete Fourier Transform rises exponentially with the sequence length. Therefore, over time new algorithms and techniques were developed to increase the performance of this transform due to its usefulness. 2.3 FAST FOURIER TRANSFORM The Fast Fourier Transform (FFT) is a family of algorithms that compute the Discrete Fourier Transform (DFT) of a sequence, and its inverse, efficiently, since the direct usage of the DFT formulation is too slow for its applications. Thus, FFT algorithms exploit the DFT matrix structure by employing a divide-and-conquer approach (Chu and George (1999)) to segment its application. Over time several variations of the algorithms were developed to improve the performance of the DFT and many aspects were rethought in the way we compute the transform. There are many algorithms and approaches on the FFT family such as the well known Cooley-Tukey Cooley and Tukey (1965), known for its simplicity and effectiveness to compute any sequence with size as a power of two, but also Rader’s algorithm Rader (1968) and Bluestein’s algorithm Bluestein (1970) to deal with prime sized sequences, and even the Split-radix FFT Yavne (1968) that recursively expresses a DFT of size N in one DFT of size N/2 and two DFTs of size N/4. Computing the DFT directly requires N2 complex multiplies and N(N−1) complex additions, by using an FFT algorithm to compute the DFT, it only requires (N/2)log (N)complex multiplies and complex additions. When we consider a DFT of a sequence of length 4, we can notice that some of the operations involved repeat themselves for the computation of the whole sequence. In Equation 7 we list the required operations to compute the DFT of a sequence with 4elements. X0=x0+x1+x2+x3 X1=x0−i·x1−x2+i·x3 X2=x0−x1+x2−x3 X3=x0+i·x1−x2−i·x3 (7) When the operations are grouped together, we have a clear perception of the repetition of the operations, as presented in Equation 8. 2.3. Fast Fourier Transform 13 X0= (x0+x2) + (x1+x3) X1= (x0−x2)−i·(x1−x3) X2= (x0+x2)−(x1+x3) X3= (x0−x2) + i·(x1−x3) (8) The FFT takes advantage of the periodic roots of unity values and the repetitive operations to generalize the grouping in Equation 8 for a sequence of size N. This algorithm itself is a multiplication of the input sequence by a sparse matrix (Heckbert (1995)). It achieves this by computing the butterfly operations at each stage and storing them to allow the reuse of the intermediate values in the next stage until the DFT is computed for every element. It is worth noting that for power of 2 sizes, the roots of unity are symmetrical and periodic over N , therefore, any twiddle factor ωk N is equal to ωk+N N and −ωk+N/2 N . For example, in figure Figure 2 we have the roots of unity for N=8, where ω0 8=−ω4 8and ω0 8=ω8 8. Figure 2: Roots of unity with N=8Source: Heckbert (1995) Section 2.3.1 and Section 2.3.2 focus on the Cooley–Tukey algorithm for sequences with power of 2 sizes. In Section 2.3.1 we describe how to generalize the FFT algorithm with decimation in time, and explain where this algorithm reduces the computational cost of the DFT. Moreover, in Section 2.3.2 we present another well-known variation of this algorithm with decimation in frequency. 2.3.1 Radix-2 Decimation-in-Time FFT The Radix-2 Decimation-in-Time FFT algorithm rearranges the computation of a DFT of size N into two DFTs of size N/2, one as a sum over the even indexed elements and the other as a sum over the odd indexed elements. 2.3. Fast Fourier Transform 14 Cooley and Tukey proved this possibility of dividing the DFT computation into two smaller DFTs by exploiting this division (Cooley and Tukey (1965)), as presented in Equation 10. Hence, it is hinted the recursive definition of this formulation on both DFT of size N/2. Xk= N−1 ∑ n=0 xn·ωkn N Xk= N/2−1 ∑ n=0 x2n·ωk(2n) N+ N/2−1 ∑ n=0 x2n+1·ωk(2n+1) N Xk= N/2−1 ∑ n=0 x2n·ωkn N/2 +ωk N N 2−1 ∑ n=0 x2n+1·ωkn N/2 (9) where ωN=ei2π N This formulation segments the full-sized DFT into two N/2 sized DFTs of the even and odd indexed elements where the latter is multiplied by a factor called twiddle ωk N . We can notice that in the N/2 sub-transform, the twiddle factors of ωk N/2 will be periodic for k=0, 1, ..., N−1 . This periodicity can also be noticed in Equation 8 within the sub-transforms of the even and odd elements. This algorithm is a Radix-2 Decimation-in-Time in the sense that elements are regrouped into 2 sub-transforms, and the decomposition reduces the time values to the frequency domain. The understanding of this algorithm comes from the recursive application of Equation 10, to compute the FFT for all the values of the input sequence. Figure 3 illustrates the composition of the sub-transforms, where the DFTs of size N/2 can be replaced with the same decomposition of even and odd sub-transforms. Figure 3: Radix-2 Decimation-in-Time FFT Source: Jones (2014) Effectively, the main building block in the FFT is the butterfly operation of 2 elements since it calculates intermediate values at each stage to compute the FFT of the whole sequence. This butterfly operation resembles the computation of a length-2 DFT but with a shifted element according to the sub-transform size 2.3. Fast Fourier Transform 15 This butterfly, without any reuse of the twiddle factor, corresponds to: A=a+b·ωk N B=a+b·ωk+N/2 N (10) where A=xkand B=xk+N/2 Yet, since the roots of unity for k and k+N/2 over N are symmetrical, as illustrated in Figure 2, we may re-express this butterfly in terms of the same twiddle factor ωk N , such that ωk+N/2 N will be equal to −ωk N (Jones (2014)). With this said, Cooley-Tukey butterfly (Chu and George (1999)) reuses the same twiddle factor, as illustrated in Figure 4. Figure 4: Cooley-Tukey butterfly The complexity work of the algorithm is distributed within the DIT approach which decomposes each DFT by 2 having log (N) stages Smith (2007). There are N complex multiplications needed for each stage of the DIT decomposition, therefore, the multiplication complexity for a N sized DFT is reduced from O(N2) to O(Nlog (N)). The splitting of the DFT into two smaller half-sized DFTs causes the original input sequence to require a special reordering to pass the even and odd numbers into, and when this algorithm is applied recursively, this reordering is always needed, so in the end, we need a special order for the input elements, fortunately, as noted by Thong (1981) this order corresponds to the bit reversed indices of the sequence, therefore, we need to apply a bit reverse permutation to the elements of the input sequence. The bit reversal of the input sequence corresponds to the permutation of swapping the position of the elements to its bit reversed index, as illustrated in Figure 5. The bit_reverse of an index depends directly on the indexing domain of the input sequence, therefore, it needs the size N , or more precisely the log (N) value, to use as a reference to reverse the bit order while maintaining the value within the sequence range. For example, for a sequence of size 16, we have some index with log 16 bits b1b2b3b4 , which corresponds to the bit reversed index as b4b3b2b1. For the DIT algorithm, we apply the bit reversal permutation to the input sequence to return a natural order result. 2.3. Fast Fourier Transform 16 Figure 5: Bit reverse permutation There are many implementations of the bit reversal, and since it is quite simple, any suitable version can be used in regard to this FFT algorithm since it is not the main bottleneck. An algorithm such as Algorithm 1 can be used for the bit_reverse function or any other efficient alternatives (Prado (2004)). Algorithm 1: Bit reverse Data: Integer i Result: Bit reversed integer i n←0 foreach i=0to log (N)−1do n←n<< 1; n←n|(i&1); i←i>> 1; end return i; In practice, Algorithm 2 demonstrates the aforesaid with an iterative DIT implementation for the forward FFT. To compute its inverse the algorithm is the same, we just need to use the twiddle factor with a positive exponent wm←exp(2πi/m). Although this algorithm is congruent with a code implementation, it is worth noting that the input sequence can either have real or complex numbers, since the arithmetic is the same for both domains the only thing that needs to be specialized is the operator overloading in the innermost loop. 2.5. Radix-4 instead of Radix-2 23 We’ve discussed multiple radix-2 approaches, however, we can explore a wide range of alternatives when we get into higher radices other than just radix-2. These FFT algorithms can use higher radix for better performance and even mixed radix Singleton (1969) for a wider range of input sequence sizes. Rearranging the Stockham algorithm to radix-4 upgrades the computation of the butterflies while reducing the number of stages which will be crucial in later sections. The radix-4 performs the work of two radix-2 iterations with less memory accesses (Bailey (1988)). Theoretically, radix-4 formulation can be twice as fast as a radix-2 Hussain et al. (2010) since it only takes half the stages with more complexity in the butterflies which are sometimes called dragonflies. Additionally, we can use less multiplications with better factorizations (Marti-Puig and Bolano (2009)). The radix-4 Stockham splits the FFT of size N into four sub-transforms of size N/4 each stage, therefore this algorithm only features log(N)/2 stages and requires the size to be power of 4 . Since this is a natural order algorithm the computation of the dragonfly includes the reordering of the elements in natural order at every stage, as illustrated in Figure 9. Figure 9: Illustration of a stage in radix-4 Stockham with each color representing a radix-4 butterfly computation The dragonfly for this algorithm is illustrated in Figure 10 and its computation involves 4 elements that compute 2stages of 2radix-2 butterflies. In Equation 15 is presented the forward radix-4 dragonfly operations and in Equation 16 its inverse. Figure 10: Radix-4 FFT butterfly structure Source: Marti-Puig and Bolano (2009) 2.5. Radix-4 instead of Radix-2 24 A= (a+c) + (b+d) B= ((a−c)−i(b−d)) ×ωp N C= ((a+c)−(b+d)) ×ω2p N D= ((a−c) + i(b−d)) ×ω3p N (15) A= (a+c) + (b+d) B= ((a−c) + i(b−d)) ×ω−p N C= ((a+c)−(b+d)) ×ω−2p N D= ((a−c)−i(b−d)) ×ω−3p N (16) With each stage, the dragonflies are calculated and the elements are reordered similarly to the radix-2 Stockham algorithm. In the end, we get the forward radix-4 Stockham algorithm Algorithm 5. The inverse algorithm is similar, but instead uses the inverse dragonfly as presented in Equation 16. 2.5. Radix-4 instead of Radix-2 25 Algorithm 5: Stockham Radix-4 Decimation-in-Frequency Forward FFT Data: Sequence pingpong0 with size Npower of 4 Result: Sequence out with size Nwith the DFT of the input foreach stage =0to log (N)/2 −1do n←1≪(((log (N)/2)−stage)∗2); s←1≪(stage ∗2); n0←0; n1←n/4; n2←n/2; n3←n1+n2; foreach i=0to N−1do p←idiv s; q←imod s; w1p=exp(−2∗π∗i∗p/n); w2p=w1p∗w1p; w3p=w1p∗w2p; if stage mod 2== 0then a←pingpong0[q+s∗(p+n0))]; b←pingpong0[q+s∗(p+n1)]; c←pingpong0[q+s∗(p+n2)]; d←pingpong0[q+s∗(p+n3)]; /*Perform dragonfly */ pingpong1[q+s∗(4∗p+0)] = a+c+b+d; pingpong1[q+s∗(4∗p+1)] = w1p∗((a−c)−(b−d)∗√−1); pingpong1[q+s∗(4∗p+2)] = w2p∗((a+c)−(b+d)); pingpong1[q+s∗(4∗p+3)] = w3p∗((a−c) + (b−d)∗√−1); else a←pingpong1[q+s∗(p+n0))]; b←pingpong1[q+s∗(p+n1)]; c←pingpong1[q+s∗(p+n2)]; d←pingpong1[q+s∗(p+n3)]; /*Perform dragonfly */ pingpong0[q+s∗(4∗p+0)] = a+c+b+d; pingpong0[q+s∗(4∗p+1)] = w1p∗((a−c)−(b−d)∗√−1); pingpong0[q+s∗(4∗p+2)] = w2p∗((a+c)−(b+d)); pingpong0[q+s∗(4∗p+3)] = w3p∗((a−c) + (b−d)∗√−1); end end end if log (N)mod 2== 0then return pingpong1; else return pingpong0; end 2.6. Two real inputs within one complex 26 2.6 TWO REAL INPUTS WITHIN ONE COMPLEX It is possible that in the case of an application there is a requirement to compute multiple real FFTs at once. Although we may simply invoke another FFT computation, we can optimize this step by computing multiple FFTs in the context of the same transform. We can take advantage of the complex input of our implementation to encode multiple real values to avoid extra explicit FFTs. Since the input has elements in a complex format, but the values are real, we only use the real component of the element, therefore, the output of the inverse also contains the complex with only the real component, as illustrated in Figure 11. Figure 11: FFT for an element without an imaginary part Based on this, we can reuse the imaginary part of the complex with a meaningful value other than 0 , as illustrated in Figure 12. Figure 12: FFT for an element with an imaginary part The result of the forward pass will be a mixed frequency value for the two original real values. Additionally, the frequency domain of both input sequences can be extracted from the mixed frequency value by performing some calculations, as demonstrated by Shima (2000). These calculations exploit the symmetry property of the frequency domain when the input sequence fills the real part of the complex and when it fills the imaginary part. With this said we may extract each FFT sequence from the mixed transform with the presented pack and unpack formulas in Equation 17 and Equation 18 respectively. When considering two real valued sequences x1 and x2 , whose transforms are X1 and X2 , and having y=x1+ix2, we can get the mixed frequencies in Yby packing the separated transforms as in Equation 17. Y(k) = X1(k) + X2(k)(17) When we have the mixed frequency components Y , we may unpack X1 and X2 as in Equation 18, where Yr and Yiare the real and imaginary parts of the complex-valued element. X1(k) = Y(k) + (Yr(N−k)−Yi(N−k)) 2 X2(k) = Y(k)−(Yr(N−k)−Yi(N−k)) 2 (18) 2.7. 2D Fourier Transform 27 2.7 2D FOURIER TRANSFORM Up until now, we only described transforming one-dimensional sequences, however, 2D FFTs are not that different from applying multiple 1D FFTs. Calculating the forward 2D FFT of a two-dimensional sequence produces the frequency domain result. For images, this frequency domain result corresponds to the change of pixel intensities in the original image. Calculating a 2D Fourier Transform requires a two-dimensional input sequence that the Forward 2D FFT converts numerical elements from real to complex domain and complex to the real domain on its inverse. When dealing with images, we might have multiple values per pixel element, therefore, we have the possibility of using it as a greyscale image if we derive the relative luminance via quantized RGB signals of the image (ITU (2002)), using only the values of one channel, or compute multiple FFTs for each channel values. With this said, when transforming the input in the forward FFT we will obtain a complex-valued output, and when applying the inverse we obtain the real-valued result. The 2D FFT is computed by performing single-dimension FFTs horizontally and vertically. In our case, we implemented this by first performing the 1D FFTs for every row and then every column, that is, 1 2D FFT of size N×N is equivalent to 2∗N 1D FFT of size N , as presented in Equation 19. This is called the row–column decomposition (Mermer et al. (2003)). Xr,c= N−1 ∑ m=0 N−1 ∑ n=0 xn,m·ωrn N·ωcn N(19) In later sections, we refer to this decomposition in the implementation as the horizontal and vertical pass, with the vertical pass being executed after the horizontal (Figure 13). Figure 13: High-level illustration of horizontal and vertical passes Additionally, the implementation of the 1D FFTs for every row and column is independent of the 2D FFT structure, hence, there is the freedom to use any of the discussed algorithms. In Chapter 3, we implement multiple FFT algorithms while reusing the same architecture of the horizontal and vertical passes. 3 IMPLEMENTATION IN GLSL The previous chapters introduced and expanded upon FFT algorithms to compute the DFT, which provides enough background to back up the GPU implementations presented in this chapter. As a result, in this chapter, we apply this background and implement these algorithms in GLSL compute shaders as a two-dimensional transform. The implementations were made using GLSL, a high-level shader language for graphics APIs such as OpenGL, and we used the compute pipeline from OpenGL to integrate an FFT implementation using compute shaders, which are general-purpose programmable shaders. As there are many aspects that can impact performance, the implementation was an iterative process that required research and testing to compose a more suitable solution. The next sections go into detail about the way every major iteration evolved into the next one and why there was the need to do it, starting with the Cooley-Tukey algorithm (Section 3.1) and then progressively improving to the Stockham algorithm (Section 3.2 and Section 3.3) 3.1 COOLEY-TUKEY The GPU implementation took as a starting point the Cooley-Tukey algorithm since it is the most popular one with a time complexity of O(Nlog (N)), and it is based on the iterative version adapted for parallel processors. Since this implementation runs multiple threads in parallel we separated the reads and writes for each processor into two different buffers in memory due lack of order between processors, therefore, the declaration of two complex pingpong buffers. layout (binding = 0, rg32f) uniform image2D pingpong0; layout (binding = 1, rg32f) uniform image2D pingpong1; Listing 3.1: Input buffer bindings The read/write control for these buffers can be achieved with a flag variable pingpong. 28 3.1. Cooley-Tukey 29 Initially, this algorithm will work with a pass-per-stage approach, where a kernel is dispatched every stage, and since there is a synchronization step at the end of the pass it is granted that all the work groups have finished off writing to the buffer when a pass ends. This case holds up for both the horizontal and vertical FFT steps. As a result, each kernel has the opportunity to work within every segment of the image, hence, the local threads can be dispatched with two dimensions. Each work group will have a total of 32 local threads, this number may vary depending on the GPU for optimal performance. An example dispatch group for this implementation could be (f f t_width/8, f f t_height/8) work groups, since 8 is the number of threads in the yaxis. By using GLSL there are some advantages to complex value operations since the addition and subtraction for vector types already have operator overloadings that work the same as in the complex domain. However, the multiplication works a bit differently so we need to provide an auxiliary function to abstract and support this operator. vec2 complex_mult(vec2 v0, vec2 v1) { return vec2(v0.x *v1.x - v0.y *v1.y, v0.x *v1.y + v0.y *v1.x); } Listing 3.2: Complex multiplication Due to the adoption of a different programming paradigm the FFT segment iteration loop doesn’t exist such as in Algorithm 2, instead, the process identifiers are used to fetch the index of the butterflies they’re gonna work on based on the work groups and threads dispatch setup, and this holds up for any implementation using compute shaders. int line = int(gl_GlobalInvocationID.x); int column = int(gl_GlobalInvocationID.y); Listing 3.3: Invocation indices Since this first approach is a dynamic implementation that invokes a pass-per-stage some stage control variables need to be fed into the shader in order to compute the correct butterfly index or control the butterfly process. uniform int pingpong; uniform int log_width; uniform int stage; uniform int fft_dir; Listing 3.4: Uniform control variables Effectively, we use these shader uniform input variables and obtain the actual index we are gonna use on the 1D FFT of the image. int group_size = 2 << stage; 3.1. Cooley-Tukey 30 int shift = 1 << stage; int idx = (line % shift) + group_size *(line / shift); Listing 3.5: FFT element index To calculate the twiddle factor we use Euler’s formula (Equation 2) and resort to the control variable fft_dir to flip the twiddle factor for the inverse if we want to reuse this shader. vec2 euler(float angle) { return vec2(cos(angle), sin(angle)); } void main() { // ... vec2 w = euler(fft_dir *2*(M_PI / group_size) *((idx % group_size) % shift)); // ... } Listing 3.6: Euler’s formula Now with the computed twiddle factor, we may proceed to compute and store the Cooley-Tukey FFT butterfly, and that’s where we need to control the reads and writes for each stage. Since the pingpong variable is toggled every pass invocation, it is used to choose with what texture we will use to look up the elements and which one to store into, and that is achieved with an if statement, just like it is used in Algorithm 4. The butterfly computation itself is simply calculated as a Cooley-Tukey DIT butterfly as presented in Section 2.3.1. if (pingpong == 0) { // Read a = imageLoad(pingpong0, ivec2(idx, column)).rg; b = imageLoad(pingpong0, ivec2(idx + shift, column)).rg; // Compute and store vec2 wb = complex_mult(w, b); vec2 raux = a + wb; imageStore(pingpong1, ivec2(idx, column), vec4(raux, 0, 0)); raux = a - wb; imageStore(pingpong1, ivec2(idx + shift, column), vec4(raux, 0, 0)); } else { // Read a = imageLoad(pingpong1, ivec2(idx, column)).rg; b = imageLoad(pingpong1, ivec2(idx + shift, column)).rg; // Compute and store 3.1. Cooley-Tukey 31 vec2 wb = complex_mult(w, b); vec2 raux = a + wb; imageStore(pingpong0, ivec2(idx, column), vec4(raux,0,0)); raux = a - wb; imageStore(pingpong0, ivec2(idx + shift, column), vec4(raux,0,0)); } Listing 3.7: Computation of the Cooley-Tukey butterfly Finally, there is only one step missing, the bit reversal of indices to have a natural order result. A bit_reverse function could be easily defined, However, there’s already a GLSL alternative which is bitfieldReverse together with bitfieldExtract (Kessenich et al.). int bit_reverse(int k) { uint br = bitfieldReverse(k); return int(bitfieldExtract(br, 32 - log_width, log_width)); } Listing 3.8: Auxiliary function that takes advantage of GLSL’s predefined utilities This auxiliary function will be conditionally used inside the branching if statement for alternating read/writes to be only applied on the first stage of transform both for the possible values of pingpong == 0 and pingpong == 1 , since the vertical pass might start reading on the first pingpong buffer. This is not the case for the horizontal pass, it is ensured that the first stage will always read from pingpong0 reinforcing that the bit reverse branching when pingpong == 0 is disposable. if (pingpong == 0) { if (stage == 0) { a = imageLoad(pingpong0, ivec2(bit_reverse(idx), column)).rg; b = imageLoad(pingpong0, ivec2(bit_reverse(idx + shift), column)).rg; } else { a = imageLoad(pingpong0, ivec2(idx, column)).rg; b = imageLoad(pingpong0, ivec2(idx + shift, column)).rg; } // ... Compute and store results } else { a = imageLoad(pingpong1, ivec2(idx, column)).rg; b = imageLoad(pingpong1, ivec2(idx + shift, column)).rg; // ... Compute and store results } Listing 3.9: Computation of the Cooley-Tukey butterfly with bit reversal 3.1. Cooley-Tukey 32 With all these steps aggregated, the shader for the horizontal pass of this FFT DIT Cooley-Tukey implementation is presented in Listing A.1, see Appendix A. For the vertical pass shader, there is, however, one extra multiplication by mult_factor in the butterfly results for the last stage when the pass is inverse. This corresponds to the normalization of the multidimensional transform similar to the normalization in the inverse DFT in Equation 5, this is demonstrated in Listing A.2. 3.1.1 All stages in one pass The previous implementation demonstrates a generic 2D FFT that may be reused for multiple FFT sizes, However, this comes at a cost of efficiency when it comes to the synchronization of the stage itself, moreover, it requires multiple uniform variables for control of the FFT that are transferred between CPU and GPU when there are updates in between stages. The ideal solution here is to port this implementation synchronization step to be kernel-wise and this detail changes how the code is structured and how it may be dispatched. There is a lot of literature on behalf of GPU compute execution synchronization (Stuart and Owens (2011)). In our case, we use execution barriers that synchronize all the threads within a work group (Kessenich et al.). The current grouping of threads doesn’t allow the use of these barrier synchronizations since the computation of one-dimensional FFT is distributed between multiple work groups, thus, using a call to barrier() inside the kernel wouldn’t fix the race conditions of several segments of the image. We could, however, change the setup of these work groups in such a way that each 1D FFT threads fit in one work group as illustrated in Figure 14. Figure 14: Difference in invocation spaces for size 256 FFT to allow barrier synchronization between local threads By restricting the invocation space to be only one dimensional we grant the possibility to use the barrier correctly but at the cost of resizability, the work group’s local size must now be restricted to half the size of the FFT. This restriction does not affect our results, where the tests were performed on images with a relatively small size. The GPU is extremely parallel hardware, therefore, it is better to dispatch smaller programs with more instances than overloading local threads in a work group to decrease its. 4.2. GLSL implementation results 39 •Power of 2 input sizes, such as 128, 256, 512, and 1024; •Complex to complex FFT, input and output buffer are complex-valued; • The benchmarks are average milliseconds of multiple executions, however, the first dispatch is not taken into account since takes extra time to set up things on the GPU. The results of the cuFFT out-of-place benchmarks are presented in Figure 16. 4.2 GLSL IMPLEMENTATION RESULTS The implementations discussed in Chapter 3 were studied and benchmarks were made to come to a conclusion about the advantages and disadvantages of using each one and how do they perform. With this in mind, we prepared an interactive test environment using Nau 3D engine (Ramires) and profiled it using an internal pass profiler. The benchmark results in this section were tested with the following hardware and software configurations: •CPU: Intel(R) Core(TM) i7-8750H @ 2.20GHz; •GPU: NVIDIA GeForce GTX 1050 Ti Max-Q; •NVIDIA driver: 511.65; •CUDA version: V11.6.124; •GLSL version: 4.60. In Section 3.1.1 we discussed how the implementation would benefit by having a unique pass that synchronized by stage instead of dispatching multiple stage passes. Accordingly, to prove this, we evaluated the difference between the Cooley-Tukey algorithm in the unique pass and in a stage-per-pass approach. Hence, in Figure 15 are presented the CPU and GPU results of the benchmarks for this test case. The GPU time corresponds to the total time spent in the GPU executing the compute shader. These benchmarks were tested in Nau3D, a generic 3D rendering engine, therefore, the CPU side of the stage-per-pass implementation takes into account the execution of a Lua script to update the state and control variables at the end of each stage. For this reason, the CPU time besides the preparation of the compute dispatch also includes some overhead for using this Lua script. 4.2. GLSL implementation results 40 Figure 15: CPU and GPU benchmarks of 2D forward FFT using stage-per-pass approach and unique pass for the Radix-2 Cooley-Tukey algorithm Sizes 128 256 512 1024 Unique pass CPU 0.06 0.052 0.049 0.06 GPU 0.05 0.15 0.598 2.698 stage-per-pass CPU 0.383 0.402 0.444 0.5 GPU 0.11 0.199 0.973 3.899 Table 1: CPU and GPU benchmarks of 2D forward FFT using stage-per-pass approach and unique pass for the Radix-2 Cooley-Tukey algorithm The results in Figure 15 and Table 1 demonstrate the approximate expected overhead of using the stageper-pass in an application. Since there is a need for the update of the following stages, the implementation synchronizes with the GPU immediately over each stage, however, the compute shader is reusable, therefore, the same shader module can be used for multiple sizes. In contrast, the CPU time of the unique pass kernel is mostly constant and this approach is optimized for its own size so most calculations are inlined by the GLSL compiler. The stage synchronization is kept inside the GPU until the kernel completes the execution. Since this type of synchronization is on the work group level, much fewer actors are involved in the synchronization step. In addition, these actors only correspond to the 1D FFT region of each row/column and not the entire 2D FFT. Due to the difference in the performance of these two approaches, the following comparisons will exclusively correspond to unique pass implementations of the investigated algorithms. As we can see in Figure 16 we plot the performance of the benchmarks of cuFFT and the unique pass algorithms in GLSL. These implementations are the single-pass approaches of the discussed algorithms in Chapter 3. Based on the findings of these benchmarks, the GLSL radix-2 implementation of the Stockham algorithm has an overall better performance comparing it to the Cooley-Tukey version. Additionally, it was also implemented with a smaller kernel. Consequently, this happens due to the removal of the data reordering process of the bit reversal done in the Cooley-Tukey version only in the first stage. Effectively, this change improves the results consistently within the test size ranges. 4.2. GLSL implementation results 41 Figure 16: Forward 2D FFT benchmarks in milliseconds of out-of-place cuFFT and unique pass algorithms in GLSL Sizes 128 256 512 1024 cuFFT 0.0359 0.0494 0.1335 0.5609 GLSL radix-2 Cooley-Tukey 0.073 0.257 1.032 2.646 GLSL radix-2 Stockham 0.049 0.135 0.545 2.341 GLSL radix-4 Stockham 0.042 0.087 0.389 1.363 Table 2: Benchmarks of the forward 2D FFT benchmarks in milliseconds of out-of-place cuFFT and unique pass algorithms in GLSL Undoubtedly, the results that come closer to the cuFFT are the ones of the radix-4 Stockham. It drastically improved the performance halfway close to the cuFFT, as presented in Figure 16 and Table 2. Yet, the complexity of the shader code was sharply increased, not only due to the higher radix alternative but also due to the additional support step in the last stage for the power of 2sizes. By choosing a radix-4 approach the number of stages reduces to half but with a lot more complex operations per dragonfly on each stage. Since this dragonfly represents the operations required in an equivalent radix-2 Stockham factorization. Although the work complexity remains the same, there are much fewer barrier synchronization events and radix-4 iterations perform the work of two radix-2 iterations with only one memory access. For each stage there’s a synchronization barrier in each local thread inside a work group, so fewer stages means fewer sync points, hence, compensating for the 70% kernel size increase of the radix-4 version. In the context of the same need for multiple FFTs, we reach a point where we need to batch executions together for cuFFT. While this can also be achieved using our GLSL implementations, we can take a step further on the usage of the same pass for double the FFTs. Accordingly, we adapted the radix-4 Stockham implementation to use vec4 instead of vec2 for the elements of the input and output buffers. Code integrity is conserved since changing to vec4 does not change operators in the compute kernel. 4.2. GLSL implementation results 42 To compare with an equivalent, the computation in cuFFT used a special setup of a 2D FFT plan to compute with a batch size of 2 . The cuFFT framework doesn’t have an explicit setup for multiple transforms in the same kernel, hence, we use batched execution for comparison. Figure 17: Forward 2D FFT benchmarks of double the transforms in milliseconds of out-of-place cuFFT and radix-4 Stockham GLSL In table Table 3, the benchmarks of Figure 17 are presented together with the results of calculating a single transform in one step. To analyze the potential gain that the implementation with twice the transforms gives, we also present the percentage of the performance gain per individual FFT. Sizes 128 256 512 1024 cuFFT Single FFT 0.035 0.049 0.133 0.560 Double FFT 0.038 0.098 0.23 1.032 Gain per single FFT 45.71% 0.0% 13.53% 7.85% GLSL radix-4 Stockham Single FFT 0.042 0.087 0.389 1.363 Double FFT 0.044 0.132 0.704 2.702 Gain per single FFT 47.61% 24.13% 9.51% 0.88% Table 3: Benchmarks of the Forward 2D FFT for cuFFT and GLSL radix-4 Stockham for single and double transforms in the same pass, together with the percentage of gain over individual transforms. As a result, we can clearly notice a larger gain for size 128 inputs, both for cuFFT and for GLSL. From there, the gain for the rest of the tested sizes is notably smaller, the most relevant gain being that of inputs of size 256 for GLSL since cuFFT did not show any difference to twice the time of a single FFT for this input size. Finally, for input sizes 512 and 1024 the cuFFT batched execution demonstrated a valuable gain while the GLSL one continued to decline. Overall, we can conclude that it is worth including multiple FFTs in the same pass for the tested sizes when multiple independent transforms are necessary. 4.3. Case of study 43 With the presented results, we can clearly notice better performance on the GPU regarding the implementation in GLSL of radix-4 Stockham with a unique pass approach that reduces the number of synchronization stages and calculates per dragonfly the equivalent 2butterflies of two stages of radix-2 Stockham. 4.3 C A S E O F S T U DY Based on the findings of Section 4.2 we measured how the GLSL implementations behaved, by analyzing the performance of a test application that converted a 2D image to the frequency domain and then reversed it to its original look. Yet, with the goal of highlighting the importance of the performance increase of these algorithms, in this section, we provide an overview of the impact of the implementation within a more demanding scenario by using an ocean rendering technique demo that heavily relies on the usage of FFT (Figure 18). Figure 18: Tensendorf waves in Nau3D Engine In this section we brief the Tensendorf waves demo in Section 4.3.1 where we describe in which way FFTs are relevant for the implementation of this rendering technique, and how we improved an existing implementation for Nau3D that used a pass-per-stage Cooley-Tukey implementation. After that we present results and how the FFT implementation improves the demo performance and by how much in Section 4.3.2. 4.3.1 Tensendorf waves The rendering of the oceans demo we used as a starting point was a real-time implementation in Nau3D of the popular article published by Tessendorf et al. (2001). In this demo there are two main stages, the generation of the height map and the actual rendering, the FFTs come into place in the generation of the height map since we need to generate it and the additional vectors used for shading. In total, there are 4 2D inverse FFTs computed for each frame, which translates to 8∗FFT_SIZE 1D FFTs in total. 4.3. Case of study 44 Regarding the results in Section 4.2, we first changed the pipeline from a pass-per-stage radix-2 Cooley-Tukey algorithm to implement radix-2 and radix-4 Stockham with synchronization within the kernel for the horizontal and vertical passes. Each pass computes multiple FFTs at a time and we take advantage of the same kernel to compute all required FFTs at the same time. Additionally, it is worth noticing that the inverse FFTs produce from the frequencies components, two usable real values. Finally, we used vec4 as described in Section 4.2 to compute multiple FFTs in the context of the SIMT instructions produced for the GPU. Although the FFTs take a big role in this demo, it also renders the ocean waves that have around 2 million vertices, hence, the performance does not only depend on the FFTs computation. For this demo, we tested its performance with sizes 512 and 1024 for the ping pong buffers. We used size 512 to test a balance of performance and wave quality and 1024 to test for better wave quality as suggested by Tessendorf et al. (2001). Similarly to Section 4.2, the radix-4 Stockham implementation for the size 512 integrates a final stage with radix-2 butterflies to be able to support the computation of this size. 4.3.2 Results Figure 19: Time spent in the CPU and GPU for the size 512 FFT horizontal and vertical passes Sizes 512 1024 Original CPU 1.805 1.946 GPU 3.241 14.086 Radix-2 Stockham CPU 0.191 0.203 GPU 2.457 13.48 Radix-4 Stockham CPU 0.249 0.2 GPU 1.632 7.621 Table 4: Time spent in the CPU and GPU for the size 512 FFT horizontal and vertical passes 4.3. Case of study 45 In Figure 19 and Table 4 we can note the performance difference the radix-2 Stockham gives when performing the horizontal and vertical pass comparing it to the radix-2 Cooley-Tukey with a pass-per-stage. Not only is there an improvement in time and resources used by eliminating CPU steps, but there is also better performance using unique pass versions. When running the application, the performance increase is noticeable and the frame rate is improved by up to 20%. Testing the demo for higher quality waves the radix-4 Stockham performance stands out in Figure 19, as predicted. When running the application, the difference in performance is noticeable, with the frame rate using Stockham radix-2 being improved by up to 20%, while the radix-4 delivers a frame rate of up to 60%. These results proved the need to implement adequate algorithms to ensure the quality and performance of the target application. 5 CONCLUSIONS AND FUTURE WORK In this thesis, we provided background on the Discrete Fourier Transform, with the goal of introducing the DIT and DIF Cooley-Tukey algorithms. From there, we described how the Stockham algorithm removes the need for the bit reversal permutation and what is the radix-4 of this algorithm. Moreover, we also discussed how a complex-valued FFT can be used to compute 2 real-valued input sequences without any additional cost, and how this transformation can be applied to the spatial domain of a 2D image. After presenting the Fast Fourier Transform and its algorithms, we processed implementing the Cooley-Tukey algorithm with a stage-per-pass approach followed by unique pass implementations of the radix-2 Cooley-Tukey and Stockham algorithm. Additionally, we also implemented the radix-4 Stockham algorithm, however, with the intent of having this algorithm support sizes of 2 instead of just power of 4 , we altered the final stage of the code to conditionally implement a final stage of radix-2 Stockham. Finally, to analyze these algorithms, we conducted a comparison of the performance of the implementations with benchmarks in a test environment and also measured implementation details such as the advantages and disadvantages of using stage-per-pass versus unique-pass and double the FFTs in GLSL by reusing the same vec4 operations. From the compared algorithms, the best results were those of radix-4 Stockham, which we used as a model to test the possible benefits in the case study. The case study demonstrated positive results, by applying the previously discussed knowledge in the real-time rendering technique with radix-4 Stockham multiple FFTs. 5.1 RESULTS Based on the findings of Chapter 4, we were able to conclude several aspects of the tests carried out. When it comes to the tested unique-pass approach, it was demonstrated to be faster than the stage-per-pass, yet, we concluded multiple advantages and disadvantages: Advantages: •Synchronization only on the GPU, and not with the CPU until the 2D FFT is fully calculated; • The work synchronizes only on the 1D FFT, and not on the whole 2D FFT, hence, it finishes the work earlier than the stage-per-pass approach; 46 5.2. Future work 47 • Unique-pass shader has more constant expressions and information compared to the stage-per-pass shader, thus it is more prone to compiler optimizations. Disadvantages: •Higher kernel size; • The implementation is specialized for a specific size, therefore, not as flexible to reuse for other sizes. While the stage-per-pass implementation is reused for multiple power of 2 sizes, the unique-pass strategies require changing the provided macro constants and recompiling as an additional shader; • Mixed radix stages do not have special invocation spaces. This was the case of the radix-4 Stockham for power of 2 sizes, we were required to use duplicate the amount of radix-2 Stockham stage performed by each thread since there were only FFT_SIZE/4 threads. Furthermore, when it comes to comparing implementations, it was found that radix-4 Stockham obtained the best results, up to 2× faster than radix-2 Cooley-Tukey. This was mainly due to the elimination of the bit reversal permutation and the reduction in the number of memory accesses by half due to the computation of 2 radix-2 butterflies of 2radix-2 stages at once. When the acquired knowledge was applied to the case study, as expected, there was an increase in performance similar to that of the tests carried out. This case study applied strategies that were described and tested, such as the use of double FFTs by using the vec4 operations and the multiple real values from the FFT frequencies. The results improved the frames per second of the application by up to 60%. 5.2 F U T U R E W O R K With the obtained results, other questions arise, which may lead to future studies to expand the work carried out. Although the radix-4 reduced the memory accesses by half, these memory loads still weigh on the performance of the compute shaders, other strategies for this problem could be studied to reduce the overhead this introduces. Since the results of the radix-4 Stockham were the fastest when using radix higher than 4, what is the most suitable configuration for dispatching a higher radix implementation while still supporting power of 2 sizes? With the size restriction of the invocation space until the kernel is done executing, the unique-pass approach needs to duplicate the number of butterflies for each lower radix stage of the implementation. Therefore, we could study the implementation of these higher radix algorithms to verify if the usage of stage passes for the lower radix stages would be more beneficial than strictly having a unique-pass approach. A study could be carried out to explore these conditions for multiple GPUs, and come to a conclusion on the best configuration on average for the tested hardware. BIBLIOGRAPHY Michael D Adams. Signals and Systems (Edition 3.0). Michael Adams, 2020. David H Bailey. A high-performance fft algorithm for vector supercomputers. The International Journal of Supercomputing Applications, 2(1):82–87, 1988. Gabriel Bengtsson. Development of stockham fast fourier transform using data-centric parallel programming, 2020. Leo Bluestein. A linear filtering approach to the computation of discrete fourier transform. IEEE Transactions on Audio and Electroacoustics, 18(4):451–455, 1970. E Oran Brigham. The fast Fourier transform and its applications. Prentice-Hall, Inc., 1988. Eleanor Chu and Alan George. Inside the FFT black box: serial and parallel fast Fourier transform algorithms. CRC press, 1999. James W Cooley and John W Tukey. An algorithm for the machine calculation of complex fourier series. Mathematics of computation, 19(90):297–301, 1965. Fynn-Jorin Flügge. Realtime gpgpu fft ocean water simulation. Technical report, 2017. Naga K Govindaraju, Brandon Lloyd, Yuri Dotsenko, Burton Smith, and John Manferdelli. High performance discrete fourier transforms on graphics processors. In SC’08: Proceedings of the 2008 ACM/IEEE conference on Supercomputing, pages 1–12. Ieee, 2008. Paul Heckbert. Fourier transforms and the fast fourier transform (fft) algorithm. Computer Graphics, 2:15–463, 1995. Roger A Horn and Charles R Johnson. Matrix analysis. Cambridge university press, 2012. Waqar Hussain, Fabio Garzia, and Jari Nurmi. Evaluation of radix-2 and radix-4 fft processing on a reconfigurable platform. In 13th IEEE Symposium on Design and Diagnostics of Electronic Circuits and Systems, pages 249–254. IEEE, 2010. ITURBT ITU. Parameter values for the hdtv standards for production and international programme exchange. Recommendation ITU-R BT, pages 709–5, 2002. Yan-Bin Jia. Polynomial multiplication and fast fourier transform. Com S, 477:577, 2020. Douglas L Jones. Digital signal processing: A user’s guide. 2014. 48 55 a = imageLoad(pingpong0, ivec2(line, idx)).rg; b = imageLoad(pingpong0, ivec2(line, idx + shift)).rg; } vec2 wb = complex_mult(w, b); vec2 raux = (a + wb) *mult_factor; imageStore(pingpong1, ivec2(line, idx), vec4(raux,0,0)); raux = (a - wb) *mult_factor; imageStore(pingpong1, ivec2(line, idx + shift), vec4(raux,0,0)); } else { if (stage == 0) { a = imageLoad(pingpong1, ivec2(line, bit_reverse(idx))).rg; b = imageLoad(pingpong1, ivec2(line, bit_reverse(idx + shift))).rg; } else { a = imageLoad(pingpong1, ivec2(line, idx)).rg; b = imageLoad(pingpong1, ivec2(line, idx + shift)).rg; } vec2 wb = complex_mult(w, b); vec2 raux = (a + wb) *mult_factor; imageStore(pingpong0, ivec2(line, idx), vec4(raux,0,0)); raux = (a - wb) *mult_factor; imageStore(pingpong0, ivec2(line, idx + shift), vec4(raux,0,0)); } } Listing A.2: FFT Radix-2 Cooley-Tukey Vertical stage pass, see Section 3.1 #version 440 #define M_PI 3.1415926535897932384626433832795 #define FFT_SIZE 256 #define LOG_SIZE 8 layout (local_size_x = FFT_SIZE/2, local_size_y = 1) in; layout (binding = 0, rg32f) uniform image2D pingpong0; layout (binding = 1, rg32f) uniform image2D pingpong1; uniform int fft_dir; vec2 complex_mult(vec2 v0, vec2 v1) { return vec2(v0.x *v1.x - v0.y *v1.y, 56 v0.x *v1.y + v0.y *v1.x); } int bit_reverse(int k) { uint br = bitfieldReverse(k); return int(bitfieldExtract(br, 32 - LOG_SIZE, LOG_SIZE)); } vec2 euler(float angle) { return vec2(cos(angle), sin(angle)); } void main() { int line = int(gl_GlobalInvocationID.x); int column = int(gl_WorkGroupID.y); int pingpong = 0; for(int stage = 0; stage < LOG_SIZE; ++stage) { int group_size = 2 << stage; int shift = 1 << stage; vec2 a, b; int idx = (line % shift) + group_size *(line / shift); vec2 w = euler(fft_dir *2*(M_PI / group_size) *((idx % group_size) % shift)); // alternate between textures if (pingpong == 0) { if (stage == 0) { a = imageLoad(pingpong0, ivec2(bit_reverse(idx), column)).rg; b = imageLoad(pingpong0, ivec2(bit_reverse(idx + shift), column)) .rg; } else { a = imageLoad(pingpong0, ivec2(idx, column)).rg; b = imageLoad(pingpong0, ivec2(idx + shift, column)).rg; } vec2 wb = complex_mult(w, b); vec2 raux = a + wb; imageStore(pingpong1, ivec2(idx, column), vec4(raux, 0, 0)); raux = a - wb; imageStore(pingpong1, ivec2(idx + shift, column), vec4(raux, 0, 0)); } else { a = imageLoad(pingpong1, ivec2(idx, column)).rg; b = imageLoad(pingpong1, ivec2(idx + shift, column)).rg; 57 vec2 wb = complex_mult(w, b); vec2 raux = a + wb; imageStore(pingpong0, ivec2(idx, column), vec4(raux,0,0)); raux = a - wb; imageStore(pingpong0, ivec2(idx + shift, column), vec4(raux,0,0)); } pingpong = (pingpong + 1) % 2; barrier(); } } Listing A.3: FFT Radix-2 Cooley-Tukey Horizontal unique pass, see Section 3.1.1 #version 440 #define M_PI 3.1415926535897932384626433832795 #define FFT_SIZE 256 #define LOG_SIZE 8 layout (local_size_x = FFT_SIZE/2, local_size_y = 1) in; layout (binding = 0, rg32f) uniform image2D pingpong0; layout (binding = 1, rg32f) uniform image2D pingpong1; uniform int fft_dir; vec2 complex_mult(vec2 v0, vec2 v1) { return vec2(v0.x *v1.x - v0.y *v1.y, v0.x *v1.y + v0.y *v1.x); } int bit_reverse(int k) { uint br = bitfieldReverse(k); return int(bitfieldExtract(br, 32 - LOG_SIZE, LOG_SIZE)); } vec2 euler(float angle) { return vec2(cos(angle), sin(angle)); } void main() { int line = int(gl_WorkGroupID.y); int column = int(gl_GlobalInvocationID.x); int pingpong = LOG_SIZE % 2; 58 for(int stage = 0; stage < LOG_SIZE; ++stage) { int group_size = 2 << stage; int shift = 1 << stage; vec2 a, b; int idx = (column % shift) + group_size *(column / shift); vec2 w = euler(fft_dir *2*(M_PI / group_size) *((idx % group_size) % shift)); float mult_factor = 1.0; if ((stage == LOG_SIZE - 1) && fft_dir == 1) { mult_factor = 1.0 / (FFT_SIZE*FFT_SIZE); } if (pingpong == 0) { if (stage == 0) { a = imageLoad(pingpong0, ivec2(line, bit_reverse(idx))).rg; b = imageLoad(pingpong0, ivec2(line, bit_reverse(idx + shift))). rg; } else { a = imageLoad(pingpong0, ivec2(line, idx)).rg; b = imageLoad(pingpong0, ivec2(line, idx + shift)).rg; } vec2 wb = complex_mult(w, b); vec2 raux = (a + wb) *mult_factor; imageStore(pingpong1, ivec2(line, idx), vec4(raux,0,0)); raux = (a - wb) *mult_factor; imageStore(pingpong1, ivec2(line, idx + shift), vec4(raux,0,0)); } else { if (stage == 0) { a = imageLoad(pingpong1, ivec2(line, bit_reverse(idx))).rg; b = imageLoad(pingpong1, ivec2(line, bit_reverse(idx + shift))). rg; } else { a = imageLoad(pingpong1, ivec2(line, idx)).rg; b = imageLoad(pingpong1, ivec2(line, idx + shift)).rg; } vec2 wb = complex_mult(w, b); vec2 raux = (a + wb) *mult_factor; imageStore(pingpong0, ivec2(line, idx), vec4(raux,0,0)); raux = (a - wb) *mult_factor; 59 imageStore(pingpong0, ivec2(line, idx + shift), vec4(raux,0,0)); } pingpong = ((pingpong + 1) % 2); barrier(); } } Listing A.4: FFT Radix-2 Cooley-Tukey Vertical unique pass, see Section 3.1.1 #version 440 #define M_PI 3.1415926535897932384626433832795 #define FFT_SIZE 256 #define LOG_SIZE 8 layout (local_size_x = (FFT_SIZE/2)/NUM_BUTTERFLIES, local_size_y = 1) in; layout (binding = 0, rg32f) uniform image2D pingpong0; layout (binding = 1, rg32f) uniform image2D pingpong1; uniform int fft_dir; vec2 complex_mult(vec2 v0, vec2 v1) { return vec2(v0.x *v1.x - v0.y *v1.y, v0.x *v1.y + v0.y *v1.x); } vec2 euler(float angle) { return vec2(cos(angle), sin(angle)); } void main() { int line = int(gl_GlobalInvocationID.x); int column = int(gl_WorkGroupID.y); int pingpong = 0; for(int stage = 0; stage < LOG_SIZE; ++stage) { int n = 1 << (LOG_SIZE - stage); int m=n>>1; int s = 1 << stage; int p = line / s; int q = line % s; vec2 wp = euler(fft_dir *2*(M_PI / n) *p); if(pingpong == 0) { 60 vec2 a = imageLoad(pingpong0, ivec2(q + s*(p + 0), column)).rg; vec2 b = imageLoad(pingpong0, ivec2(q + s*(p + m), column)).rg; vec2 res = (a + b); imageStore(pingpong1, ivec2(q + s*(2*p + 0), column), vec4(res,0,0)); res = complex_mult(wp,(a - b)); imageStore(pingpong1, ivec2(q + s*(2*p + 1), column), vec4(res,0,0)); } else { vec2 a = imageLoad(pingpong1, ivec2(q + s*(p + 0), column)).rg; vec2 b = imageLoad(pingpong1, ivec2(q + s*(p + m), column)).rg; vec2 res = (a + b); imageStore(pingpong0, ivec2(q + s*(2*p + 0), column), vec4(res,0,0)); res = complex_mult(wp,(a - b)); imageStore(pingpong0, ivec2(q + s*(2*p + 1), column), vec4(res,0,0)); } pingpong = (pingpong + 1) % 2; barrier(); } } Listing A.5: FFT Radix-2 Stockham Horizontal unique pass, see Section 3.2 #version 440 #define M_PI 3.1415926535897932384626433832795 #define FFT_SIZE 256 #define LOG_SIZE 8 layout (local_size_x = FFT_SIZE/2, local_size_y = 1) in; layout (binding = 0, rg32f) uniform image2D pingpong0; layout (binding = 1, rg32f) uniform image2D pingpong1; uniform int fft_dir; vec2 complex_mult(vec2 v0, vec2 v1) { return vec2(v0.x *v1.x - v0.y *v1.y, v0.x *v1.y + v0.y *v1.x); } vec2 euler(float angle) { return vec2(cos(angle), sin(angle)); } 61 void main() { int line = int(gl_WorkGroupID.y); int column = int(gl_GlobalInvocationID.x); int pingpong = LOG_SIZE % 2; for(int stage = 0; stage < LOG_SIZE; ++stage) { int n = 1 << (LOG_SIZE - stage); int m=n>>1; int s = 1 << stage; float mult_factor = 1.0; if ((stage == LOG_SIZE-1) && fft_dir == 1) { mult_factor = 1.0 / (FFT_SIZE*FFT_SIZE) ; } int p = column / s; int q = column % s; vec2 wp = euler(fft_dir *2*(M_PI / n) *p); if(pingpong == 0) { vec2 a = imageLoad(pingpong0, ivec2(line, q + s*(p + 0))).rg; vec2 b = imageLoad(pingpong0, ivec2(line, q + s*(p + m))).rg; vec2 res = (a + b) *mult_factor; imageStore(pingpong1, ivec2(line, q + s*(2*p + 0)), vec4(res,0,0)); res = complex_mult(wp,(a - b)) *mult_factor; imageStore(pingpong1, ivec2(line, q + s*(2*p + 1)), vec4(res,0,0)); } else { vec2 a = imageLoad(pingpong1, ivec2(line, q + s*(p + 0))).rg; vec2 b = imageLoad(pingpong1, ivec2(line, q + s*(p + m))).rg; vec2 res = (a + b) *mult_factor; imageStore(pingpong0, ivec2(line, q + s*(2*p + 0)), vec4(res,0,0)); res = complex_mult(wp,(a - b)) *mult_factor; imageStore(pingpong0, ivec2(line, q + s*(2*p + 1)), vec4(res,0,0)); } pingpong = (pingpong + 1) % 2; barrier(); } } Listing A.6: FFT Radix-2 Stockham Vertical unique pass, see Section 3.2 #version 440 62 #define M_PI 3.1415926535897932384626433832795 #define FFT_SIZE 256 #define LOG_SIZE 8 // log2(FFT_SIZE) #define HALF_LOG_SIZE 4 // log2(FFT_SIZE) / 2 layout (local_size_x = FFT_SIZE/4, local_size_y = 1) in; layout (binding = 0, rg32f) uniform image2D pingpong0; layout (binding = 1, rg32f) uniform image2D pingpong1; uniform int fft_dir; vec2 complex_mult(vec2 v0, vec2 v1) { return vec2(v0.x *v1.x - v0.y *v1.y, v0.x *v1.y + v0.y *v1.x); } vec2 euler(float angle) { return vec2(cos(angle), sin(angle)); } void main() { int line = int(gl_GlobalInvocationID.x); int column = int(gl_WorkGroupID.y); int pingpong = 0; for(int stage = 0; stage < HALF_LOG_SIZE; ++stage) { int n = 1 << (HALF_LOG_SIZE - stage)*2; int s = 1 << stage*2; int n0 = 0; int n1 = n/4; int n2 = n/2; int n3 = n1 + n2; int p = line / s; int q = line % s; vec2 w1p = euler(2*(M_PI / n) *p*fft_dir); vec2 w2p = complex_mult(w1p,w1p); vec2 w3p = complex_mult(w1p,w2p); if(pingpong == 0) { vec2 a = imageLoad(pingpong0, ivec2(q + s*(p + n0), column)).rg; vec2 b = imageLoad(pingpong0, ivec2(q + s*(p + n1), column)).rg; vec2 c = imageLoad(pingpong0, ivec2(q + s*(p + n2), column)).rg; vec2 d = imageLoad(pingpong0, ivec2(q + s*(p + n3), column)).rg; 63 vec2 apc = a + c; vec2 amc = a - c; vec2 bpd = b + d; vec2 jbmd = complex_mult(vec2(0,1), b - d); imageStore(pingpong1, ivec2(q + s*(4*p + 0), column), vec4(apc + bpd, 0,0)); imageStore(pingpong1, ivec2(q + s*(4*p + 1), column), vec4( complex_mult(w1p, amc + jbmd*fft_dir), 0,0)); imageStore(pingpong1, ivec2(q + s*(4*p + 2), column), vec4( complex_mult(w2p, apc - bpd ), 0,0)); imageStore(pingpong1, ivec2(q + s*(4*p + 3), column), vec4( complex_mult(w3p, amc - jbmd*fft_dir), 0,0)); } else { vec2 a = imageLoad(pingpong1, ivec2(q + s*(p + n0), column)).rg; vec2 b = imageLoad(pingpong1, ivec2(q + s*(p + n1), column)).rg; vec2 c = imageLoad(pingpong1, ivec2(q + s*(p + n2), column)).rg; vec2 d = imageLoad(pingpong1, ivec2(q + s*(p + n3), column)).rg; vec2 apc = a + c; vec2 amc = a - c; vec2 bpd = b + d; vec2 jbmd = complex_mult(vec2(0,1), b - d); imageStore(pingpong0, ivec2(q + s*(4*p + 0), column), vec4(apc + bpd, 0,0)); imageStore(pingpong0, ivec2(q + s*(4*p + 1), column), vec4( complex_mult(w1p, amc + jbmd*fft_dir), 0,0)); imageStore(pingpong0, ivec2(q + s*(4*p + 2), column), vec4( complex_mult(w2p, apc - bpd ), 0,0)); imageStore(pingpong0, ivec2(q + s*(4*p + 3), column), vec4( complex_mult(w3p, amc - jbmd*fft_dir), 0,0)); } pingpong = (pingpong + 1) % 2; barrier(); } } Listing A.7: FFT Radix-4 Stockham Horizontal unique pass, see Section 3.3 #version 440 #define M_PI 3.1415926535897932384626433832795 64 #define FFT_SIZE 256 #define LOG_SIZE 8 // log2(FFT_SIZE) #define HALF_LOG_SIZE 4 // log2(FFT_SIZE) / 2 layout (local_size_x = (FFT_SIZE/4)/NUM_BUTTERFLIES, local_size_y = 1) in; layout (binding = 0, rg32f) uniform image2D pingpong0; layout (binding = 1, rg32f) uniform image2D pingpong1; uniform int fft_dir; vec2 complex_mult(vec2 v0, vec2 v1) { return vec2(v0.x *v1.x - v0.y *v1.y, v0.x *v1.y + v0.y *v1.x); } vec2 euler(float angle) { return vec2(cos(angle), sin(angle)); } void main() { int line = int(gl_WorkGroupID.y); int column = int(gl_GlobalInvocationID.x); int pingpong = HALF_LOG_SIZE % 2; for(int stage = 0; stage < HALF_LOG_SIZE; ++stage) { int group_size = 2 << stage; int shift = 1 << stage; int n = 1 << ((HALF_LOG_SIZE - stage)*2); int s = 1 << (stage*2); int n0 = 0; int n1 = n/4; int n2 = n/2; int n3 = n1 + n2; float mult_factor = 1.0; if((stage == HALF_LOG_SIZE - 1) && fft_dir == 1) { mult_factor = 1.0 / (FFT_SIZE*FFT_SIZE) ; } int p = column / s; int q = column % s; vec2 w1p = euler(2*(M_PI / n) *p*fft_dir); vec2 w2p = complex_mult(w1p,w1p);