Method and system for computer-aided processing of data samples using an N-point radix-p-Fast Fourier Transform (FFT).

Abstract

Method for computer-aided processing of data samples using an N-point radix p-fast Fourier transform (FFT) with a total of l transformation stages, wherein the output of each transformation stage i, with i = 1...,l, and an iteration range Ri, where Ri = Ri-1 / p with R0 = N, were computed in p-groups and Ri / Vsize vector iterations with Vsize as vector width, and wherein each register stores a vector of data values ​​computed by Vsize / DT, where DT is the data type, and wherein up to the penultimate transformation stage l-1 the transformation was performed in natural order, the transformation of the last transformation stage l comprising the steps: a) sequential loading (LRS) of the vector registers of the penultimate transformation stage i-1 by a vectorized operation and adhering to a pre-computed reorder index that differs from the natural order;b) Applying a radix-p butterfly (ABS) to the arranged vector registers; and c) Storing the output (STLS) of the radix-p butterfly, wherein the indices of the data values ​​of the last transformation stage l are specified in natural order, wherein at least two loaded vector registers are combined (CLVS) prior to applying the radix-p butterfly (ABS) to form wider vector registers for higher data throughput.

Description

[0001] The present invention relates to a method and a system for computer-aided processing of data samples using an N-point radix-p Fast Fourier Transform, FFT.

[0002] The Fast Fourier Transform (FFT) is a widely used algorithm for calculating the Discrete Fourier Transform (DFT) of a sequence or its inverse (IDFT). It is a cornerstone of digital signal processing and is widely used in various scientific, engineering, and technological fields. The FFT can be used in the following different areas: 1. Signal processing: The FFT is central to the processing of audio, video, and telecommunications signals. It enables efficient signal filtering, analysis, and transformation, enabling various applications such as noise reduction, equalization, and compression. 2. Medical imaging: In medical imaging procedures such as MRI and CT scans, FFT helps reconstruct and enhance images. It is used to analyze frequency components within the image and improve both clarity and detail. 3. Aerospace and radar systems: FFT plays a key role in synthetic aperture radar (SAR) technology and other radar systems. It aids in target acquisition, imaging, and tracking, enhancing the capabilities of modern defense and space research systems. 4. Weather forecasting: FFT provides important insights into weather patterns and helps in making accurate forecasts by analyzing meteorological data. 5. Financial analysis: In the financial sector, FFT is a tool for analyzing trends and patterns in various financial instruments, facilitating the development of trading algorithms and risk management strategies. 6. Music and audio processing: FFT enables the analysis of audio signals for music production, including pitch correction, equalization, and spectral analysis. It is an indispensable tool for both professional sound engineers and hobbyists. 7. Scientific research and computational chemistry: FFT helps in solving partial differential equations (PDEs) and analyzing complex data sets in research that includes areas such as quantum mechanics and molecular dynamics simulations. 8. Seismology: By analyzing seismic waves, FFT helps in the study of the Earth's interior and the prediction and analysis of earthquakes. 9. Network Analysis: FFT analyzes network traffic, enabling better management and security monitoring in modern communication systems. 10. Machine learning and data analysis: FFT is also used in machine learning for feature extraction and data preprocessing to enable efficient model training.

[0003] The widespread use of the FFT is due to its efficiency. Conventional DFT calculations require O(N 2) operations. In contrast, the FFT considerably reduces this complexity to O(N log(N)), making it a powerful tool for modern computations (see reference [P1]). Ultimately, it enables real-time processing and analysis in various applications. Through continuous research and development, its applicability is constantly expanding, making the FFT an indispensable element in the ever-expanding world of technology and science.

[0004] Numerous FFT algorithms (e.g., as described in references [P1, P2, P3, P4, P5, P6, P7]) can be found in digital signal processing. Due to its popularity, the Cooley-Tukey algorithm (see [P7]) is included in most FFT libraries optimized for specific hardware architectures. Despite the effectiveness of conventional approaches, certain drawbacks limit the full optimization potential of the FFT.

[0005] Conventional FFT implementations face challenges such as: - Requirement of a reordering step: Many typical FFT algorithms require a reordering step to rearrange the data into a specific sequence. This additional step increases the computational effort and complicates the design and execution of the algorithm. - Incomplete vectorization: Existing FFT methods often achieve only partial vectorization, especially in the reordering stages. Scalar operations frequently occur in these stages, which hinders full SIMD utilization and limits optimization options. - Implementation complexity: Reordering and partial vectorization lead to a more complex and error-prone implementation process. This complexity can impair efficiency, especially when adapting the algorithm to different hardware environments. - Optimization challenges: While various radix-based FFT implementations exist, optimal performance tuned for modern hardware has yet to be achieved. Conventional techniques often require complicated tuning and may not deliver the best possible performance. Radix algorithms typically implement the FFT by splitting the transformation into subproblems based on a specific basis, or radix.

[0006] The radix approach is a widely used solution, but it faces the problem of reordering or bit reversal, which can lead to additional computational overhead and complexity. This includes challenges such as: - In-place reordering and decoupling of computation: Radix algorithms often enable in-place computation by reordering within the existing data structure. This in-place technique can eliminate the need for subsequent reordering by using a bit-flip algorithm for the initial reordering of the elements. Furthermore, special radix algorithms can decouple computation from reordering, allowing computation to occur in the natural order and integrating reordering into the computation steps. - Specialized implementations and mixed-radix algorithms: Some specialized radix algorithms minimize or eliminate reordering by structuring the computation to conform to the natural order of the data. While this may increase the complexity of indexing or additional computations, it avoids the explicit reordering step. In mixed-radix algorithms, the problem of reordering is further mitigated by using different radixes, allowing for a flexible computation structure and reducing or eliminating the need for reordering. - Hardware capabilities and vectorization: The use of modern hardware with vectorized instructions is an efficient response to the challenge of reordering. Hardware-based shuffling can be more powerful than software-based solutions and corresponds to the innovative approach that fully vectorizes the FFT algorithm, eliminating the need for a separate reordering stage.

