Computational Complexity of Radix-2, Radix-4 and Bluestein Algorithms Implementation of the Discrete Fourier Transform (DFT)

1. Introduction

The Discrete Fourier Transform (DFT) is a cornerstone of modern signal processing, enabling the transition between time-domain and frequency-domain representations of discrete signals. While the DFT is computationally intensive for large data sets, the development of Fast Fourier Transform (FFT) algorithms has revolutionized this field by significantly reducing computational complexity. Among the various FFT algorithms [1], Radix-2, Radix-4, and Bluestein [2] stand out for their unique advantages and applicability to different scenarios [3].

This paper investigates the computational complexity of these three prominent FFT algorithms, which have been implemented and analyzed using MATLAB [4]. Radix-2, known for its simplicity and efficiency for power-of-two data lengths, remains a widely adopted algorithm in many practical applications [3]. On the other hand, Radix-4 enhances this efficiency by processing four data points at a time, reducing the number of stages and multiplicative operations. However, both algorithms are constrained by their dependence on specific sequence lengths.

The Bluestein algorithm, often referred to as the chirp-z transform method, addresses this limitation by enabling FFT computations for arbitrary sequence lengths, including prime numbers. By reformulating the DFT as a convolution problem, the Bluestein algorithm leverages FFTs to compute convolutions efficiently. While this approach introduces additional computational overhead, it eliminates the need for zero-padding, which can distort frequency-domain results in certain applications [5].

To minimize computational complexity and assess the practical performance of these algorithms, MATLAB implementations were developed with optimization strategies to reduce the number of multiplications during runtime. These implementations facilitated a comparative analysis of the three algorithms, focusing on their multiplication and addition requirements across various data lengths.

This study aims to provide a comprehensive understanding of the trade-offs involved in selecting an FFT algorithm for specific applications. By evaluating the strengths and limitations of Radix-2, Radix-4, and Bluestein algorithms [2], this research offers valuable insights into their suitability for diverse signal processing tasks. Additionally, the findings contribute to the ongoing development and optimization of FFT algorithms, ensuring their continued relevance in addressing the computational demands of modern signal processing applications.

2. Complexity of Radix-2 DIT/DIF Algorithm

A MATLAB program has been written and is in Appendix 1, to compute the FFT of a set of discrete data: $\left\{x_n\right\},n=1,2,3,\dots ,N$ where $N=2^n$ and n is an integer. A modification of the supplied Radix-2 program was made in order to reduce the number of multiplications involved and is highlighted in the attached Radix-2 program [3,6].

The final loop in the program made use of the property ${\displaystyle e^{-j\frac{2\pi (i-1)}{N}}=-e^{-j\frac{2\pi \left(\frac{N}{2}i+1\right)}{N}}}$, which allowed the number of multiplications in the final loop to be halved. The program uses a recursive procedure which calls up the FFT of the odd and even set of data points [7].

In the special case when $N=2$, there are no multiplications to evaluate and only 2 additions. In the general case when $N=2^n$, the recursive procedure calls up the Fourier transform of two sets of data, each of length $N/2=2^{n-1}$. The final loop in the procedure contains $N/2$ multiplications and N additions. Using $N=2^n$, $n=1,2,3,\dots$, let $M\left(n\right)$ represent the number of computed multiplications and $A\left(n\right)$ represent the number of computed additions. It can be deduced that $M\left(n\right)$ and $A\left(n\right)$ satisfy the difference equations:

$$M(n+1)=M\left(n\right)+2^n,n=1,2,3,4,\dots \mathrm{where}M\left(1\right)=0$$

(1)

$$A(n+1)=A\left(n\right)+2^n,n=1,2,3,4,\dots \mathrm{where}A\left(1\right)=2$$

(2)

The solutions to these difference equations can be obtained by taking a z-transform of the equation. It can be shown that the formulas for $M\left(N\right)$ and $A\left(N\right)$ [7,8] are given by:

$$M\left(n\right)={\displaystyle \frac{N}{2}}[{log}_{2}N-1]$$

(3)

$$A\left(n\right)=N{log}_{2}N$$

