Hexagonal Fast Fourier Transform

Motivation

The Hexagonal Fast Fourier Transform (HFFT) is a variant of the Fast Fourier Transformation that is developed to utilize the advantages of hexagonal sampling. It has been proven that the hexagonal grid serves as an optimal sampling lattice for isotropically band-limited two-dimensional signals and has a sampling efficiency which is 13.4% more when compared to the sampling efficiency obtained from rectangular sampling. Several other advantages of hexagonal sampling include consistent connectivity, higher symmetry, greater angular resolution, and equidistant neighbouring pixels. Sometimes, more than one of these advantages compound together thereby increasing the efficiency by 50% in terms of computation and storage when compared to rectangular sampling. Despite of all these advantages of hexagonal sampling over the rectangular sampling, its application has been limited because of its inability to result in the orthogonal rows and columns that can be transformed independently as is done with rectangular samples. For this reason, an efficient FFT with a separable Fourier kernel has been developed by utilizing the Array Set Addressing (ASA) scheme and is termed as Hexagonal Fast Fourier Transform.

Preliminaries

Array Set Addressing

Figure 1. The figure illustrates representation of a hexagonally sampled data in the form of rectangular arrays using ASA coordinate system.

The ASA coordinate system has been proposed based on a simple fact that a hexagonal grid can be considered as a combination of two interleaved rectangular arrays. From this, it is easy to address each individual array using familiar integer-valued row and column indices. Now, to select one among these two rectangular arrays, a single binary coordinate is used. From this, a full address for any point in the hexagonal grid can be uniquely represented by three coordinates.

[math]\displaystyle{ (a,r,c) \in \{ 0,1 \} * n_1 * n_2 }[/math]

where the coordinates a,r and c represent the array, row and column respectively. The hexagonal lattice can be separated into rectangular arrays by considering alternate rows as one array and the remaining rows as the other. Figure 1 can be used to illustrate the process of Array Set Addressing, where alternate rows of the hexagonally sampled data are mapped to 0-array and 1-array rectangular arrays.

Hexagonal Discrete Fourier Transform (HDFT)

The Hexagonal DFT has been developed by Mersereau and it has been converted to an ASA representation by Rummelt. Let [math]\displaystyle{ x(a, r, c) }[/math] be a two-dimensional hexagonally sampled signal and let both arrays be of size [math]\displaystyle{ nXm }[/math]. Let, [math]\displaystyle{ X(b, s, d) }[/math] be the Fourier transform of x. The HDFT equation for the forward transform as shown in is given by

[math]\displaystyle{ X(b, s, d) = \sum_{a} \sum_{r} \sum_{c} x(a, r, c)E(.) }[/math]

where

[math]\displaystyle{ E(.) = exp[-j\pi(\frac{(a+2c)(b+2d)}{2m} + \frac{(a+2r)(b+2s)}{n})] }[/math]

Note that the above equation is separable and hence can be expressed as

[math]\displaystyle{ X(b, s, d) = f_{0}(b, s, d) + W(.) f_{1}(b, s, d) }[/math]

where

[math]\displaystyle{ W(.) = exp[-j\pi(\frac{b+2d}{2m} + \frac{b+2s}{n})] }[/math]

and

[math]\displaystyle{ g_{a}(b, r, d) = \sum_{c} x(a, r, c) exp(-j2\pi \frac{(c)(b+2d)}{2m} ) }[/math]

[math]\displaystyle{ f_{a}(b, s, d) = \sum_{r} g_{a}(b, r, d) exp(-j2\pi\frac{(r)(b+2s)}{n}) }[/math]

Existing Approaches

The linear transforms [math]\displaystyle{ g_{a} }[/math] and [math]\displaystyle{ f_{a} }[/math] are similar to the rectangular Fourier kernel where a linear transform is applied along each dimension of the 2-D rectangular data. Note that each of these equations, discussed above, is a combination of four rectangular arrays that serve as precursors to the HDFT. Two, out of those four rectangular [math]\displaystyle{ g_{a} }[/math] terms contribute to the sub-array of HFFT. Now, by switching the binary coordinate, we have four different forms of equations. In, the three out of those four expressions have been evaluated using the Non-Standard Transforms (NSTs) while one expression is computed using any correct and applicable FFT algorithm.

[math]\displaystyle{ g_{a}(0, r, d) = \sum_{c} x(a, r, c) exp(-j2\pi \frac{(c)(d)}{m}) }[/math]

[math]\displaystyle{ g_{a}(1, r, d) = \sum_{c} x(a, r, c) exp(-j2\pi \frac{(c)(2d+1)}{2m}) }[/math]

[math]\displaystyle{ f_{a}(0, s, d) = \sum_{r} g_{a}(a, r, d) exp(-j2\pi \frac{(r)(2s)}{n}) }[/math]

[math]\displaystyle{ f_{a}(1, s, d) = \sum_{r} g_{a}(a, r, d) exp(-j2\pi \frac{(r)(2s+1)}{n}) }[/math]

Looking at the second expression, [math]\displaystyle{ g_{a}(1,r,d) }[/math], we see that it is nothing more than a standard DFT with a constant offset along the rows of rectangular sub-arrays of a hexagonally-sampled image [math]\displaystyle{ x(a, r, c) }[/math]. This expression is nothing more than a circular rotation of the DFT. Note that the shift must happen in integer number of samples for the property to hold. This way, the function [math]\displaystyle{ g_{a} }[/math] can be computed using the standard DFT, in same number of operations, without introducing an NST.

If we have a look at 0-array [math]\displaystyle{ f_{a} }[/math], the expression will always be symmetric about half its spatial period. Because of this, it is enough to compute only half of it. We find that this expression is the standard DFT of the columns of [math]\displaystyle{ g_{a} }[/math], which is decimated by a factor of 2 and then is duplicated to span the space of r for the identical second period of the complex exponential. Mathematically,

[math]\displaystyle{ X_{even}[k] = \sum_{n=0}^{N-1} x[n]e^{-\tfrac{2j\pi}{N}2kn} }[/math]

[math]\displaystyle{ = \sum_{n=0}^{\tfrac{N}{2}-1} x[n]e^{-\tfrac{2j\pi}{N/2}kn} + \sum_{n=\tfrac{N}{2}}^{N-1} x[n]e^{-\tfrac{2j\pi}{N/2}kn} }[/math]

[math]\displaystyle{ = \sum_{n=0}^{\tfrac{N}{2}-1} x[n]e^{-\tfrac{2j\pi}{N/2}kn} + \sum_{n=0}^{\tfrac{N}{2}-1} x[n+\frac{N}{2}]e^{-\tfrac{2j\pi}{N/2}kn} }[/math]

[math]\displaystyle{ = \sum_{n=0}^{\tfrac{N}{2}-1} (x[n]+x[n+\tfrac{N}{2}])e^{-\tfrac{2j\pi}{N/2}kn} }[/math]

The expression for the 1-array [math]\displaystyle{ f_{a} }[/math] is equivalent to the 0-array expression with a shift of one sample. Hence, the 1-array expression can be expressed as columns of the DFT of [math]\displaystyle{ g_{a} }[/math] decimated by a factor of two, starting with second sample providing a constant offset needed by 1-array, and then doubled in space to span the range of s. Thus, the method developed by James B. Birdsong and Nicholas I. Rummelt in is able to successfully compute the Hexagonal Fast Fourier Transform (HFFT) using the standard FFT routines unlike the previous work in.