Vector-radix FFT algorithm explained

The vector-radix FFT algorithm, is a multidimensional fast Fourier transform (FFT) algorithm, which is a generalization of the ordinary Cooley–Tukey FFT algorithm that divides the transform dimensions by arbitrary radices. It breaks a multidimensional (MD) discrete Fourier transform (DFT) down into successively smaller MD DFTs until, ultimately, only trivial MD DFTs need to be evaluated.[1]

The most common multidimensional FFT algorithm is the row-column algorithm, which means transforming the array first in one index and then in the other, see more in FFT. Then a radix-2 direct 2-D FFT has been developed,[2] and it can eliminate 25% of the multiplies as compared to the conventional row-column approach. And this algorithm has been extended to rectangular arrays and arbitrary radices,[3] which is the general vector-radix algorithm.

Vector-radix FFT algorithm can reduce the number of complex multiplications significantly, compared to row-vector algorithm. For example, for a

NM

element matrix (M dimensions, and size N on each dimension), the number of complex multiples of vector-radix FFT algorithm for radix-2 is
2M-1
2M

NMlog2N

, meanwhile, for row-column algorithm, it is
MNM
2

log2N

. And generally, even larger savings in multiplies are obtained when this algorithm is operated on larger radices and on higher dimensional arrays.[3]

Overall, the vector-radix algorithm significantly reduces the structural complexity of the traditional DFT having a better indexing scheme, at the expense of a slight increase in arithmetic operations. So this algorithm is widely used for many applications in engineering, science, and mathematics, for example, implementations in image processing,[4] and high speed FFT processor designing.[5]

2-D DIT case

As with the Cooley–Tukey FFT algorithm, the two dimensional vector-radix FFT is derived by decomposing the regular 2-D DFT into sums of smaller DFT's multiplied by "twiddle" factors.

A decimation-in-time (DIT) algorithm means the decomposition is based on time domain

x

, see more in Cooley–Tukey FFT algorithm.

We suppose the 2-D DFT is defined

X(k1,k2)=

N1-1
\sum
n1=0
N2-1
\sum
n2=0

x[n1,n2]

k1n1
W
N1
k2n2
W
N2

,

where

k1=0,...,N1-1

,and

k2=0,...,N2-1

, and

x[n1,n2]

is an

N1 x N2

matrix, and

WN=\exp(-j2\pi/N)

.

For simplicity, let us assume that

N1=N2=N

, and the radix-

(r x r)

is such that

N/r

is an integer.

Using the change of variables:

ni=rpi+qi

, where

pi=0,\ldots,(N/r)-1;qi=0,\ldots,r-1;

ki=ui+viN/r

, where

ui=0,\ldots,(N/r)-1;vi=0,\ldots,r-1;

where

i=1

or

2

, then the two dimensional DFT can be written as:[6]

X(u1+v1N/r,u2+v2

r-1
N/r)=\sum
q1=0
r-1
\sum
q2=0

\left[

N/r-1
\sum
p1=0
N/r-1
\sum
p2=0

x[rp1+q1,rp2+q2]

p1u1
W
N/r
p2u2
W
N/r

\right]

q1u1+q2u2
W
N
q1v1
W
r
q2v2
W
r

,

The equation above defines the basic structure of the 2-D DIT radix-

(r x r)

"butterfly". (See 1-D "butterfly" in Cooley–Tukey FFT algorithm)

When

r=2

, the equation can be broken into four summations, and this leads to:[1]

X(k1,k2)=S00(k1,k2)+S01(k1,k2)

k2
W
N

+S10(k1,k2)

k1
W
N

+S11(k1,k2)

k1+k2
W
N
for

0\leqk1,k2<

N
2
,

where

Sij(k1,k2)=\sum

N/2-1
n1=0
N/2-1
\sum
n2=0

x[2n1+i,2n2+j]

n1k1
W
N/2
n2k2
W
N/2
.

The

Sij

can be viewed as the

N/2

-dimensional DFT, each over a subset of the original sample:

S00

is the DFT over those samples of

x

for which both

n1

and

n2

are even;

S01

is the DFT over the samples for which

n1

is even and

n2

is odd;

S10

is the DFT over the samples for which

n1

is odd and

n2

is even;

S11

is the DFT over the samples for which both

n1

and

n2

are odd.

Thanks to the periodicity of the complex exponential, we can obtain the following additional identities, valid for

0\leqk1,k2<

N
2
:
Xl(k
1+N
2

,k2r)=S00(k1,k2)+S01(k1,k2)

k2
W
N

-S10(k1,k2)

k1
W
N

-S11(k1,k2)

k1+k2
W
N
;

Xl(k1,k

2+N
2

r)=S00(k1,k2)-S01(k1,k2)

k2
W
N

+S10(k1,k2)

k1
W
N

-S11(k1,k2)

k1+k2
W
N
;
Xl(k
1+N
2
,k
2+N
2

r)=S00(k1,k2)-S01(k1,k2)

k2
W
N

-S10(k1,k2)

k1
W
N

+S11(k1,k2)

k1+k2
W
N
.