(4)

The values for $M\left(N\right)$ and $A\left(N\right)$ have been tabulated as shown in the Table 1 below:

This data is consistent with Cooley and Tukey algorithm [1] formula ${\displaystyle M\left(n\right)=\frac{N}{2}{log}_2\left(N\right)}$ for large values of N. The above results were programmed with MATLAB and are listed in Appendix 1

3. Complexity of Radix-4 DIT/DIF Algorithm

A MATLAB program has been written and is in Appendix 2, to compute the FFT of a set of discrete data: ${\displaystyle \left\{x_n\right\},n=1,2,3,\dots ,N}$ where $N=4^n$ and n is an integer.

The program uses a recursive call of a function that computes the FFT of the odd and even data points. The original program was modified to reduce the loop length by a factor of 4 by making use of the property ${\displaystyle e^{-j\frac{2\pi (i-1)}{N}}=-e^{-j\frac{2\pi \left(\frac{N}{4}i+1\right)}{N}}}$. Hence, the number of multiplications in the loop was halved [3].

In the special case when $N=4$, there are no multiplications to evaluate and only 8 additions. In the general case when $N=4^n$, let $M\left(n\right)$ represent the number of computed multiplications in the program.

The program has 4 function calls of the FFT of data which is half the original length and therefore contains $4\times M(n-1)$ multiplications. The for loop in the algorithm contains one multiplication which is repeated $N/4$ times (i.e., $4^{n-1}$ times) [5,9]. Hence, the number of multiplications in the program can be computed from the sequence:

$$M(n+1)=M\left(n\right)+4^n,n=1,2,3,4,\dots ,\mathrm{where}M\left(1\right)=0$$

(5)

The solution to this difference equation can be shown to be given by [10,11]:

$$M\left(n\right)={\displaystyle \frac{3N}{4}}[{log}_{4}N-1]$$

(6)

The number of additions in the program can be computed from the sequence:

$$A(n+1)=A\left(n\right)+4^n,n=1,2,3,4,\dots ,\mathrm{where}A\left(1\right)=0$$

(7)

The solution to this difference equation can be shown to be given by:

$$A\left(n\right)=2N{log}_{4}N=2N{\displaystyle \frac{{log}_{2}N}{{log}_{2}4}}=N{log}_{2}N$$

(8)

This sequence is tabulated for various values of n and N as shown in Table 2 below:

Compared to Radix-2, the Radix-4 algorithm, for large N, has reduced the number of multiplications by $62.5\%$, but the number of additions has remained unchanged [1,3].

4. Complexity of Bluestein

The Bluestein algorithm [2] is a digital filtering approach to determine the DFT of a set of data and is valid for any value of N [7,8].

Consider an input sequence $x\left(n\right)$, $0\le n\le N-1$, which is input to a digital filter with an impulse response