[0007] Radix-based FFTs, such as Radix-2, are based on a divide-and-conquer method that enables significant computational efficiency. This strategy requires a continuous subdivision of the problem, systematically processing pairs of data points and larger groups until the entire dataset is covered. In the initial stages of the Radix-2 FFT, the algorithm pairs neighboring points. In subsequent stages, the algorithm combines these pairs into larger groups, doubling in size with each step until the entire dataset has been processed. During this subdivision, the least significant bits in the indices are processed, resulting in the indices being arranged in a bit-reversed order (see [P1]). The bit-reversed ordering is not a mere artifact, but a structural consequence of the efficient computational approach.Although this order is mathematically clever, it does not match the original sequence, so a reordering step is required.

[0008] The bit-reversed order in radix-based FFTs arises directly from the divide-and-conquer technique that gives these algorithms their efficiency. The handling of this order varies between implementations and reflects a broader consideration of efficiency, complexity, and the specific requirements of the application.

[0009] The procedure described above is illustrated by Table 1, which shows the characteristics of a 4096-point radix-4 FFT. This example assumes systems offering vector support (SIMD) from 128 bits to 2048 bits, or possibly more. Table 1 shows that the first five stages (i = 1...5) are conveniently vectorizable because their iteration range exceeds the vector length for the float16, float32, or float64 data format. Table 1: 4096-point radix-4 FFT features Level i Iterationsb area Vector (register) Description 1 1024 128 - 2048 bits for F32 and F64 FFT transformation process 2 256 128 - 2048 bits for F32 and F64 3 64 128 - 2048 bits for F32 and F64 4 16 128 - 512 bits for F32128 - 1024 bits for F64 5 4 128 bits for F32128 - 256 bits for F64 6 1 - Rearrangement (final stage) - - Reorganization process

[0010] In the final stage, i = 6, vectorization is not possible. Therefore, scalar operations are required for further processing. Furthermore, the "reordering" process cannot be processed with vectors using SIMD operations.

[0011] Fig.Figure 1 illustrates the difficulties involved in vectorizing the final transformation stage, stage 6 in this example, and the subsequent reordering stage of radix-based FFTs. These difficulties arise from the inherent computational patterns of these stages and the properties of the SIMD architecture. In the final step of an FFT, operations must be performed on groups of data points that match the radix of the FFT (e.g., four adjacent elements in radix-4 FFTs).

[0012] In Fig. 1A the groups of data points are ν1Re,ν2Re,ν3Re,ν4Re and ν1Im,ν2Im,ν3Im,ν4Im , where the elements of each group in the left column Out2 / In3 are labeled Re[n], ..., Im[n], where n represents an index 0 ... 7 of the data values. As in Fig.As shown in Figure 1, the indices of the input In3, which correspond to the output Out3 of the previous stage, are in natural order (NO for short). However, SIMD operations are designed for maximum efficiency when dealing with contiguous data blocks. In Fig. In Figure 1A, the middle column labeled RV-B and IV-B illustrates that a uniform operation (a set of identical instructions) is not possible because the mathematical operation applied to the loaded elements of a group is not uniform for both the real-valued butterfly operation RV-B and the imaginary-valued butterfly operation IV-B. Furthermore, after applying the butterfly operations, the elements in the same groups or vectors are stored in bit-reversed order (BRO for short), as can be seen from the right column Out3.

[0013] Since these data groups (ie vector registers ν1Re,ν2Re,ν3Re,ν4Re ) are not normally arranged sequentially in memory, this misalignment poses a significant challenge to the simple application of SIMD. Furthermore, the reordering stage, which takes the data from output Out3 as input and reorders the data in bit-reversed order, is very complex. This subsequent stage involves a non-uniform operation, which contradicts the uniformity of SIMD operations, in which the same instruction is executed across multiple data points simultaneously. The reordering process requires knowledge of the position of each data item, which is different from the SIMD model, in which operations on multiple data points are performed in parallel.

[0014] The object of the present invention is to provide an improved method and system for computer-aided processing of data samples using an N-point radix-p Fast Fourier Transformation (FFT) that avoids the problems described above.

[0015] These objects are achieved by a method according to the features of method claim 1, a computer program product according to the features of claim 10 and a system according to the features of claim 11. Preferred embodiments are set out in the independent claims.

[0016] According to a first alternative of a first aspect of the invention, a method for computer-aided processing of data samples using an N-point radix-p Fast Fourier Transform (FFT) with a total of l transformation stages is proposed. The output of each transformation stage i, with i=1...,l, and an iteration range R i , where Ri = R i-1 / p with R0 = N, were divided into p-groups and R i / V size Vector iterations calculated, where V size is the vector width. Each register stores a vector of data values ​​defined by V size / DT, where DT is the data type. Up to the penultimate transformation level l-1, the transformation was performed in natural order. The transformation of the last transformation level l comprises the following steps: a) sequential loading of the vector registers of the penultimate transformation stage i-1 by a vectorized operation and observing a pre-calculated reordering index that differs from the natural order; b) applying a radix-p butterfly to the arranged vector registers; and c) Storing the output of the radix-p butterfly, where the indices of the data values ​​of the last transformation stage l are specified in natural order.