2-D DIF case

Similarly, a decimation-in-frequency (DIF, also called the Sande–Tukey algorithm) algorithm means the decomposition is based on frequency domain

X

, see more in Cooley–Tukey FFT algorithm.

Using the change of variables:

ni=pi+qiN/r

, where

pi=0,\ldots,(N/r)-1;qi=0,\ldots,r-1;

ki=rui+vi

, where

ui=0,\ldots,(N/r)-1;vi=0,\ldots,r-1;

where

i=1

or

2

, and the DFT equation can be written as:[6]

X(ru1+v1,ru2+v2)=\sum

N/r-1
p1=0
N/r-1
\sum
p2=0

\left[

r-1
\sum
q1=0
r-1
\sum
q2=0

x[p1+q1N/r,p2+q2N/r]

q1v1
W
r
q2v2
W
r

\right]

p1v1+p2v2
W
N
p1u1
W
N/r
p2u2
W
N/r

,

Other approaches

The split-radix FFT algorithm has been proved to be a useful method for 1-D DFT. And this method has been applied to the vector-radix FFT to obtain a split vector-radix FFT.[6] [7]

In conventional 2-D vector-radix algorithm, we decompose the indices

k1,k2

into 4 groups:

\begin{array}{lcl} X(2k1,2k2)&:&even-even\\ X(2k1,2k2+1)&:&even-odd\\ X(2k1+1,2k2)&:&odd-even\\ X(2k1+1,2k2+1)&:&odd-odd \end{array}

By the split vector-radix algorithm, the first three groups remain unchanged, the fourth odd-odd group is further decomposed into another four sub-groups, and seven groups in total:

\begin{array}{lcl} X(2k1,2k2)&:&even-even\\ X(2k1,2k2+1)&:&even-odd\\ X(2k1+1,2k2)&:&odd-even\\ X(4k1+1,4k2+1)&:&odd-odd\\ X(4k1+1,4k2+3)&:&odd-odd\\ X(4k1+3,4k2+1)&:&odd-odd\\ X(4k1+3,4k2+3)&:&odd-odd\end{array}

That means the fourth term in 2-D DIT radix-

(2 x 2)

equation,

S11(k1,k2)

k1+k2
W
N
becomes:[8]

A11(k1,k2)

k1+k2
W
N

+A13(k1,k2)

k1+3k2
W
N

+A31(k1,k2)

3k1+k2
W
N

+A33(k1,k2)

3(k1+k2)
W
N

,

where

Aij(k1,k2)=\sum

N/4-1
n1=0
N/4-1
\sum
n2=0

x[4n1+i,4n2+j]

n1k1
W
N/4
n2k2
W
N/4

The 2-D N by N DFT is then obtained by successive use of the above decomposition, up to the last stage.

It has been shown that the split vector radix algorithm has saved about 30% of the complex multiplications and about the same number of the complex additions for typical

1024 x 1024

array, compared with the vector-radix algorithm.[7]

Notes and References

  1. Book: Dudgeon. Dan. Russell. Mersereau. Multidimensional Digital Signal Processing. September 1983. Prentice Hall. 0136049591. 76.
  2. Rivard. G.. Direct fast Fourier transform of bivariate functions. IEEE Transactions on Acoustics, Speech, and Signal Processing. 25. 3. 250–252. 10.1109/TASSP.1977.1162951. 1977.
  3. Book: Harris. D.. McClellan. J.. Chan. D.. Schuessler. H.. ICASSP '77. IEEE International Conference on Acoustics, Speech, and Signal Processing . Vector radix fast Fourier transform . 2. 548–551. 10.1109/ICASSP.1977.1170349. 1977.
  4. Buijs. H.. Pomerleau. A.. Fournier. M.. Tam. W.. Implementation of a fast Fourier transform (FFT) for image processing applications. IEEE Transactions on Acoustics, Speech, and Signal Processing. Dec 1974. 22. 6. 420–424. 10.1109/TASSP.1974.1162620.
  5. Book: Badar. S.. Dandekar. D.. 2015 International Conference on Industrial Instrumentation and Control (ICIC) . High speed FFT processor design using radix −4 pipelined architecture . 1050–1055. 10.1109/IIC.2015.7150901. 2015. 978-1-4799-7165-7. 11093545 .
  6. Chan. S. C.. Ho. K. L.. Split vector-radix fast Fourier transform. IEEE Transactions on Signal Processing. 40. 8. 2029–2039. 10.1109/78.150004. 1992ITSP...40.2029C. 1992.
  7. Book: Pei. Soo-Chang. Wu. Ja-Lin. ICASSP '87. IEEE International Conference on Acoustics, Speech, and Signal Processing . Split vector radix 2D fast Fourier transform . 12. April 1987. 1987–1990. 10.1109/ICASSP.1987.1169345. 118173900 .
  8. Wu. H.. Paoloni. F.. On the two-dimensional vector split-radix FFT algorithm. IEEE Transactions on Acoustics, Speech, and Signal Processing. Aug 1989. 37. 8. 1302–1304. 10.1109/29.31283.