$h\left(n\right)=\left\{\begin{array}{ccc}0& \mathrm{if}& n<0\\ e^{-j\pi n^2/N}& \mathrm{if}& 0\le n\le 2N-1\\ 0& \mathrm{if}& n\ge 2N\end{array}\right.$

     The output of the filter in the range $N\le n\le 2N-1$ can be shown to be given by:

$$y\left(n\right)=\sum _{r=0}^{N-1}{x}_{r}h(n-r),N\le n\le 2N-1$$

(9)

By substituting $\mathrm{k}=\mathrm{n}-\mathrm{N}$ into (9), the output of the filter can be expressed in the form:

$$\begin{array}{c}\hfill y\left(k\right)=e^{j\pi \left(k^2/N\right)}{e}^{j\pi N}\sum _{r=0}^{N-1}\left\{x_r{e}^{j\pi \left(r^2/N\right)}\right\}{e}^{-j(2\pi /N)rk},0\le k\le N-1\end{array}$$

(10)

This equation implies that $\mathrm{y}\left(\mathrm{n}\right)$ is the weighted output of the N point DFT of the sequence $x_r{e}^{j\pi \left(r^2/N\right)}$. This result suggests that the DFT of a sequence can be obtained by first pre-multiplying the input sequence by the factor of ${\mathrm{e}}^{-\mathrm{j}\pi \left({\mathrm{n}}^2/N\right)}$ and inputting this sequence to a digital filter [10,11,12]. The output sequence can be obtained multiplying the filter output by a factor of $e^{-j\pi \left({(n-N)}^2/N\right)}{e}^{-j\pi N}$. This is called Bluestein’s algorithm [2]. It should be noted that the filter output can be determined from the filter input by using a convolution of the impulse filter response and the input data [13]. A MATLAB program has been written to determine the DFT of a random set of data which uses the built-in MATLAB convolution function [7]. However, the complexity of the built-in MATLAB convolution function is not known. Therefore, a MATLAB program was written to evaluate the convolution using the result that the convolution of a function can be determined from the inverse transform of the transform products of the impulse response and the input data [6].

The convolution calculation in our MATLAB program requires three calculations of the FFT transform of a signal, which is of length $2N$ [5,9]. However, we need not include the FFT transform of the impulse response in our complexity calculation, since this can be precalculated and is independent of the data [6].

As seen in the modified Bluestein MATLAB program, the number of multiplications can be summarized as follows:

Bluestein Multiplications: $\left\{\begin{array}{c}\mathrm{N}\mathrm{Scaling}\mathrm{Inputs}\\ 2\mathrm{FFTs}\mathrm{each}\mathrm{of}\mathrm{length}2\mathrm{N}\\ 2\mathrm{N}\mathrm{FFT}\mathrm{products}\\ \mathrm{N}\mathrm{Scaling}\mathrm{Outputs}\end{array}\right.$     

Since the FFT calculations in the Bluestein Algorithm used the Radix-4 algorithm [14,15], then the total number of multiplications computed, $M\left(N\right)$, is given by the formula:

$$\begin{array}{ccc}\hfill M\left(N\right)& =4N+2\times {\displaystyle \frac{3\left(2N\right)}{4}}\left[{log}_4\left(2N\right)-1\right]\hfill & \hfill =4N+{\displaystyle \frac{3N}{2}}\left[{log}_4\left(2N\right)-1\right]\end{array}$$

(11)

For large N, this formula tends to ${\displaystyle M\left(N\right)=\frac{3N}{2}{log}_4\left(N\right)}$ ${\displaystyle =\frac{3N}{4}{log}_2\left(N\right)}$. This is an increase of $50\%$ over the Radix-2 calculation [16].

In the case of the addition complexity of the Bluestein algorithm in our program, this count is equal to the count of 2 FFT’s of length $2N$, using the Radix-4 algorithm [5].

$$\mathrm{i}.\mathrm{e}.,A\left(n\right)=2\times \left(2N\right){log}_2\left(2N\right)=4N^2{log}_2\left(2N\right)$$

(12)

For large N, this value is four times that of the Radix-2 calculation.

It should be noted that the Bluestein algorithm is less efficient than the Radix-2 and Radix-4 algorithms [17], but has the advantage that the value of N does not need to be a power of 2 [1,9]. The above results were programmed with MATLAB and are listed in Appendix 2.

5. Theoretical and Mathematical Model of Bluestein Algorithm Using Three FFTs for Convolution

The Bluestein algorithm [2,3], also known as the chirp z-transform method, is a versatile approach for computing the Discrete Fourier Transform (DFT) of sequences where the length N is not restricted to powers of two [1]. This algorithm reformulates DFT as a convolution, enabling its computation via fast Fourier transforms (FFT), making it particularly efficient when N is not suitable for traditional radix-based FFT algorithms [6].

5.1. Formulation of the DFT

Let $x\left(n\right)$ be an input sequence of length N. The DFT of $x\left(n\right)$ is defined as:

$$X\left(k\right)=\sum _{n=0}^{N-1}x\left(n\right)e^{-j{\displaystyle \frac{2\pi }{N}}kn},0\le k\le N-1$$

(13)

5.2. Bluestein’s Reformulation

Bluestein’s [2] key insight was to rewrite the complex exponential term:

$$e^{-j2\pi kn/N}=e^{-j\pi kn/N}·e^{j\pi k^2/N}·e^{-j\pi n^2/N}$$

(14)

This allows us to express the DFT as:

$$X\left(k\right)=e^{j\pi k^2/N}\sum _{n=0}^{N-1}\left[x\left(n\right)e^{-j\pi n^2/N}\right]·\left[e^{-j\pi {(k-n)}^2/N}\right]$$

(15)

5.3. Convolution Formulation

The definition of the two sequences is:

$$a\left(n\right)=x\left(n\right)e^{-j\pi n^2/N},n=0,1,\dots ,N-1$$

(16)

$$b\left(n\right)=e^{-j\pi n^2/N},n=-(N-1),\dots 0,\dots N-1$$

(17)

By using the (16) and (17), the DFT can now be expressed as a convolution:

$$X\left(k\right)=e^{j\pi k^2/N}\sum _{n=0}^{N-1}a\left(n\right)b(k-n)$$

(18)

5.4. Implementation by Using 3 FFTs

To compute this convolution efficiently, we use the convolution theorem, which states that convolution in the time domain is equivalent to multiplication in the frequency domain [13,16]. The $a\left(n\right)$ and $b\left(n\right)$ can be extended with zeros to length $M\ge 2N-1$. Typically, M is chosen as the next power of 2 greater than or equal to $2N-1$ for efficient FFT computation. This leads to the following steps:

The FFTs of the padded sequences can be computed as:

$$A\left(k\right)=\mathrm{FFT}\left\{a\right(n\left)\right\}$$

(19)

$$B\left(k\right)=\mathrm{FFT}\left\{b\right(n\left)\right\}$$

(20)

Multiply the FFTs pointwise by using the (19) and (20):

$$C\left(k\right)=A\left(k\right)·B\left(k\right)$$

(21)

Then the (21) can compute the inverse FFT of the product:

$$c\left(n\right)=\mathrm{IFFT}\left\{C\right(k\left)\right\}$$

(22)

By extracting the first N points of $c\left(n\right)$ from the (22) and multiply by the chirp:

$$X\left(k\right)=e^{j\pi k^2/N}c\left(k\right),k=0,1,\dots ,N-1$$

(23)

5.5. Complexity Analysis

The complexity of this algorithm is dominated by the three FFT operations (two forward and one inverse) of length M [12]. Each FFT has a complexity of $O(MlogM)$. Since M is chosen to be the next power of 2 greater than $2N-1$, we have $M<4N$. Therefore, the overall complexity of the Bluestein algorithm is $O(NlogN)$ [8].

More precisely, if we use a radix-2 FFT algorithm, the number of complex multiplications is approximately:

$$3·{\displaystyle \frac{M}{2}}{log}_{2}M+2N$$

(24)

The first term accounts for the three FFTs, and the additional $2N$ multiplications come from the pre- and post-processing steps [13].

5.6. Advantages and Considerations

The Bluestein algorithm [2] allows computation of DFTs for any length N, including prime lengths. It maintains the $O\left(N{log}_{2}N\right)$ complexity of the FFT for all N [5]. The algorithm introduces some overhead due to zero-padding and additional multiplications, making it less efficient than specialized FFT algorithms for highly composite lengths. The flexibility in choosing M allows for potential optimizations in specific hardware implementations [14].

In conclusion, the Bluestein algorithm provides a versatile approach to computing DFTs of arbitrary lengths while maintaining the efficiency of the FFT [7]. Its reformulation of the DFT as a convolution, computed via three FFTs, makes it a valuable tool in signal processing applications where the input length is not a power of two or a highly composite number [15].

6. Discussion and Conclusion

The computational complexity analysis of the Radix-2, Radix-4, and Bluestein algorithms, as presented in this study, provides key insight into the trade-offs [7] inherent in selecting an appropriate algorithm for Discrete Fourier Transform (DFT) computations [14]. Each algorithm has its advantages and limitations, which depend on the specific application requirements and sequence length N [5].

The Radix-2 algorithm, fundamental in FFT computations, achieves an efficient computational complexity of $O\left(N{log}_{2}N\right)$ for sequence lengths that are powers of two. Its simplicity and iterative nature make it suitable for general-purpose applications where the input size N aligns with its constraints. However, when N deviates from a power of two, zero-padding becomes necessary, which can introduce computational overhead [16,18].

The Radix-4 algorithm extends the efficiency of Radix-2 by grouping data into quartets rather than pairs. This approach reduces the number of multiplicative operations required, achieving a complexity of ${\displaystyle \frac{3}{4}}N{log}_{4}N$, translating to significant computational savings for large N. However, the trade-off lies in the additional programming complexity and the restriction that N must be a power of four [17]. For applications where N meets this criterion, the Radix-4 algorithm provides a compelling choice with superior performance.

The Bluestein algorithm [2] emerges as a versatile solution for DFT computations with arbitrarily arbitrary N, including prime lengths or lengths that are not suitable for radix-based FFTs [17]. By reformulating the DFT as a convolution, Bluestein uses the FFT to compute this convolution efficiently [6]. However, the flexibility of the algorithm comes at the cost of increased computational complexity, particularly for large N. Using three FFTs of length $2N$, its complexity is approximately $9N{log}_2\left(2N\right)$, which is higher than the Radix-2 and Radix-4 algorithms for highly composite N. However, the algorithm’s ability to handle non-standard N without zero padding makes it indispensable in scenarios where input length restrictions cannot be relaxed [3].

A comparative analysis reveals that while the Radix-2 and Radix-4 algorithms excel in computational efficiency for their specific N constraints, they lack the flexibility of the Bluestein approach [17]. In contrast, Bluestein [2] provides a robust alternative for arbitrary N but at the expense of increased computational overhead and implementation complexity.

In practice, the choice of algorithm depends on the nature of the application. For general-purpose computations where N aligns with powers of two or four, Radix-2 and Radix-4 offer optimal performance [7]. For applications in digital signal processing and related fields, where N is dictated by the data acquisition process or hardware limitations, the flexibility of the Bluestein algorithm becomes essential despite its higher complexity.

Another consideration is the hardware implementation of these algorithms. Radix-based FFTs benefit from streamlined hardware designs because of their structured nature [17]. In contrast, reliance on the Bluestein algorithm on convolution computations introduces additional design complexity, which can impact real-time applications.

In conclusion, this study highlights the importance of understanding the computational and implementation trade-offs associated with each FFT algorithm. The Radix-2 and Radix-4 algorithms remain highly efficient for sequences with lengths that match their structural constraints, while the Bluestein algorithm provides unparalleled flexibility for arbitrary lengths [10]. Future work may explore hybrid approaches that combine the strengths of these algorithms or optimize hardware implementations for specific application domains, ensuring that the computational requirements of modern signal processing tasks are met effectively.

6.1. Future Work

Future research could focus on several promising directions to enhance the utility and efficiency of FFT algorithms. One potential avenue is the optimization of these algorithms for specific hardware architectures, such as Graphics Processing Units (GPUs) and Field-Programmable Gate Arrays (FPGAs). Using parallel processing capabilities and hardware-specific designs could significantly reduce computational overhead and improve execution speeds.

Another area of exploration is the development of hybrid FFT algorithms that combine the strengths of the Radix-2, Radix-4, and Bluestein approaches [17]. Such hybrids could adapt dynamically to input characteristics, providing a balance between computational efficiency and flexibility.

Furthermore, the application of machine learning techniques to predict the most efficient FFT method for given data lengths and constraints could prove valuable. Adaptive systems that select algorithms based on real-time input characteristics and available computational resources could revolutionize the usage of FFT in dynamic environments.

Finally, extending the application of these algorithms to emerging domains, such as quantum computing and big data analytics, presents exciting opportunities. Quantum FFTs could leverage quantum parallelism, while big data applications may benefit from distributed and cloud-based FFT implementations that handle massive datasets effectively.

These directions underscore the enduring relevance and potential of FFT research, ensuring its continued impact on signal processing and related fields.