[0017] According to a second alternative of the first aspect of the invention, a method for computer-aided processing of data samples using an N-point radix-p Fast Fourier Transform (FFT) with a total of l transformation stages is proposed. The output of each transformation stage i, with i=1...,l, and an iteration range R i , where R i = R i-1 / p with R0 = N, were divided into p-groups and R i / V size Vector iterations calculated, where V size is the vector width. Each register stores a vector of data values ​​defined by V size / DT, where DT is the data type. Up to the penultimate transformation level l-1, the transformation was performed in natural order. The transformation of the last transformation level l comprises the following steps: a) Index-based loading and construction of pre-computed, i.e. required, vectors having a reordering index different from the natural order for the last transformation stage l; b Applying a radix-p butterfly (ABS) to the ordered vector registers; and c) Storing the output (STLS) of the radix-p butterfly, where the indices of the data values ​​of the last transformation stage l are given in natural order.

[0018] The proposed approach provides a transformative solution that seamlessly integrates reordering into the vectorization of the final step of the FFT. By eliminating the reordering step and fully vectorizing, this new FFT algorithm overcomes the limitations of existing methods and introduces a streamlined and efficient method that is adaptable to various applications.

[0019] In a further preferred embodiment of the first alternative, prior to applying the radix-p butterfly, at least two loaded vector registers are combined into wider vector registers to achieve higher data throughput. Before applying the radix-p butterfly, it should be understood that this step is performed between steps a) and b). In particular, a transposition may be performed to arrange the vectors according to the transformation of the last stage l, thereby providing a predetermined reordering index. In some embodiments, the transposition (PTS) may be a partial transposition. For example, a transposition or a partial transposition is performed on each of the first half of the combined vector registers and a second half of the combined vector registers to arrive at further combined vector registers with 2p data values ​​having a predetermined reordering index.

[0020] In a further preferred embodiment of the second alternative, reorder indices can be processed as input for the construction of index vectors in index-based load instructions.

[0021] According to a further embodiment, the radix-p based FFT is a radix-2 algorithm or a radix-4 algorithm, ie p=2 or p=4.

[0022] According to another embodiment, maintaining the pre-computed reordering index comprises reordering operations of the indices to achieve the natural order after the radix-p butterfly has been applied. In particular, the reordering operation comprises reordering in blocks of uniform arithmetic operations, allowing the reuse of vector registers for each iteration. Further reordering may comprise reordering of indices by using sequential or index-based load instructions to load p vector registers at predetermined indices. It may be advantageous to perform a p × p vectorized transpose to create input information for step b).

[0023] According to another embodiment, performing the partial transposition comprises transposing parts of the pre-arranged and combined vectors. For example, in the case of 256-bit vectors and the F32 data type, the partial transposition applies to the upper and lower 128 bits of the vector group.

[0024] According to a further embodiment, steps a) to c) are performed for each individual iteration.

[0025] According to a second aspect, a computer program product is proposed comprising instructions which, when a computer executes the program, cause the computer to carry out the steps of the method of one or more preferred embodiments.

[0026] According to a third aspect, a system for computer-aided processing of data samples using an N-point Radix-p Fast Fourier Transform, FFT, with a total of l transform stages is proposed, wherein the output of each transform stage i, with i=1...,l, and an iteration range R i , where R i = R i-1 / p with R0 = N, in p-groups and R i / V size Vector iterations were calculated, where V size is the vector width, and where each register stores a vector of data values ​​defined by V size / DT, where DT is the data type, and where up to the penultimate transformation stage l-1 the transformation was performed in natural order, the system comprising a processor configured to carry out one of the methods or one or more preferred embodiments thereof.

[0027] The approach described in this paper, which combines reordering with vectorization, provides a solution to this challenge and is consistent with ongoing efforts to optimize FFT calculations.

[0028] The invention is described in detail with reference to embodiments shown in the accompanying drawings. Fig. Figure 1A illustrates in a diagram the challenges of vectorizing the final stage and the reordering phase of an FFT implementation; Fig. Figure 1B shows the procedure for generating input vectors for the final transformation stage of the FFT, for both float32 and float64 data types; Fig. 2 a flowchart for the last stage of the radix-based FFTs by generating input vectors for the last transformation stage using the example of a radix-4 based 64-point FFT AVX256-bit Float 32; Fig. 3 is a table illustrating the last transformation stage and showing a bit-reversed output resulting from the penultimate transformation stage; Fig. 4 a table illustrates the final transformation stage using scalar optimization; Fig. 5 a table illustrates the final transformation stage using vectorization; Fig. 6 a table illustrating the final transformation stage, which involves vectorization and sorting; Fig. Figure 7 illustrates sequential SIMD loading to improve read performance; Fig. Figure 8 illustrates the 4x4 vectorized transposition to form the required vectors for the final transformation stage; Fig.Figure 9 shows a table illustrating the first iteration of the final transformation stage of a 64-point radix-4 FFT vectorized with 256 / 512 bits; Fig. Figure 10 is a table illustrating the second iteration of the final transformation stage of a 64-point radix-4 FFT vectorized with 256 / 512 bits; Fig. Figure 11 illustrates sequential SIMD loading to improve read performance for vectors 0 to 7; Fig. Figure 12 is a diagram illustrating the generation of larger vectors to improve transformation performance and throughput; Fig. Figure 13 is a table illustrating the final transformation stage of a 64-point radix-4 FFT vectorized to 512 bits for iteration 1; Fig. Figure 14 is a table illustrating the final transformation stage of a 64-point radix-4 FFT vectorized to 512 bits for iteration 2; Fig.15 a flowchart for the last stage of the radix-based FFTs by generating input vectors for the last transformation stage using the example of a radix-4 based 64-point FFT - ARM SVE 256-bit Float 32; Fig. 16 a flowchart for the last stage of the radix-based FFTs by generating input vectors for the last transformation stage using the example of a radix-4 based 64-point FFT AVX128-bit Float 32; Fig. 17 a flowchart for the last stage of the radix-based FFTs by generating input vectors for the last transformation stage using the example of a radix-4 based 64-point FFT AVX256-bit Float 64; Fig. 18 a flowchart for the final stage of the radix-based FFTs by generating input vectors for the final transformation stage using the example of a radix-4-based 64-point FFT AVX512-bit Float 64; and Fig.19 a flowchart for the final stage of the radix-based FFTs by generating input vectors for the final transformation stage using the example of a radix-4 based 64-point FFT - ARM NEON 128-bit Float32.

