Hexagonal fast Fourier transform

<h2 id="preliminaries">Preliminaries</h2>
<h3>Hexagonal Efficient Coordinate System (HECS)</h3>
Main article: <a href="/facts/Hexagonal_Efficient_Coordinate_System_(HECS)/8anYKQvX">Hexagonal Efficient Coordinate System (HECS)</a>

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.<a class="footnote-ref" id="fnref:8" href="#fn:8">8</a> It is easy to address each individual array using familiar integer-valued row and column indices and the individual arrays are 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
        )
        ∈
        {
        0
        ,
        1
        }
        ×
        
          Z
        
        ×
        
          Z
        
      
    
    {\displaystyle (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.

<h3>Hexagonal discrete Fourier transform</h3>
The hexagonal discrete Fourier transform (HDFT) has been developed by Mersereau<a class="footnote-ref" id="fnref:9" href="#fn:9">9</a> and it has been converted to an HECS representation by Rummelt.<a class="footnote-ref" id="fnref:10" href="#fn:10">10</a> Let 
 
 
 
 x
 (
 a
 ,
 r
 ,
 c
 )
 
 
 {\displaystyle x(a,r,c)}
 
 be a two-dimensional hexagonally sampled signal and let both arrays be of size 
 
 
 
 n
 ×
 m
 
 
 {\displaystyle n\times m}
 
. Let, 
 
 
 
 X
 (
 b
 ,
 s
 ,
 d
 )
 
 
 {\displaystyle X(b,s,d)}
 
 be the <a href="/facts/Fourier_transform/XXhKJtc8">Fourier transform</a> of x. The HDFT equation for the forward transform as shown in <a class="footnote-ref" id="fnref:11" href="#fn:11">11</a> is given by

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

where

E
        (
        ⋅
        )
        =
        exp
        ⁡
        
          [
          
            −
            j
            π
            
              (
              
                
                  
                    
                      (
                      a
                      +
                      2
                      c
                      )
                      (
                      b
                      +
                      2
                      d
                      )
                    
                    
                      2
                      m
                    
                  
                
                +
                
                  
                    
                      (
                      a
                      +
                      2
                      r
                      )
                      (
                      b
                      +
                      2
                      s
                      )
                    
                    n
                  
                
              
              )
            
          
          ]
        
      
    
    {\displaystyle 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
        (
        ⋅
        )
        
          f
          
            1
          
        
        (
        b
        ,
        s
        ,
        d
        )
      
    
    {\displaystyle X(b,s,d)=f_{0}(b,s,d)+W(\cdot )f_{1}(b,s,d)}

where

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

and

g
          
            a
          
        
        (
        b
        ,
        r
        ,
        d
        )
        =
        
          ∑
          
            c
          
        
        x
        (
        a
        ,
        r
        ,
        c
        )
        exp
        ⁡
        
          (
          
            −
            j
            2
            π
            
              
                
                  (
                  c
                  )
                  (
                  b
                  +
                  2
                  d
                  )
                
                
                  2
                  m
                
              
            
          
          )
        
      
    
    {\displaystyle 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
        )
        =
        
          ∑
          
            r
          
        
        
          g
          
            a
          
        
        (
        b
        ,
        r
        ,
        d
        )
        exp
        ⁡
        
          (
          
            −
            j
            2
            π
            
              
                
                  (
                  r
                  )
                  (
                  b
                  +
                  2
                  s
                  )
                
                n
              
            
          
          )
        
      
    
    {\displaystyle f_{a}(b,s,d)=\sum _{r}g_{a}(b,r,d)\exp \left(-j2\pi {\frac {(r)(b+2s)}{n}}\right)}

<h2 id="hexagonal-fast-fourier-transform-hfft">Hexagonal fast Fourier transform (HFFT)</h2>
The linear transforms 
 
 
 
 
 g
 
 a
 
 
 
 
 {\displaystyle g_{a}}
 
 and 
 
 
 
 
 f
 
 a
 
 
 
 
 {\displaystyle f_{a}}
 
 are similar to the rectangular Fourier kernel where a linear transform is applied along each dimension of the 2-D rectangular data.<a class="footnote-ref" id="fnref:12" href="#fn:12">12</a> 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 
 
 
 
 
 g
 
 a
 
 
 
 
 {\displaystyle g_{a}}
 
 terms contribute to the sub-array of HFFT. Now, by switching the binary coordinate, we have four different forms of equations. In,<a class="footnote-ref" id="fnref:13" href="#fn:13">13</a> the three out of those four expressions have been evaluated using what the author called "non-standard transforms (NSTs)" (shown below) while one expression is computed using any correct and applicable FFT algorithm.

g
          
            a
          
        
        (
        0
        ,
        r
        ,
        d
        )
        =
        
          ∑
          
            c
          
        
        x
        (
        a
        ,
        r
        ,
        c
        )
        exp
        ⁡
        
          (
          
            −
            j
            2
            π
            
              
                
                  (
                  c
                  )
                  (
                  d
                  )
                
                m
              
            
          
          )
        
      
    
    {\displaystyle 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
        )
        =
        
          ∑
          
            c
          
        
        x
        (
        a
        ,
        r
        ,
        c
        )
        exp
        ⁡
        
          (
          
            −
            j
            2
            π
            
              
                
                  (
                  c
                  )
                  (
                  2
                  d
                  +
                  1
                  )
                
                
                  2
                  m
                
              
            
          
          )
        
      
    
    {\displaystyle 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
        )
        =
        
          ∑
          
            r
          
        
        
          g
          
            a
          
        
        (
        a
        ,
        r
        ,
        d
        )
        exp
        ⁡
        
          (
          
            −
            j
            2
            π
            
              
                
                  (
                  r
                  )
                  (
                  2
                  s
                  )
                
                n
              
            
          
          )
        
      
    
    {\displaystyle 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
        )
        =
        
          ∑
          
            r
          
        
        
          g
          
            a
          
        
        (
        a
        ,
        r
        ,
        d
        )
        exp
        ⁡
        
          (
          
            −
            j
            2
            π
            
              
                
                  (
                  r
                  )
                  (
                  2
                  s
                  +
                  1
                  )
                
                n
              
            
          
          )
        
      
    
    {\displaystyle f_{a}(1,s,d)=\sum _{r}g_{a}(a,r,d)\exp \left(-j2\pi {\frac {(r)(2s+1)}{n}}\right)}

Looking at the second expression, 
 
 
 
 
 g
 
 a
 
 
 (
 1
 ,
 r
 ,
 d
 )
 
 
 {\displaystyle g_{a}(1,r,d)}
 
, we see that it is nothing more than a standard <a href="/facts/Discrete_Fourier_transform/uK9Tjsbl">discrete Fourier transform</a> (DFT) with a constant offset along the rows of rectangular sub-arrays of a hexagonally-sampled image 
 
 
 
 x
 (
 a
 ,
 r
 ,
 c
 )
 
 
 {\displaystyle x(a,r,c)}
 
.<a class="footnote-ref" id="fnref:14" href="#fn:14">14</a> This expression is nothing more than a <a href="/facts/Permutation/JSw9ukmv">circular rotation</a> of the DFT. Note that the shift must happen in the <a href="/facts/Integer/PUwwrawx">integer number</a> of samples for the property to hold. This way, the function 
 
 
 
 
 g
 
 a
 
 
 
 
 {\displaystyle g_{a}}
 
 can be computed using the standard DFT, in same number of operations, without introducing an NST.
If we have a look at 0-array 
 
 
 
 
 f
 
 a
 
 
 
 
 {\displaystyle f_{a}}
 
, the expression will always be symmetric about half its <a href="/facts/Wavelength/kcJU0ymz">spatial period</a>. 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 
 
 
 
 
 g
 
 a
 
 
 
 
 {\displaystyle 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.<a class="footnote-ref" id="fnref:15" href="#fn:15">15</a> Mathematically,

X
                  
                    even
                  
                
                [
                k
                ]
              
              
                
                =
                
                  ∑
                  
                    n
                    =
                    0
                  
                  
                    N
                    −
                    1
                  
                
                x
                [
                n
                ]
                
                  e
                  
                    −
                    
                      
                        
                          
                            2
                            j
                            π
                          
                          N
                        
                      
                    
                    2
                    k
                    n
                  
                
              
            
            
              
              
                
                =
                
                  ∑
                  
                    n
                    =
                    0
                  
                  
                    
                      
                        
                          N
                          2
                        
                      
                    
                    −
                    1
                  
                
                x
                [
                n
                ]
                
                  e
                  
                    −
                    
                      
                        
                          
                            2
                            j
                            π
                          
                          
                            N
                            
                              /
                            
                            2
                          
                        
                      
                    
                    k
                    n
                  
                
                +
                
                  ∑
                  
                    n
                    =
                    
                      
                        
                          N
                          2
                        
                      
                    
                  
                  
                    N
                    −
                    1
                  
                
                x
                [
                n
                ]
                
                  e
                  
                    −
                    
                      
                        
                          
                            2
                            j
                            π
                          
                          
                            N
                            
                              /
                            
                            2
                          
                        
                      
                    
                    k
                    n
                  
                
              
            
            
              
              
                
                =
                
                  ∑
                  
                    n
                    =
                    0
                  
                  
                    
                      
                        
                          N
                          2
                        
                      
                    
                    −
                    1
                  
                
                x
                [
                n
                ]
                
                  e
                  
                    −
                    
                      
                        
                          
                            2
                            j
                            π
                          
                          
                            N
                            
                              /
                            
                            2
                          
                        
                      
                    
                    k
                    n
                  
                
                +
                
                  ∑
                  
                    n
                    =
                    0
                  
                  
                    
                      
                        
                          N
                          2
                        
                      
                    
                    −
                    1
                  
                
                x
                
                  [
                  
                    n
                    +
                    
                      
                        
                          N
                          2
                        
                      
                    
                  
                  ]
                
                
                  e
                  
                    −
                    
                      
                        
                          
                            2
                            j
                            π
                          
                          
                            N
                            
                              /
                            
                            2
                          
                        
                      
                    
                    k
                    n
                  
                
              
            
            
              
              
                
                =
                
                  ∑
                  
                    n
                    =
                    0
                  
                  
                    
                      
                        
                          N
                          2
                        
                      
                    
                    −
                    1
                  
                
                
                  (
                  
                    x
                    [
                    n
                    ]
                    +
                    x
                    
                      [
                      
                        n
                        +
                        
                          
                            
                              N
                              2
                            
                          
                        
                      
                      ]
                    
                  
                  )
                
                
                  e
                  
                    −
                    
                      
                        
                          
                            2
                            j
                            π
                          
                          
                            N
                            
                              /
                            
                            2
                          
                        
                      
                    
                    k
                    n
                  
                
              
            
          
        
      
    
    {\displaystyle {\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
 
 
 
 
 {\displaystyle 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
 
 
 
 
 {\displaystyle 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 in <a class="footnote-ref" id="fnref:16" href="#fn:16">16</a> is able to successfully compute the HFFT using the standard FFT routines unlike the previous work in.<a class="footnote-ref" id="fnref:17" href="#fn:17">17</a>

<h2 id="references">References</h2>

<ol>
<li id="fn:1">James B. Birdsong, Nicholas I. Rummelt, "The Hexagonal Fast Fourier Transform", 2016 IEEE International Conference on Image Processing (ICIP), pp. 1809–1812, doi:10.1109/ICIP.2016.7532670 <a href="/wiki/Doi_(identifier)" target="_blank">/wiki/Doi_(identifier)</a> <a href="#fnref:1" class="footnote-back-ref">↩</a></li>
<li id="fn:2">D. P. Petersen and D. Middleton, Dec. 1962, "Sampling and reconstruction of wave-number-limited functions in n-dimensional Euclidean spaces", Inf. Control, vol. 5, no. 4, pp. 279–323 <a href="#fnref:2" class="footnote-back-ref">↩</a></li>
<li id="fn:3">R. M. Mersereau, June 1979, "The Processing of Hexagonally Sampled Two-Dimensional Signals", Proceedings of the IEEE, vol. 67, no. 6, pp. 930–949 <a href="#fnref:3" class="footnote-back-ref">↩</a></li>
<li id="fn:4">X. He and W. Jia, 2005, "Hexagonal structure for intelligent vision", in Proc. 1st Int. Conf. Information and Communication Technologies, pp. 52–64 <a href="#fnref:4" class="footnote-back-ref">↩</a></li>
<li id="fn:5">W. E. Snyder, 1999, H. Qi, and W. Sander, "A coordinate system for hexagonal pixels", in Proc. SPIE Medical Imaging: Image Processing, vol. 3661, pp. 716–727 <a href="#fnref:5" class="footnote-back-ref">↩</a></li>
<li id="fn:6">R. M. Mersereau, June 1979, "The Processing of Hexagonally Sampled Two-Dimensional Signals", Proceedings of the IEEE, vol. 67, no. 6, pp. 930–949 <a href="#fnref:6" class="footnote-back-ref">↩</a></li>
<li id="fn:7">Nicholas I. Rummelt and Joseph N. Wilson "Array set addressing: enabling technology for the efficient processing of hexagonally sampled imagery," Journal of Electronic Imaging 20(2), 023012 (1 April 2011). https://doi.org/10.1117/1.3589306 <a href="https://doi.org/10.1117/1.3589306" target="_blank">https://doi.org/10.1117/1.3589306</a> <a href="#fnref:7" class="footnote-back-ref">↩</a></li>
<li id="fn:8">Nicholas I. Rummelt, 2010, Array Set Addressing: Enabling Efficient Hexagonally Sampled Image Processing, Ph.D. thesis, University of Florida <a href="#fnref:8" class="footnote-back-ref">↩</a></li>
<li id="fn:9">R. M. Mersereau, June 1979, "The Processing of Hexagonally Sampled Two-Dimensional Signals", Proceedings of the IEEE, vol. 67, no. 6, pp. 930–949 <a href="#fnref:9" class="footnote-back-ref">↩</a></li>
<li id="fn:10">Nicholas I. Rummelt, 2010, Array Set Addressing: Enabling Efficient Hexagonally Sampled Image Processing, Ph.D. thesis, University of Florida <a href="#fnref:10" class="footnote-back-ref">↩</a></li>
<li id="fn:11">Nicholas I. Rummelt, 2010, Array Set Addressing: Enabling Efficient Hexagonally Sampled Image Processing, Ph.D. thesis, University of Florida <a href="#fnref:11" class="footnote-back-ref">↩</a></li>
<li id="fn:12">James B. Birdsong, Nicholas I. Rummelt, "The Hexagonal Fast Fourier Transform", 2016 IEEE International Conference on Image Processing (ICIP), pp. 1809–1812, doi:10.1109/ICIP.2016.7532670 <a href="/wiki/Doi_(identifier)" target="_blank">/wiki/Doi_(identifier)</a> <a href="#fnref:12" class="footnote-back-ref">↩</a></li>
<li id="fn:13">Nicholas I. Rummelt, 2010, Array Set Addressing: Enabling Efficient Hexagonally Sampled Image Processing, Ph.D. thesis, University of Florida <a href="#fnref:13" class="footnote-back-ref">↩</a></li>
<li id="fn:14">James B. Birdsong, Nicholas I. Rummelt, "The Hexagonal Fast Fourier Transform", 2016 IEEE International Conference on Image Processing (ICIP), pp. 1809–1812, doi:10.1109/ICIP.2016.7532670 <a href="/wiki/Doi_(identifier)" target="_blank">/wiki/Doi_(identifier)</a> <a href="#fnref:14" class="footnote-back-ref">↩</a></li>
<li id="fn:15">James B. Birdsong, Nicholas I. Rummelt, "The Hexagonal Fast Fourier Transform", 2016 IEEE International Conference on Image Processing (ICIP), pp. 1809–1812, doi:10.1109/ICIP.2016.7532670 <a href="/wiki/Doi_(identifier)" target="_blank">/wiki/Doi_(identifier)</a> <a href="#fnref:15" class="footnote-back-ref">↩</a></li>
<li id="fn:16">James B. Birdsong, Nicholas I. Rummelt, "The Hexagonal Fast Fourier Transform", 2016 IEEE International Conference on Image Processing (ICIP), pp. 1809–1812, doi:10.1109/ICIP.2016.7532670 <a href="/wiki/Doi_(identifier)" target="_blank">/wiki/Doi_(identifier)</a> <a href="#fnref:16" class="footnote-back-ref">↩</a></li>
<li id="fn:17">Nicholas I. Rummelt, 2010, Array Set Addressing: Enabling Efficient Hexagonally Sampled Image Processing, Ph.D. thesis, University of Florida <a href="#fnref:17" class="footnote-back-ref">↩</a></li>
</ol>

Hexagonal fast Fourier transform open-in-new

Hexagonal fast Fourier transform