Hexagonal fast Fourier transform

The hexagonal fast Fourier transform (HFFT) is a tool in image and signal processing which uses fast Fourier transform (FFT) routines to compute the discrete Fourier transform (DFT) of images captured with hexagonal sampling. Hexagonal sampling's application is limited due to the lack of an efficient coordinate system. The existence of a separable Fourier kernel for a hexagonally sampled image allows the use of existing FFT routines to efficiently compute the DFT of such an image.

Preliminaries

Hexagonal Efficient Coordinate System (HECS)

Main article: Hexagonal Efficient Coordinate System

Representation of hexagonally sampled data as a pair of rectangular arrays using the HECS coordinate system

The Hexagonal Efficient Coordinate System (formerly known as Array Set Addressing (ASA)) was developed based on the fact that a hexagonal grid can be represented as a combination of two interleaved rectangular arrays. One can address each individual array by using integer-valued row and column indices, and the individual arrays can be distinguished by a single binary coordinate. Therefore, a full address for any point in the hexagonal grid can be uniquely represented by three coordinates:

$(a,r,c)\in \{0,1\}\times \mathbb {Z} \times \mathbb {Z}$

where the coordinates a, r and c represent the array, row and column respectively. The figure shows how the hexagonal grid is represented by two interleaved rectangular arrays in HECS coordinates.

Hexagonal discrete Fourier transform

The hexagonal discrete Fourier transform (HDFT) has been developed by Mersereau and it has been converted to an HECS representation by Rummelt. Let $x(a,r,c)$ be a two-dimensional hexagonally sampled signal and let both arrays be of size $n\times m$. Let, $X(b,s,d)$ be the Fourier transform of x. The HDFT equation for the forward transform as shown in is given by

$X(b,s,d)=\sum _{a}\sum _{r}\sum _{c}x(a,r,c)E(\cdot )$

where

$E(\cdot )=\exp \left[-j\pi \left({\frac {(a+2c)(b+2d)}{2m}}+{\frac {(a+2r)(b+2s)}{n}}\right)\right]$

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

$X(b,s,d)=f_{0}(b,s,d)+W(\cdot )f_{1}(b,s,d)$

where

$W(\cdot )=\exp \left[-j\pi \left({\frac {b+2d}{2m}}+{\frac {b+2s}{n}}\right)\right]$

and

$g_{a}(b,r,d)=\sum _{c}x(a,r,c)\exp \left(-j2\pi {\frac {(c)(b+2d)}{2m}}\right)$

$f_{a}(b,s,d)=\sum _{r}g_{a}(b,r,d)\exp \left(-j2\pi {\frac {(r)(b+2s)}{n}}\right)$

Hexagonal fast Fourier transform (HFFT)

The linear transforms $g_{a}$ and $f_{a}$ are similar to the rectangular Fourier kernel where a linear transform is applied along each dimension of the 2-D rectangular data. 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 $g_{a}$ terms contribute to the sub-array of HFFT. By switching the binary coordinate, we have four different forms of equations. Three out of those four expressions have been evaluated using "non-standard transforms (NSTs)" (shown below) while one expression is computed using any correct and applicable FFT algorithm.

$g_{a}(0,r,d)=\sum _{c}x(a,r,c)\exp \left(-j2\pi {\frac {(c)(d)}{m}}\right)$

$g_{a}(1,r,d)=\sum _{c}x(a,r,c)\exp \left(-j2\pi {\frac {(c)(2d+1)}{2m}}\right)$

$f_{a}(0,s,d)=\sum _{r}g_{a}(a,r,d)\exp \left(-j2\pi {\frac {(r)(2s)}{n}}\right)$

$f_{a}(1,s,d)=\sum _{r}g_{a}(a,r,d)\exp \left(-j2\pi {\frac {(r)(2s+1)}{n}}\right)$

The second expression, $g_{a}(1,r,d)$, is a standard discrete Fourier transform (DFT) with a constant offset along the rows of rectangular sub-arrays of a hexagonally-sampled image $x(a,r,c)$. This expression is nothing more than a circular rotation of the DFT. Note that the shift must happen in the integer number of samples for the property to hold. This way, the function $g_{a}$ can be computed using the standard DFT, in same number of operations, without introducing an NST.

Since 0-array $f_{a}$ will always be symmetric about half its spatial period, it is enough to compute only half of it. This expression is the standard DFT of the columns of $g_{a}$, 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,

${\begin{aligned}X_{\text{even}}[k]&=\sum _{n=0}^{N-1}x[n]e^{-{\tfrac {2j\pi }{N}}2kn}\\[5pt]&=\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}\\[5pt]&=\sum _{n=0}^{{\tfrac {N}{2}}-1}x[n]e^{-{\tfrac {2j\pi }{N/2}}kn}+\sum _{n=0}^{{\tfrac {N}{2}}-1}x\left[n+{\tfrac {N}{2}}\right]e^{-{\tfrac {2j\pi }{N/2}}kn}\\[5pt]&=\sum _{n=0}^{{\tfrac {N}{2}}-1}\left(x[n]+x\left[n+{\tfrac {N}{2}}\right]\right)e^{-{\tfrac {2j\pi }{N/2}}kn}\end{aligned}}$

The expression for the 1-array $f_{a}$ 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 $g_{a}$ 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 is able to successfully compute the HFFT using the standard FFT routines.