[0029] According to the present invention, the method is based on processing data samples using an N-point radix-p Fast Fourier Transform (FFT), which is well known to those skilled in the art. The N-point radix-p FFTs are based on SIMD operations to ensure efficient processing.

[0030] SIMD, an acronym for Single Instruction, Multiple Data, is a parallel processing architecture within a CPU that significantly increases computational efficiency by performing the same operation on multiple data points simultaneously. This contrasts with scalar operations, which process a single data point per operation. The effectiveness of SIMD is particularly evident in data-intensive tasks such as digital signal processing, multimedia applications, and scientific computing.

[0031] The key terms of SIMD are: - SIMD registers: The heart of SIMD capabilities are the registers, which range from 128 bits to 2048 bits in advanced architectures. The size of these registers determines the amount of data that can be processed in parallel. For example, a 128-bit SIMD register can process four 32-bit or two 64-bit operations simultaneously, while a 2048-bit register can handle sixteen 128-bit operations simultaneously. - Efficiency over scalar operations: SIMD's parallel processing capability contrasts with the scalar operations performed by traditional ALU (Arithmetic Logic Unit) components in CPUs. While ALUs process one data element per instruction, SIMD can process multiple data elements. This parallelism enables significant speed improvements for tasks with repetitive and parallelizable operations, such as matrix multiplications or vector operations, by reducing the number of instruction cycles. Because of these differences, SIMD and ALU units have different roles in a CPU. The ALU performs basic arithmetic and logical operations, usually with scalar values. In contrast, SIMD units are designed for "wide" operations and utilize data-level parallelism.This distinction makes SIMD units particularly effective for tasks with high data parallelism, while ALUs are more suitable for general calculations where such parallelism is lacking. - Application-specific optimization: The effectiveness of SIMD depends heavily on the nature of the task. Applications with repetitive and regular calculations can achieve significant performance gains with SIMD. However, these benefits are not necessarily present for tasks with irregular data patterns or those that require extensive branching.

[0032] Previously, the final stage of an N-point radix-p FFT could not be processed using SIMD operations because sequential loading of data elements and unified operations were not possible. The proposed method overcomes these disadvantages.

[0033] Vectorizable loops are programming loops that can be optimized through vectorization. Operations within the loop are performed simultaneously on multiple data points using SIMD instructions. This technique improves computational efficiency, especially for data-intensive tasks.

[0034] The most important concepts for vectorizable loops are: - Unified operations: The loop performs the same operation on several independent data elements. - Data independence: No iteration of the loop depends on the results of another, so parallel execution is possible. - Regular data structures: Efficient vectorization often requires storing data in regular, contiguous memory blocks.

[0035] Vector-vector multiplication is a classic example of a vectorizable loop. Consider two equal-sized arrays (vectors) A and B, where each element of the resulting array C is the product of the corresponding elements in A and B. In this case, each multiplication operation is independent of the others. The loop iterates over arrays A and B, multiplying each pair of elements. This independence makes it an ideal candidate for vectorization, where a SIMD processor can calculate multiple products simultaneously, significantly speeding up the operation.

[0036] The Fast Fourier Transform (FFT) is a widely used algorithm that efficiently computes the Discrete Fourier Transform (DFT) of a sequence x(n). The following formulas are used to directly compute the DFT: XF(k)=∑0N−1x(n)WNkn n=0,1,2,…,N−1 WNkn=exp(−j2πknN)=cos(2πknN)+j sin(2πknN)

[0037] The direct calculation of an N-point DFT requires almost O(N 2) complex arithmetic operations. An arithmetic operation involves multiplication and addition. However, by developing efficient algorithms, such as radix-2 or radix-4, the complexity of FFT processing can be significantly reduced. In general, fast algorithms reduce the computational complexity of an N-point DFT to about N log2 (N) complex arithmetic operations. Other advantages include smaller memory requirements and lower computational errors due to finite-bit arithmetic (multiplication / division and addition / subtraction are implemented with finite word lengths for practicality). Fast algorithms have contributed to DFT implementation by DSP chips. FFT algorithms can be implemented in several ways, with the decimation-in-time (DIT) 1 and decimation-in-frequency (DIF) 2 methods being the most commonly used for the radix-2 and radix-4 algorithms.Both methods lead to the same result, but differ in their approach to decomposing and processing the sequence.

[0038] The Radix-2 Decimation-In-Frequency (DIF) Fast Fourier Transform (FFT) is an efficient algorithm for calculating the Discrete Fourier Transform (DFT) of a sequence. It is particularly suitable for sequences whose length is a power of two. The key features of the Radix-2 DIF FFT can be outlined as follows: - Frequency decimation: The algorithm decomposes the original DFT problem (size) by systematically decomposing the frequency spectrum. First, the entire spectrum is divided into two halves, representing the lower and upper halves of the frequency range. This process is repeated recursively, dividing the spectrum into smaller frequency bins. - Butterfly operations: At each stage of the decomposition, the Radix-2 DIF FFT algorithm employs a computational element called a "butterfly operation." These operations include basic arithmetic calculations such as addition, subtraction, and multiplication with complex twiddle factors. WNkn. The twiddle factors are complex exponential functions that are fundamental to the Fourier transform. - In-place computation: The algorithm uses the same memory space for input and output data, significantly reducing memory requirements. It uses small local computation buffers for efficiency and caching.

[0039] The number of transformation stages of a radix-2 FFT can be given by 2 l = N → l levels. The complex-valued radix-2 butterfly is given in equations (3) and (4). xSi(n)=xSi−1(n)+xSi−1(n+N2) xSi(n+N2)=(xSi−1(n)+xSi−1(n+N2))WNkn for n = 0, 1, 2, ..., N i - 1.

[0040] x Si (n) represents the transformation of level i for i = 1, 2, ..., l. x So (n) are the complex input samples, given by: n(n)=re(n)+j im(n)

[0041] N i is the iteration domain of stage i. Equation (5) represents the complex natural format, in which the real and imaginary components are stored sequentially as pairs in memory. The invention applies analogously to the split-complex format, in which the real and imaginary parts are stored sequentially in separate memory areas.

[0042] As already explained, radix FFT algorithms are based on a divide-and-conquer approach, which allows for significant computational efficiency. This strategy requires a continuous subdivision of the problem, systematically processing pairs of data points and larger groups until the entire dataset is covered.

[0043] In the first stages of the radix-2 FFT, the algorithm pairs neighboring points. In subsequent stages, the algorithm combines these pairs into larger groups, doubling in size with each step until it has processed the entire data set. This division involves processing the least significant bits in the indices, resulting in a result in which the indices are arranged in a bit-reversed order (see Fig.1). The bit-reversed order is not a mere artifact, but a structural consequence of the efficient computational approach. While mathematically ingenious, it does not match the original sequence, so a reordering step is required.

[0044] There are various strategies to address the challenge of bit reordering in radix-based FFTs. Some algorithms incorporate this reordering into the computation and process it as part of the normal algorithm flow. Others perform it as a discrete step at the end of the transformation process.

[0045] The solution described below vectorizes the final transformation stage of the FFT algorithm and simultaneously integrates the reordering stage, eliminating the need for a separate reordering process. This approach is explained in detail in Table 2, where NO stands for natural ordering and BRO for bit-reversed ordering. X[s] with s = 1...7 represents data values, where s is the index of the input data values ​​(left column). Table 2: Last transformation stage - bit-inverted output NO Butterfly operations BRO 0 x[0] + x[1] 0 1 x[0] - x[1] 4 2 x[2] + x[3] 2 3 x[2] - x[3] 6 4 x[4] + x[5] 1 5 x[4] - x[5] 5 6 x[6] + x[7] 3 7 x[6] - x[7] 7

[0046] For an 8-point FFT, the input is taken to BRO, resulting in an output in the natural order NO. Furthermore, this method orders arithmetic operations in uniform groups of four, as shown in Table 3. This arrangement facilitates the implementation of vectored additions and subtractions. A notable improvement in this method is that only two vectors are needed to perform the final transformation and reordering steps in a single step. This optimization significantly increases data throughput and FLOPS. However, one challenge remains: the inability to use efficient sequential vectored load instructions to construct the input vectors, since index-based load instructions are comparatively slow and therefore undesirable on x86 and ARM Neon systems.On SVE, index-based load instructions are efficient and can replace transposition functions. A solution is proposed that reorders the data accordingly after two sequential loads. This method is described in . Fig. 1B is shown as an example and relies only on fast unpacking operations. Table 3: Last transformation stage - natural order BRO Butterfly operations NO 0 x[0] + x[1] 0 4 x[4] + x[5] 1 2 x[2] + x[3] 2 6 x[6] + x[7] 3 1 x[0] - x[1] 4 3 x[2] - x[3] 6 5 x[4] - x[5] 5 7 x[6] - x[7] 7

[0047] Table 3 shows that this strategic reordering of the data allows to continue to use efficient loading methods while achieving the desired vectorization and reordering in the final FFT transformation stage.

[0048] Like the radix-2 approach, the radix-4 DIF-FFT algorithm provides an effective method for computing the DFT. This algorithm is particularly well-suited for processing sequences whose length is a power of four. A key distinguishing feature of the radix-4 algorithm lies in its butterfly operations: At each stage of the decomposition process, the butterfly operations involve basic arithmetic calculations, as shown in equation (6): [xSi(n)xSi(n+N4)xSi(n+N2)xSi(n+3N4)]=[1111WNkn(1−j−1j)WN2n(1− 11−1)WN3n(1j−1−j)][xSi−1(n)xSi−1(n+N4)xSi−1(n+N2)xSi−1(n+3N4)]

[0049] These butterfly operations in radix-4 are more complex than those in radix-2 because the fewer stages required for the decomposition involve more data, as radix-4 processes four data points simultaneously, compared to radix-2's two. This property makes radix-4 inherently more efficient for sequences of appropriate length (powers of four), as the same transformation can be achieved with fewer operations. The efficiency of the radix-4 FFT thus arises from this strategic combination of advanced butterfly operations and a reduced number of decomposition stages. The number of transformation stages of a radix-4 FFT can be reduced by 4 l = N → l stages can be determined.

[0050] The described method refers to the radix-4 FFT algorithm. In a 16-point radix-4 FFT, the output is initially presented in bit-reversed order (BRO), which mirrors the situation in radix-2 FFT, where a reordering step is required to achieve the natural order of the output. This similarity in the output format of radix-2 and radix-4 FFTs requires additional focus on the reordering process to ensure efficient data processing and output representation.

[0051] To further shed light on the mechanics of this new reordering process, it is instructive to examine the final transformation stage of a 64-point FFT. The 64-point FFT is only an illustrative example that presents a more complex scenario and highlights the challenges and solutions associated with managing larger data sets. By examining this stage, one can better understand how the reordering process is set up and optimized for a larger FFT, such as the 64-point FFT, compared to the simpler case of a 16-point FFT.

[0052] The reordering process in the final stage of the 64-point FFT is crucial for aligning the output in a naturally ordered sequence. This step is important to ensure that the FFT output is presented in a format that meets the requirements of conventional data processing and analysis. By focusing on this stage, we gain insights into the practical application of the reordering process and its impact on the overall efficiency and effectiveness of the FFT algorithm, especially in more complex and data-intensive scenarios.

[0053] Fig. Figure 2 illustrates the basic steps of the proposed method in the context of a 64-point FFT for processing the final transformation stage of complex-valued data based on Float32. The proposed approach in Fig.2 uses eight 128-bit registers R0,..., R7. The BRO look-up table (LUT) is used to load the vectors required for the transformation (represented by LRS). The highlighted cell of each register (R0, ..., R7) represents the base index for the load instructions. This loading scheme deviates from the natural order. The next step, CLVS, is a pairwise combination of the vectors R0 with R4 to form V0, R1 with R5 to form V1, R2 with R6 to form V2, and R3 with R7 to form V3 to create a larger vector for higher data throughput and more FLOPS in the butterfly kernel (denoted by INLS). Subsequently, the partial transpose PTS is performed on the halves of the generated structure, which is then shifted to the radix butterfly (ABS). Finally, the data is sequentially stored back into memory (see OULTS) in NO using the designated memory indices (highlighted cells in the first row).This efficient handling of output ensures that the entire process is streamlined and unnecessary memory operations are avoided.

[0054] In the described case (Float32 data and 256-bit vectors) of 64-point FFTs, a total of two iterations are required to complete the operation, where Fig. 2 shows one of these iterations.

[0055] The present solution eliminates a reordering stage and simultaneously vectorizes the final step. This approach results in a fully vectorized FFT algorithm that discards the scalar overhead and achieves optimal SIMD utilization rates.

[0056] In detail, the process for providing the pre-calculated indices used in the LRS step of Fig. 2 are loaded, carried out by a series of specific operations which are carried out with reference to the Fig.3 to 14. The starting point for this procedure is the penultimate transformation stage .i = l - 1

[0057] The output of the penultimate transformation stage l - 1 is in bit-reversed order BRO, which hinders the vectorization of the last transformation stage. This limitation arises because butterfly operations, which are essential for this stage, require processing adjacent samples, which deviates from conventional SIMD capabilities and leads to inefficiency. Furthermore, the non-uniformity of arithmetic operations, e.g., + - - +, complicates the implementation of SIMD instructions such as vectorized addition or subtraction, making them infeasible or requiring significant overhead. The first optimization consists in reordering the operations to achieve a natural order. The result is shown in the table in Fig. 3 shown

[0058] The Fig. The modified method shown in Figure 4 leads to natural order results. However, it introduces inefficiencies due to redundant load instructions (as indicated by the lines at BRO = 0, 2, 1, 3) and the exclusive use of scalar operations for butterfly computations. This approach requires 64 iterations to solve the problem (see the last column). The following optimization strategy consists of reordering and consolidating similar iterations to improve efficiency. This consolidation aims to reduce the total number of iterations and minimize redundancy, thereby improving computational efficiency.

[0059] This approach significantly improves the conversion and reordering processes by using 128-bit vectors, which is achieved by reordering into blocks of uniform arithmetic operations (in Fig.5 by the block UAO). Furthermore, vectors are efficiently reused multiple times in each iteration, optimizing data processing (highlighted by RV0, RV1, RV2, RV3). This leads to a significant reduction in the total number of iterations to only four. In the following phase, the focus is on further improving write performance. This is achieved by reordering the structure again, a strategic move aimed at streamlining the data output process and maximizing overall computational efficiency. The result is shown in Fig. 6 (marked by the block RIO).

[0060] To generate input vectors, reordering indices from the look-up table are used, with the focus on optimal efficiency ( Fig.6). This is achieved by using sequential load instructions, while index-based load instructions are deliberately avoided due to their relative inefficiency.

[0061] However, some systems can efficiently process index-based load instructions. In such cases, the index vectors used to load the vectors are formed by mapping BRO indexes. This creates the required data structure for the final transformation stage. Fig.Figure 15 shows an example for ARM SVE systems. The first step is to map the BRO indices from the look-up table (BMI). A predictor pg sets the register or vector width. In the given example, the vector length is set to 256 bits. Then, the vectors v0, v1, v2, and v3 are generated using the predictors pg, the input data, and the pre-built index vectors from the BRO look-up table (LUT). The four vectors are shifted into the radix kernel (ABS) and then stored in memory (OUTLS) in the same form as previously described.

[0062] As in Fig. 7, we load four vectors with indices 0, 32, 16 and 48. Then a 4x4 vectorized transposition is performed ( Fig.8). This step is crucial because it constructs the vectors for the final transformation stage, ensures streamlined processing, and increases the overall efficiency of the computation.

[0063] Starting with the float32 data type, performance can be further increased by merging 128-bit vectors into larger 256-bit vectors. This approach leverages the ability to efficiently process more data in a single operation. As highlighted in boxes CRV0, CRV1, CRV2, and CRV3, the vectors x[0]x[0] and x[8]x[8] can be combined into a single 256-bit vector (the combined vector from Fig. 2). The efficiency gains of this method, which include a lower number of iterations and higher data throughput, are shown in the results in Fig.9 (for iteration 1) and 10 (for iteration 2). These tables demonstrate the improved data processing efficiency achieved by this vector combination strategy.

[0064] The one in the Fig. The merge or combination process described in Figures 8 to 10 involves the sequential loading of eight 128-bit vectors for Float32 (f32) data types or eight 256-bit vectors for Float64 (f64) data types using the specified indices ( Fig. 11). These vectors are then combined into larger vectors, resulting in 256-bit vectors for f32 or 512-bit vectors for f64.

[0065] A partial transposition generates the necessary vectors, such as the boxes labeled CRV0, CRV1, CRV2, CRV3 in Fig.9. Combining and transposing vectors significantly improves performance, resulting in a remarkable reduction in the required iterations to only two. This streamlined approach optimizes processing efficiency through effective management of data types and vector sizes.

[0066] Similarly, we can optimize performance with even larger registers, such as the Fig. 13 and Fig. 14. This approach requires only a single iteration and provides the best performance. The four vectors can be generated by four partial transpositions.

[0067] Fig. Figures 16 to 18 show further examples for AVX128-bit Float 32, AVX256-bit Float 64 and AVX512-bit Float 64. Fig. Figure 19 shows another example of a 64-point FFT - ARM SVE 256-bit, illustrating the techniques described above.

[0068] As described, the inventive method vectorizes the final FFT transformation stage and eliminates the reordering stage by integrating it directly into it. This approach simplifies the overall process and makes it more efficient. By addressing the challenges typically associated with the final stage and the reordering stage of the FFT, an algorithm can be implemented that is fully vectorized and optimized for modern hardware. Such advances can have a significant impact on various real-time applications and open up new avenues for research and development in this field.

[0069] The inventors used various hardware configurations, including ground-based servers, to evaluate the implementation. Single- and double-precision one-dimensional complex transforms were used in the evaluation. Double-precision FFTs enable more accurate calculations, albeit at the expense of computational speed and memory.

[0070] These evaluations demonstrate that the proposed FFT implementation represents an effective and efficient alternative to existing FFT libraries. The generic code design enabled efficient testing on x86 and ARM platforms, demonstrating its flexibility and potential for broader applications. References [P1] KR Rao, DN Kim, and J.-J. Hwang. 2011. Fast Fourier Transform - Algorithms and Applications. Springer Science+Business Media BV, Dordrecht. https: / / doi.org / 10.1007 / 978-1-4020-6629-0 [P2] L. Bluestein. 1970. A linear filtering approach to the computation of discrete Fourier transform. IEEE Transactions on Audio and Electroacoustics 18, 4 (1970), 451-455. https: / / doi.org / 10.1109 / TAU. 1970.1162132 [P3] P. Duhamel and Henk Hollmann. 1984. ‚Split radix‘ FFT algorithm. Electronics Letters 20 (02 1984), 14 - 16. https: / / doi.org / 10.1049 / el: 19840012 [P4] D. Kolba and T. Parks. 1977. A prime factor FFT algorithm using high-speed convolution. IEEE Transactions on Acoustics, Speech, and Signal Processing 25, 4 (1977), 281-294. https: / / doi.org / 10.1109 / TASSP.1977.1162973 [P5] Zhihao Li, Haipeng Jia, Yunquan Zhang, Tun Chen, Liang Yuan, Luning Cao, and XiaoWang. 2019. AutoFFT: A Template-Based FFT Codes Auto-Generation Framework for ARM and X86 CPUs. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (Denver, Colorado) (SC '19). Association for Computing Machinery, New York, NY, USA, Article 25,15 pages. https: / / doi.org / 10.1145 / 3295500.3356138 [P6] Zhuo Qian and Martin Margala. 2016. Low-Power Split-Radix FFT Processors Using Radix-2 Butterfly Units. IEEE Transactions on Very Large Scale Integration (VLSI) Systems 24, 9 (2016), 3008-3012. https: / / doi.org / 10.1109 / TVLSI.2016.2544838 [P7] Sergio Vázquez, Margarita Amor, and Basilio B. Fraguela. 2019. Portable and efficient FFT and DCT algorithms with the Heterogeneous Butterfly Processing Library. J. Parallel and Distrib. Comput. 125 (2019), 135-146. https: / / doi.org / 10.1016 / j.jpdc.2018.11.011 QUOTES CONTAINED IN THE DESCRIPTION

[0000] This list of documents submitted by the applicant was generated automatically and is included solely for the convenience of the reader. This list is not part of the German patent or utility model application. The DPMA assumes no liability for any errors or omissions. Cited non-patent literature

[0000] KR Rao, DN Kim, and J.-J. Hwang. 2011. Fast Fourier Transform - Algorithms and Applications. Springer Science+Business Media BV, Dordrecht. https: / / doi.org / 10.1007 / 978-1-4020-6629-0

[0070] L. Bluestein. 1970. A linear filtering approach to the computation of discrete Fourier transform. IEEE Transactions on Audio and Electroacoustics 18, 4 (1970), 451-455. https: / / doi.org / 10.1109 / TAU. 1970.1162132

[0070] P. Duhamel and Henk Hollmann. 1984. ‚Split radix‘ FFT algorithm. Electronics Letters 20 (02 1984), 14 - 16. https: / / doi.org / 10.1049 / el: 19840012

[0070] D. Kolba and T. Parks. 1977. A prime factor FFT algorithm using high-speed convolution. IEEE Transactions on Acoustics, Speech, and Signal Processing 25, 4 (1977), 281-294. https: / / doi.org / 10.1109 / TASSP.1977.1162973

[0070] Zhihao Li, Haipeng Jia, Yunquan Zhang, Tun Chen, Liang Yuan, Luning Cao, and XiaoWang. 2019. AutoFFT: A Template-Based FFT Codes Auto-Generation Framework for ARM and X86 CPUs. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (Denver, Colorado) (SC '19). Association for Computing Machinery, New York, NY, USA, Article 25,15 pages. https: / / doi.org / 10.1145 / 3295500.3356138

[0070] Zhuo Qian and Martin Margala. 2016. Low-Power Split-Radix FFT Processors Using Radix-2 Butterfly Units. IEEE Transactions on Very Large Scale Integration (VLSI) Systems 24, 9 (2016), 3008-3012. https: / / doi.org / 10.1109 / TVLSI.2016.2544838

[0070] Sergio Vázquez, Margarita Amor, and Basilio B. Fraguela. 2019. Portable and efficient FFT and DCT algorithms with the Heterogeneous Butterfly Processing Library. J. Parallel and Distrib. Comput. 125 (2019), 135-146. https: / / doi.org / 10.1016 / j.jpdc.2018.11.011

[0070]

Claims

[1] Method for the computer-aided processing of data samples using an N-point Radix-p-Fast Fourier Transform, FFT, with a total of l transformation stages, where the output of each transformation stage i, with i=1...,l, and an iteration range R i , where R i = R i-1 / p with R0 = N, in p-groups and R i / V size Vector iterations with V size as vector width, and where each register stores a vector of data values ​​defined by V size / DT is calculated, where DT is the data type, and where up to the penultimate transformation level l-1 the transformation was performed in natural order, where the transformation of the last transformation level l comprises the steps: a) sequential loading (LRS) of the vector registers of the penultimate transformation stage i-1 by a vectorized operation and respecting a pre-calculated reordering index that differs from the natural order; b) applying a radix-p butterfly (ABS) to the arranged vector registers; and c) Storing the output (STLS) of the radix-p butterfly, where the indices of the data values ​​of the last transformation stage l are given in natural order. [2] Method for the computer-aided processing of data samples using an N-point Radix-p-Fast Fourier Transform, FFT, with a total of l transformation stages, where the output of each transformation stage i, with i=1...,l, and an iteration range R i , where R i = R i-1 / p with R0 = N, in p-groups and R i / V size Vector iterations with V sizeas vector width, and where each register stores a vector of data values ​​defined by V size / DT is calculated, where DT is the data type, and where up to the penultimate transformation level l-1 the transformation was performed in natural order, where the transformation of the last transformation level l comprises the steps: a) Index-based loading and construction of precomputed vectors that have a reordering index different from the natural order for the last transformation stage l; b Applying a radix-p butterfly (ABS) to the ordered vector registers; and c) Storing the output (STLS) of the radix-p butterfly, where the indices of the data values ​​of the last transformation stage l are given in natural order. [3] The method of claim 1, wherein prior to applying the radix-p butterfly (ABS), at least two loaded vector registers are combined (CLVS) to form wider vector registers for higher data throughputs; [4] A method according to claim 3, wherein a transposition (PTS) is performed to arrange the vectors according to the transformation of the last stage l, thereby obtaining a predetermined reordering index. [5] The method of claim 4, wherein the transposition (PTS) is a partial transposition. [6] The method of claim 2, wherein reorder indices are processed as input for the construction of index vectors in index-based load instructions. [7] A method according to any one of the preceding claims, wherein the radix-p based FFT is a radix-2 algorithm or a radix-4 algorithm. [8] Method according to one of the preceding claims, wherein steps a) to c) are carried out for each individual iteration. [9] A method according to any one of the preceding claims, wherein maintaining the pre-computed reordering index comprises reordering operations of the indices to achieve the natural ordering after the radix-p butterfly has been applied. [10] The method of claim 9, wherein the reordering operation comprises reordering into blocks of uniform arithmetic operations that allows reuse of vector registers for each iteration. [11] The method of claim 9 or 10, wherein further reordering comprises reordering indices by using sequential load instructions or index-based load instructions to load p vector registers at given indices. [12] A method according to any one of claims 9 to 11 in conjunction with claim 3, wherein a p × p vectorized transposition is performed to create input information for step b). [13] A method according to any one of claims 9 to 12 in conjunction with claim 3, wherein performing the partial transposition comprises transposing parts of the pre-arranged and combined vectors. [14] The method of claim 13, wherein the interleaving of data values ​​is repeated, the number of iterations depending on the number N of N-points and the radix p. [15] A computer program product comprising instructions which, when executed by a computer, cause the computer to carry out the steps of the method according to any one of claims 1 to 14. [16] System for computer-aided processing of data samples using an N-point Radix-p Fast Fourier Transform, FFT, with a total number of l transform stages, where the output of each transform stage i, with i=1...,l, and an iteration range R i , where R i = R i-1 / p with R0 = N, in p-groups and R i / V size Vector iterations were calculated, where V size is the vector width, and where each register stores a vector of data values ​​defined by V size / DT is calculated, where DT is the data type, and where up to the penultimate transformation stage l-1 the transformation has been performed in natural order, the system comprising a processor configured to perform the following steps to perform the transformation of the last transformation stage l: a) sequential loading (LRS) of the vector registers of the penultimate transformation stage i-1 by a vectorized operation and respecting a pre-calculated reordering index that differs from the natural order; b) applying a radix-p butterfly (ABS) to the arranged vector registers; and c) Storing the output (STLS) of the radix-p butterfly, where the indices of the data values ​​of the last transformation stage l are given in natural order. [17] System for computer-aided processing of data samples using an N-point Radix-p-Fast Fourier Transform, FFT, with a total number of l transform stages, where the output of each transform stage i, with i=1...,l, and an iteration range R i , where R i = R i-1 / p with R0 = N, in p-groups and R i / V size Vector iterations were calculated, where V sizeis the vector width, and where each register stores a vector of data values ​​defined by V size / DT is calculated, where DT is the data type, and where up to the penultimate transformation level l-1 the transformation has been performed in natural order, the system comprising a processor configured to perform the following steps to perform the transformation of the last transformation level l: a) Index-based loading and construction of precomputed vectors that have a reordering index different from the natural order for the last transformation stage l; b Applying a radix-p butterfly (ABS) to the ordered vector registers; and c) Storing the output (STLS) of the radix-p butterfly, where the indices of the data values ​​of the last transformation stage l are given in natural order. [18] The system of claim 16 or 17, wherein the processor is configured to perform a method according to any one of claims 2 to 14.