Discrete two dimensional Fourier transform in polar coordinates part II: numerical computation and approximation of the continuous transform

The theory of the continuous two-dimensional (2D) Fourier Transform in polar coordinates has been recently developed but no discrete counterpart exists to date. In the first part of this two-paper series, we proposed and evaluated the theory of the 2D Discrete Fourier Transform (DFT) in polar coordinates. The theory of the actual manipulated quantities was shown, including the standard set of shift, modulation, multiplication, and convolution rules. In this second part of the series, we address the computational aspects of the 2D DFT in polar coordinates. Specifically, we demonstrate how the decomposition of the 2D DFT as a DFT, Discrete Hankel Transform and inverse DFT sequence can be exploited for coding. We also demonstrate how the proposed 2D DFT can be used to approximate the continuous forward and inverse Fourier transform in polar coordinates in the same manner that the 1D DFT can be used to approximate its continuous counterpart.


INTRODUCTION
The Fourier Transform (FT) is a powerful analytical tool and has proved to be invaluable in many disciplines such as physics, mathematics and engineering. The development of the Fast Fourier Transform (FFT) algorithm (Cooley & Tukey, 1965), which computes the Discrete Fourier Transform (DFT) with a fast algorithm, firmly established the FT as a practical tool in diverse areas, most notably signal and image processing.
In two dimensions, the FFT can still be used to compute the DFT in Cartesian coordinates. However, in many applications such as photoacoustics (Xu, Feng & Wang, 2002) and tomography (Scott et al., 2012;Fahimian et al., 2013;Lee et al., 2008;Lewitt & Matej, 2003), it is often necessary to compute the FT in polar coordinates. Moreover, for functions that are naturally described in polar coordinates, a discrete version of the 2D FT in polar coordinates is needed. There have been some attempts to calculate the FT in polar coordinates, most notably through the Hankel transform, since the zeroth order Hankel transform is known to be a 2D FT in polar coordinates for rotationally symmetric functions. However, prior work has focused on numerically approximating the continuous transform. This stands in contrast to the FT, where the DFT can stand alone as an orthogonal transform, independent of the existence of its continuous counterpart.
The idea of a Polar FT has been previously investigated, where the spatial function is in Cartesian coordinates but its FT is computed in polar coordinates (Averbuch et al., 2006;Abbas, Sun & Foroosh, 2017;Fenn, Kunis & Potts, 2007). FTs have been proposed for non-equispaced data, referred to as Unequally Spaced FFT (USFFT) or Non-Uniform FFT (NUFFT) (Dutt & Rokhlin, 1993;Fourmont, 2003;Dutt & Rokhlin, 1995;Potts, Steidl & Tasche, 2001;Fessler & Sutton, 2003). A recent book gives a unified treatment of these topics (Plonka et al., 2018). Previous work has also considered the implications of using a polar grid (Stark, 1979;Stark & Wengrovitz, 1983). Although the above references demonstrate that the computation of a discrete 2D FT on a polar grid has previously been considered in the literature, there is to date no discrete 2D FT in polar coordinates that exists as a transform in its own right, with its own set of rules of the actual manipulated quantities.
In part I of this two part series, we proposed an independent discrete 2D FT in polar coordinates, which has been defined to be discrete from first principles . For a discrete transform, the values of the transform are only given as entries in a vector or matrix and the transform manipulates a set of discrete values. To quote Bracewell (Bracewell, 1999), "we often think of this as though an underlying function of a continuous variable really exists and we are approximating it. From an operational viewpoint, however, it is irrelevant to talk about the existence of values other than those given and those computed (the input and output). Therefore, it is desirable to have a mathematical theory of the actual quantities manipulated". Hence, in our previous paper , standard operational 'rules' of shift, modulation and convolution rules for this 2D DFT in polar coordinates were demonstrated. The operational rules were demonstrated via the key properties of the proposed discrete kernel of the transform. However, using the discrete kernel may not be the most effective way to compute the transform. Furthermore, while the 2D DFT in polar coordinates was demonstrated to have properties and rules as a standalone transform independent of its relationship to any continuous transform, an obvious application of the proposed discrete transform is to approximate its continuous counterpart.
Hence, the goal of this second part of this two-part paper series is to propose computational approaches to the computation of the previously proposed 2D DFT in polar coordinates and also to validate its effectiveness to approximate the continuous 2D FT in polar coordinates.
The outline of the paper is as follows. "Definition of the Discrete 2D FT in Polar Coordinates" states the proposed definition of the discrete 2D FT in polar coordinates. The motivation of this definition and the transform rules (multiplication, convolution, shift etc.) are given in the first part of this two-part paper. The transform exists in its own right and manipulates discrete quantities that do not necessarily stem from sampling an underlying continuous quantity. Nevertheless, the motivation for the definition of the transform is based on an implied underlying discretization scheme. "Discrete Transform to approximate the continuous transform" introduces the implied underlying discretization scheme where we show the connection between discrete samples of the continuous functions and the discrete transform, should it be desirable to interpret the transform in this manner. Here, the connection between using the proposed 2D DFT and sampled vales of the continuous functions is explained. The proposed 2D DFT was motivated by a specific sampling scheme (Discrete Transform to approximate the continuous transform) which can be plotted and analyzed for "grid coverage"-how much of the 2D plane is covered and at which density. Thus, "Discretization Points and Sampling Grid" analyzes the proposed discretization points and their implication on the sampling grid for density and coverage of the grid. The insights gained from this section will be useful in interpreting the results of approximating the continuous transform with the discrete transform. "Numerical Computation of the Transform" introduces numerical computation schemes whereby the interpretation of the proposed 2D transform as a sequence of 1D DFT, 1D Discrete Hankel Transform (DHT) and 1D inverse DFT (IDFT) is exploited. "Numerical Evaluation of the 2D DFT in Polar Coordinates to Approximate the Continuous FT" then investigates the ability of the proposed 2D DFT to approximate the continuous transform in terms of precision and accuracy. Three test functions for which closed-form continuous transforms are known are analyzed. Finally, "Summary and Conclusion" summarizes and concludes the paper.

DEFINITION OF THE DISCRETE 2D FT IN POLAR COORDINATES
The 2D-DFT in polar coordinates has been defined in the first part of this two-paper series as the discrete transform that takes the matrix (or double-subscripted series) f pk to the matrix (double-subscripted series) F ql such that f pk ! F qm is given by (1) where p; k; q; m; n, N 1 , and N 2 are integers such that ÀM n M, where 2M þ 1 ¼ N 2 1 m; k; N 1 À 1 and ÀM p; q M. Unless otherwise stated, in the remainder of the paper it shall be assumed that p; k; q; m; n, N 1 , and N 2 are within these stated ranges.
Similarly, for the inverse transform we propose In Eqs. (1) and (2), E AE qm;pk are the kernels of the transformation. These can be chosen as the "non-symmetric" form given by Here, J n z ð Þ is the nth order Bessel function of the first kind and j nk denotes the kth zero of the nth Bessel function. The subscript (+ or −) indicated the sign on the i AE and on the exponent containing the p variable; the q variable exponent then takes the opposite sign. From a matrix point of view, both f pk and F ql are N 2 Â N 1 À 1 ð Þsized matrices. The form of the kernel in Eq. (3) arises naturally from discretization of the continuous transform, but does not lead to the expected Parseval relationship. A possible symmetric kernel is discussed in the first part of this two-part paper and Parseval relationships are discussed further there .

DISCRETE TRANSFORM TO APPROXIMATE THE CONTINUOUS TRANSFORM
In this section, relationships between discretely sampled values of the function and its continuous 2D FT are presented in the case of a space-limited or band-limited function. These relationships were derived in the first part of the paper and are repeated here to demonstrate how they form the basis for the using the discrete transform to approximate the continuous transform at specified sampling points.

Space-limited functions
Consider a function in the space domain f ðr; uÞ which is space limited to r 2 0; R ½ . This implies that the function is zero outside of the circle bounded by r 2 0; R ½ . An approximate relationship between sampled values of the continuous function and sampled values of its continuous forward 2D transform F r; c ð Þ has been derived in the first part of the two-part paper as Similarly, an approximate relationship between sampled values of the continuous forward transform F r; c ð Þ and sampled values of the continuous original function f ðr; uÞ was shown to be given by In Eqs. (4) and (5), f(r, θ) is the original function in 2D space and Fðr; cÞ is the 2D FT of the function in polar coordinates.
To evaluate if the 2D DFT as proposed in Eqs. (1) and (2) can be used to approximate sampled values of f(r, θ) and Fðr; cÞ, the process is as follows. For the forward transform, we start with the continuous f(r, θ), evaluate it at the sampling points and then assign this value to f pk via Then, F qm is calculated from the 2D DFT scaled by 2pR 2 , Eq.
(1), that is The factor of 2pR 2 is necessary so that the evaluation in Eq. (7) matches the expression in Eq. (4). To evaluate if the proposed 2D DFT can be used to approximate the continuous transform, the question becomes how well F qm calculated from the 2D DFT in Eq. (7) approximates F j qm -the values of the continuous 2D FT evaluated on the sampling grid.
To evaluate the inverse 2D DFT, the process is similar. We start with the continuous Fðr; cÞ, evaluate it at the sampling points and assign this value to F qm via Now, f pk is calculated from a scaled version of the inverse 2D DFT, Eq.
(2) that is To evaluate if the proposed transform can approximate the continuous transform, the question becomes how well f pk calculated from Eq. (9) approximates f j pk R j pN1 ; 2pp

Band-limited functions
The process for band-limited functions follows the same process as outlined in the previous section, with the exception that the sampling points and scaling factors are slightly different as they are now given in terms of the band limit rather than the space limit. Now consider functions in the frequency domain F q; c ð Þ with an effective band limit r 2 0; W r Â Ã . That is, we suppose that the 2D FT F r; c ð Þ of f ðr; uÞ is band-limited, meaning that F r; c ð Þ is zero for r ! W r ¼ 2pW. The variable W r is written in this form since W would typically be quoted in units of Hz (cycles per second) if using temporal units or cycles per meter if using spatial units. Therefore, the multiplication by 2p ensures that the final units are in s À1 or m À1 . The approximate relationship between sampled values of the continuous 2D FT F r; c ð Þ and sampled values of the original continuous function f r; u ð Þ was derived in the first part of the paper and is given by N 2 (10) Similarly, the inverse relationship between sampled values of F r; c ð Þ and sampled values of f ðr; uÞ was shown to be given by The relationships in Eqs. (10) and (11) give relationships between the sampled values of the original function and sampled values of its 2D FT.
To evaluate the forward 2D DFT, we start with f r; u ð Þ, evaluate it at the (bandlimited specific) sampling points and assign this value to f pk via Then, F qm is calculated from the discrete transform scaled by 2p To evaluate if the proposed 2D DFT can be used to approximate the continuous transform, the question is how well F qm calculated from Eq. (13) approximates  (11), were the motivating definition of a 2D DFT in polar coordinates, defined in the first part of this two-part paper. In the context of this second part of the two-part paper, they are also the relationships that permit the use of the discrete transform to approximate the continuous transform at the specified sampling points. They are also the relationships that permit the examination of whether the discrete quantities f pk and F qm calculated via the proposed 2D DFT are in fact reasonable approximations to the sampled values of the continuous functions, as stated in the objectives of the paper.

DISCRETIZATION POINTS AND SAMPLING GRID
The transforms defined in Eqs. (1) and (2) can be applied to any matrix f pk to yield its forward transform F qm , which can then be transformed backwards by using the inverse transform. However, if these same discrete transforms are to be used for the purpose of approximating a continuous 2D FT, then these transforms need to be applied to the specific sampled values of the continuous functions in both space and frequency domains, as shown in Eqs. (6), (8) and (12). The relationships in Eqs. (4) and (10) define the sampling points that need to be used and it is noted that the points are defined differently based on whether we start with the assumption of a space or band limited function. These specific sampling points imply a specific sampling grid for the function. In this section, the sampling grid (its coverage and density in 2D) is analyzed.

Sampling points
For a space-limited function, we assume that the original function of interest is defined over continuous r; u ð Þ space where 0 r R and 0 u 2p. The discrete sampling spaces used for radial and angular sampling points in regularr space r; u ð Þ andṽ frequency r; c ð Þ space are defined as and For a band limited function, the function is assumed band-limited to 0 r W r , 0 c 2p. The sampling space used for radial and angular sampling points in regularṽ frequency space r; c ð Þ andr space r; u ð Þ for a bandlimited function is defined as and Clearly, the density of the sampling points depends on the numbers of points chosen, that is on N 1 and N 2 . Also clear is the fact that the grid is not equispaced in the radial variable. The sampling grid for a space-limited function are plotted below to enable visualization. In the first instance, the polar grids are plotted for the case R ¼ 1, N 1 ¼ 16 and N 2 ¼ 15. These are shown in space (r space) and frequency (ρ space) in Figs. 1 and 2 respectively. It should be noted that although we refer the grids in this article as polar grids, they are not true polar grids in the sense of equispaced sampling in the radial and angular coordinates.
Clearly, the grids in Figs. 1 and 2 are fairly sparse, but the low values of N 2 and N 1 have been chosen so that the structure of the sampling points can be easily seen. It can be observed that there is a hole at the center area in both domains which is caused by the special sampling points. For higher values of the N 2 and N 1 , the grid becomes fairly dense, obtaining good coverage of both spaces, but details are harder to observe. To demonstrate, the polar grids are plotted for the case R = 1, N 1 ¼ 96 and N 2 ¼ 95. These are shown in Figs. 3 and 4 respectively. From Figs. 3 and 4, by choosing higher values of N 1 and N 2 , the sampling grid becomes denser, however there is still a gap in the center area. The sampling grids for band-limited functions are not plotted here since the sample grid for a band-limited function has the same shape as with space limited function but the domains are reversed.

Sample grid analysis
From part I of the paper, it was shown that the 2D-FT can be interpreted as a DFT in the angular direction, a DHT in the radial direction and then an IDFT in the angular direction. Hence, the sample size in the angular direction could have been decided by the Nyquist sampling theorem (Shannon, 1984), which states that where f s is the sample frequency and f max is the highest frequency or band limit. In the radial direction, the necessary relationship for the DHT is given by Baddour & Chouinard (2015) where W r is the effective band-limit, R is the effective space limit and j nN is the Nth zero of J n r ð Þ. For the 2D FT, since ÀM p M, the order of the Bessel zero ranges from ÀM to M, the required relationship becomes The relationships j nN ¼ j ÀnN and j 0N 1 < j AE1N 1 < j AE2N 1 < … < j AEMN 1 are valid (Lozier, 2003), hence Eq. (20) can be written as It is pointed out in Baddour (2019) and Guizar-Sicairos & Gutiérrez-Vega (2004) that the zeros of J n z ð Þ are almost evenly spaced at intervals of p and that the spacing becomes exactly p in the limit as z ! 1. The reader unfamiliar with Bessel functions is directed to references (Bracewell, 1999;Lozier, 2003). In fact, it is shown in Dutt & Rokhlin (1993) that a simple asymptotic form for the Bessel function is given by Therefore, an approximation to the Bessel zero, j nk is given by Hence, Eq. (21) can be written to choose N 1 approximately as where the reader is reminded that the units of W is m −1 (the space equivalent of Hz). N 1 =R is the spatial sampling frequency and we see that Eq. (24) effectively makes the same statement as Eq. (18), as it should. Intuitively, more sample points lead to more information captured, which gives an expectation that increasing N 1 or N 2 individually will give a better sampling grid coverage. However, it can be seen from Figs. 1-4 that there is a gap in the center of the sample grid. From Eqs. (14) and (15), the area of the gap in the center is related to the ranges of p and k, that is N 2 and N 1 . In the sections below, it is assumed that the sampling theorems are already satisfied (that is, an appropriate space and band limit is selected) and the relationship between N 2 , N 1 and the size of the gap will be discussed.

Space-limited function
In this section, it is assumed that the function is a space limited function, defined in r 2 ½0; R. The sampling points are defined as Eq. (14) in the space domain and Eq. (15) in the frequency domain. In the following, a relationship between N 2 , N 1 and the area of the gap in both domains is discussed.

Sample grid in the space domain
In the space domain, the effective limit in the space domain, R, is fixed. To analyze how the values of N 2 and N 1 affect the coverage of the grid in space domain, consider the following definition of 'grid coverage' where r denotes the average radius of the gap (the hole in the middle of the grid). A r as defined in Eq. (25) is a measure of the "grid coverage" since it gives a percentage of how much of the original space limited domain area is captured by the discrete grid. For example, if the average radius of the center gap is zero, then A r would be 100%, that is, complete coverage. Based on the observation of Figs. 1 and 3, the relationship r 01 < r AE11 < r AE21 < r AEM1 is valid. Therefore, from Eq. (14), the average radius of the gap is given by Hence, Eq. (25) for grid coverage can be written as Table 1 shows the different values of grid coverage A r as the values of N 1 and N 2 are changed. From Table 1, it can be seen that increasing N 1 (sample size in the radial direction) tends to increase the grid coverage. Since the effective space limit R is fixed, from Eq. (21) it follows that increasing N 1 actually increases the effective band limit. However, increasing N 2 (sample size in angular direction) will result in a bigger gap in the center of the grid, which then decreases the coverage.

Sample grid in the frequency domain
Similarly, coverage of the grid in the frequency domain is defined as where r denotes the average radius of the gap. Since Then, it follows that Eq. (28) for frequency domain grid coverage can be written as From Eq. (30), it can be observed that the sample grid coverage in the frequency domain is affected by R, W r and M. Since N 2 ¼ 2M þ 1, in order to get a better grid coverage with a fixed W r , R and N 2 can be adjusted. Table 2 shows the grid coverage A r for different values of R and N 2 .
From Table 2, the conclusion for the frequency domain is that when the effective band limit is fixed, increasing R (effective space limit) tends to increase the coverage in the frequency domain, while increasing N 2 (sample size in the angular direction) decreases the coverage. However, from Eq. (21) it should be noted that to satisfy the sampling theorem, increasing R with fixed W r requires an increase in N 1 at the same time.

Band-limited function
In this section, we suppose that the function is an effectively band limited function, defined on r 2 ½0; W p . The sampling points are defined as in Eq. (16) in the space domain and as in the frequency domain. In this subsection, the relationship between N 2 , N 1 and the area of the gap in both domains is discussed.

Sampling grid in the space domain
The same definition of grid coverage in the space domain will be used as in Eq. (25). Since the sampling points of a band-limited function are given by Eqs. (16) and (17), the average radius of the gap can be defined as Therefore, the coverage of the grid in space domain can be written as It can be observed that the grid coverage in the space domain of a band-limited function is the same as the grid coverage in the frequency domain of space limited function.

Sample grid in frequency domain
The coverage of the grid in the frequency domain of a band limited function is defined by Eq. (28). With sampling points defined in Eq. (17), the average radius of the gap can be defined as The coverage of the grid in frequency domain can be written as It can be observed that the grid coverage in the frequency domain of a band-limited function is the same as the grid coverage in the space domain of a space limited function.

Conclusion
Based on the discussion above, the following conclusions can be made: 1. Increasing N 2 (angular direction) tends to decrease the sampling grid coverage in both domains. Increasing N 1 (radial direction) tends to increase the sampling coverage in the space domain for a space-limited function and in the frequency domain for a frequency-limited function. So, if a signal changes sharply in the angular direction such that large values of N 2 are needed, a large value of N 1 is also needed to compensate for the effect of increasing N 2 on the grid coverage.
2. For a space-limited function, if there is a lot of energy at the origin in the space domain, a larger value of N 1 will be required to ensure that the sampling grid gets as close to the origin as possible in the space domain. If the function has a lot of energy at the origin in the frequency domain, a large value for both N 1 and R will be required to ensure adequate grid coverage.
3. For a band-limited function, if there is a lot of energy at the origin in the frequency domain, a large value of N 1 will be needed to ensure that the sample grid gets as close to the origin as possible in the frequency domain. If the function has a lot of energy at the origin in the space domain, large values for both N 1 and W r are required.

NUMERICAL COMPUTATION OF THE TRANSFORM
We have already demonstrated in part I of the paper that the discrete 2D FT in polar coordinates can be interpreted as a DFT, DHT and then IDFT. This interpretation is quite helpful in coding the transform and in exploiting the speed of the FFT (Fast Fourier Transform) in implementing the computations. In this section, we explain how the speed of Matlab's (Mathworks 2018) built-in code (or similar software) can be exploited to implement the 2D DFT in polar coordinates.

Forward transform
The values f pk can be considered as the entries in a matrix. To transform f pk ! F qm , the operation is performed as a sequence of steps which are a 1D DFT (column-wise), followed by a scaled 1D DHT (row-wise), finally followed by a 1D IDFT (column-wise).
The reader is reminded that the range of indices is given by m; k ¼ 1 . . . N 1 À 1 and n; p; q ¼ ÀM . . . M, where 2M þ 1 ¼ N 2 . These steps can be summarized succinctly by rewriting Eq. (1) as where the DHT is defined in Baddour & Chouinard (2015) via the transformation matrix Matlab code for the DHT is described in Baddour & Chouinard (2017). The inverse 2D DFT can be similarly interpreted, as shown in "Inverse Transform".

Inverse transform
The steps of the inverse 2D DFT are the reverse of the steps outlined above for the forward 2D DFT. For p ¼ ÀM . . . M and k ¼ 1 . . . N 1 À 1, Eq. (2) this can be expressed as This parallels the steps taken for the continuous case, with each continuous operation (Fourier series, Hankel transform) replaced by its discrete counterpart (DFT, DHT).
Therefore, for both forward and inverse 2D-DFT, the sequence of operations is a DFT of each column of the starting matrix, followed by a DHT of each row, a term-by-term scaling, followed by an IDFT of each column. This is a significant computational improvement because by interpreting the transform this way, the Fast Fourier Transform (FFT) can be used, which reduces the computational time quite significantly in comparison with a direct implementation of the summation definitions in Eqs. (1) and (2).

Interpretation of the sampled forward transform in Matlab terms
To use the built-in Matlab function fft, a few operations are required. First, we define Matlab-friendly indices p 0 ¼ p þ ðM þ 1Þ and n 0 ¼ n þ ðM þ 1Þ so that p; n ¼ ÀM .
. That is, the primed variables range from 1 . . . 2M þ 1 rather than ÀM . . . M. Hence, if the matrix f with entries f p 0 k is defined, where p 0 ¼ 1 . . . N 2 ; k ¼ 1 . . . N 1 À 1, then the first step in which is a column-wise DFT can be written as the Matlab-defined DFT as f pk e À2piðp 0 À1ÀMÞðn 0 À1ÀMÞ N 2 The overbar denotes a DFT. The definition of DFT in Matlab is actually given by the relationship Since the relationship P N 2 p 0 ¼1 f p 0 k e À2piðp 0 À1Þðn 0 À1ÀMÞ N 2 ¼ P N 2 p 0 ¼1 f pk e À2piðp 0 À1ÀMÞðn 0 À1ÀMÞ N 2 is valid, we can sample the original function to obtain the discrete f pk values, put them in the matrix f p 0 k then shift the matrix f p 0 k by M þ 1 along the column direction. In Matlab, the function circshift A; K; dim ð Þcan be used, which circularly shifts the values in array A by K positions along dimension dim. Inputs K and dim must be scalars. Specifically, dim ¼ 1 indicates the columns of matrix A and dim ¼ 2 indicates the rows of matrix A. Hence, Eq. (38) can be written as In matrix operations, this is equivalent to stating that each column of f p 0 k is DFT'ed to yield f n 0 k . The second step in Eq. (35) is a DHT of order n, transforming f n 0 k ! f n 0 l so that the k subscript is Hankel transformed to the l subscript. The overhat denotes a DHT. In order to relate the order n to the index n 0 , we need to shift f n 0 k by ÀðM þ 1Þ along column direction so that the order ranges from -M to M.
By using the Hankel transform matrix defined in Baddour & Chouinard (2015), Eq. (41) can be rewritten aŝ In matrix operations, this states that each row of f n 0 k is DHT'ed to yield f n 0 l . These are now scaled to give the Fourier coefficients of the 2D DFT f n 0 l ! F n 0 l . In order to proceed to an IDFT in the next step, it is necessary to shift the matrix by M þ 1 along the column direction after scaling This last step is a 1D IDFT for each column of F n 0 l to obtain F ql . Using 2M þ 1 ¼ N 2 , and q 0 ¼ q þ 1 þ M, this can be written as

Interpretation of the sampled inverse transform in Matlab terms
Similar to the forward transform, matlab-friendly indices q 0 ¼ q þ ðM þ 1Þ and n 0 ¼ n þ ðM þ 1Þ are also defined. Hence, if the matrix F with entries F q 0 l is defined, where q 0 ¼ 1 . . . N 2; l ¼ 1 . . . N 1 À 1, it then follows that the first 1D DFT step in Eq. (37) can be written as the Matlab-defined DFT as F n 0 l ¼ F ql e Àiðn 0 ÀMÀ1Þ 2pðq 0 À1ÀMÞ N 2 If the original function can be sampled as F ql and then put into matrix F q 0 l , then we need an circshift operation. So Eq. (45) can be written as F n 0 l ¼ fft circshiftðF q 0 l ; M þ 1; 1Þ; N 2 ; 1 À Á Subsequently, a DHT of order n is required, transforming F n 0 l ! F n 0 l so that the l subscript is Hankel transformed to the k subscript. To achieve this, circshift is also needed here.
This is followed by a scaling operation to obtain F n 0 k ! f n 0 k and then a circshift by ðM þ 1Þ so that This last step is a 1D IDFT for each column of f n 0 k to get f p 0 k . Using 2M þ 1 ¼ N 2 , and p 0 ¼ p þ 1, Eq. (37) can be written as In conclusion, in this section, by using the interpretation of the kernel as sequential DFT, DHT and IDFT operations, Matlab (or similar software) built-in code can be used to efficiently implement the 2D DFT algorithm in polar coordinates.

NUMERICAL EVALUATION OF THE 2D DFT IN POLAR COORDINATES TO APPROXIMATE THE CONTINUOUS FT
In this section, the 2D DFT is evaluated for its ability to estimate the continuous FT at the selected special sampling points in the spatial and frequency domains.

Accuracy
In order to test accuracy of the 2D-DFT and 2D-IDFT to calculate approximate the continuous counterpart, the dynamic error is proposed as a metric. The dynamic error is defined as Guizar-Sicairos & Gutiérrez-Vega (2004) EðvÞ ¼ 20 log 10 CðvÞ À DðvÞ j j max DðvÞ j j where CðvÞ is the continuous forward or inverse 2D-FT and DðvÞ is the value obtained from the discrete counterpart. The dynamic error is defined as the ratio of the absolute error to the maximum amplitude of the discrete function, calculated on a log scale. Therefore, a large negative value represents an accurate discrete transform. The dynamic error is used instead of the percentage error in order to avoid division by zero.

Precision
The precision of the algorithm is an important evaluation criterion, which is tested by sequentially performing a pair of forward and inverse transforms and comparing the result to the original function. High precision indicates that numerical evaluation of the transform does not add much error. An average of the absolute error between the original function and the calculated counterpart at each sample point is used to measure the precision. It is given by where f is the original function and f Ã is the value obtained after sequentially performing a forward and then inverse transform. An ideal precision would result in the absolute error being zero.

Test functions
In this section, three test functions are chosen to evaluate the ability of the discrete transform to approximate the continuous counterpart. The first test case is the circularly symmetric Gaussian function. Given that it is circularly symmetric and that the Gaussian is continuous and smooth, the proposed DFT is expected to perform well. The second test case is "Four-term sinusoid and Sinc" function, which is not symmetric in the angular direction and suffers a discontinuity in the radial direction. The third test function presents a more challenging test function, the "Four-term sinusoid and Modified exponential" function. In this case, the test function is not circularly symmetric and it explodes at the origin (approaches infinity at the origin). Given that as shown above, the sampling grid cannot cover the area around the origin very well, a function that explodes at the origin should give more error and should provide a reasonable test case for evaluating the performance of the discrete transform. The test functions are chosen to test specific aspects of the performance of the discrete transform but also because a closed-form expression for both the function and its transform are available. This then allows a numerical evaluation of the error between the quantities computed with the 2D DFT and the quantities obtained by evaluating (sampling) the continuous (forward or inverse) transform at the grid points.

Gaussian
The first function chosen for evaluation is a circular symmetric function which is Gaussian in the radial direction. Specifically, the function in the space domain is given by where a is some real constant. Since the function is circularly symmetric, the 2D-DFT is a zeroth-order Hankel Transform (Poularikas, 2010) and is given by Àr 2 4a 2 (53) The graphs for the original function and its continuous 2D-DFT (which is also a Gaussian) are plotted with a ¼ 1 and shown in Fig. 5. From Fig. 5, the function is circular symmetric and fairly smooth in the radial direction. Moreover, the function can be considered as either an effectively space-limited function or an effectively band-limited function. For the purposes of testing it, it shall be considered as a space-limited function and Eqs. (14) and (15) will be used to proceed with the forward and inverse transform in sequence.
To perform the transform, the following variables need to be chosen: N 2 , R and N 1 . In the angular direction, since the function in the spatial domain is circularly symmetric, N 2 can be chosen to be small. Thus, N 2 ¼ 15 is chosen.
In the radial direction, from plotting the function, it can be seen that the effective space limit can be taken to be R ¼ 5 and the effective band limit can be taken to be W r ¼ 10. From Eq. (21), j 0N 1 ! R Á W r ¼ 50. Therefore, N 1 ¼ 17 is chosen (we could also have obtained a rough estimate of N 1 from Eq. (24)). However, most of the energy of the function in both the space and frequency domains is located in the center near the origin. Based on the discussion in "Conclusion", relatively large values of R and W r are needed. The effective space limit R ¼ 40 and effective band-limit W p ¼ 30 are thus chosen, which gives j 0N 1 ! R Á W r ¼ 1200. Therefore N 1 ¼ 383 is chosen in order to satisfy this constraint. Both cases discussed here (N 1 ¼ 17 and N 1 ¼ 383) are tested in following.

Forward transform
Test results with R ¼ 5, N 1 ¼ 17 are shown in Figs. 6 and 7. Figure 6 shows the sampled continuous forward transform and the discrete forward transform. Figure 7 shows the error between the sampled values of the continuous transform and the discretely calculated values. From Fig. 7, it can be observed that the error gets bigger at the center, which is as expected because the sampling grid shows that the sampling points can never attain the origin. The maximum value of the error is E max ¼ À0:9115 dB and this occurs at the center. The average error is E avg: ¼ À30:4446 dB.
Error test results with R ¼ 40, N 1 ¼ 383 are shown in Fig. 8. Similar to the previous case, the error gets larger at the center, as expected. However, the maximum value of the error is E max ¼ À8:3842 dB and this occurs at the center. The average value of the error is E avg: ¼ À63:8031 dB. Clearly, the test with R ¼ 40, N 1 ¼ 383 gives a better approximation, which verifies the discussion in "Conclusion".
With R ¼ 40, Table 3 shows the errors (max and average error) with respect to different value of N 1 and N 2 . The trends as functions of N 1 and N 2 are shown as plots in Figs. 9 and 10. From Fig. 9, it can be seen that when N 1 individually (N 2 is fixed at N 2 ¼ 15) is less than the minimum of 383 obtained from the sampling theorem, increasing N 1 will lead to smaller errors, as expected. When N 1 is bigger than the sampling-theorem threshold of 383, increasing N 1 still decreases the error which verifies the discussion about sample grid coverage in "Conclusion". Increasing N 1 tends to increase the sample grid coverage and capture more information at the center area and thus leads to smaller errors. From Fig. 10, increasing N 2 alone (i.e., without a corresponding increase in N 1 ) leads to larger errors, both Error max and Error average . Although at first counterintuitive, this result is actually reasonable because the function is radially symmetric which implies that N 2 ¼ 1 should be sufficient based on the sampling theorem for the angular direction. Therefore, increasing N 2 will not lead to a better approximation. Moreover, from the discussion of the sample grid coverage in "Conclusion", the sampling grid coverage in both domains gets worse when N 2 gets bigger because more information from the center is lost. This problem can be solved by increasing N 1 at the same time, but it could be computationally time consuming. Therefore, choosing N 2 properly is very important from the standpoint of accuracy and computational efficiency.

Inverse transform
Test results for the inverse transform with R ¼ 5, N 1 ¼ 17 are shown in Figs. 11 and 12. Figure 11 shows the sampled continuous inverse transform and discrete inverse transform and Fig. 12 shows the error between the sampled continuous and discretely calculated values. Similar to the case for the forward transform, the error gets larger at the center, which is as expected because the sampling grid shows that the sampling points never attain the center. The maximum value of the error is E max ¼ 3:1954 dB and this occurs at the center. The average of the error is E avg: ¼ À25:7799 dB.
Error test results for the inverse transform with R ¼ 40, N 1 ¼ 383 are shown in Fig. 13. In this case, the maximum value of the error is E max ¼ À12:2602 dB and this occurs at the center. The average of the error is E avg: ¼ À98:0316 dB. Table 4 shows the errors with respect to different value of N 1 and N 2 , from which Figs. 14 and 15 demonstrate the trend. From Fig. 15 it can be observed that increasing N 1 tends to improve the result but not significantly. This could be explained by the discussion for R ¼ 40, N 1 ¼ 383 that with fixed R and W r , increasing N 1 will not allow the sampling grid in the frequency domain to get any closer to the origin to capture more information. From Fig. 15, increasing N 2 (with fixed N 1 ¼ 383) leads to a worse approximation which verifies the discussion for R ¼ 40, N 1 ¼ 383.  Performing sequential 2D-DFT and 2D-IDFT results in e ¼ 4:1656 Â e À17 where e is calculated with Eq. (51). Therefore, performing sequential forward and inverse transforms does not add much error.

Four-term sinusoid & Sinc function
The second function chosen for evaluation is given by f ðr; uÞ ¼ sinðarÞ ar ½3 sinðuÞ þ sinð3uÞ þ 4 cosð10uÞ þ 12 sinð15uÞ (54) which is a sinc function in the radial direction and a four-term sinusoid in the angular direction. The graphs for the original function and the magnitude of its continuous 2D-FT with a ¼ 5 are shown in Fig. 16. From Fig. 16, the function can be considered as a bandlimited function. Therefore Eqs. (16) and (17) In the angular direction, the highest frequency term in the function in the space domain is 12sinð15uÞ. From the sampling theorem, the sampling frequency should be at least twice that of the highest frequency present in the signal. Thus, N 2 ¼ 41 is chosen in order to go a little past the minimum requirement of 31. In the radial direction, from the graphs of the original function and its 2D-FT, it can be assumed that f ðr; uÞ is space-limited at R ¼ 15 and band-limited at W r ¼ 30. However, since most of the energy in the space domain is located at the origin, a relatively large band limit should be chosen based on the discussion in "Conclusion". Therefore, W r ¼ 90, N 1 ¼ 430 are chosen.

Forward transform
The error results for the forward 2D-DFT of Four-term sinusoid & Sinc function with W r ¼ 90, N 1 ¼ 430 are shown in Fig. 17. The discrete transform does not approximate the continuous transform very well. This is expected because the function in the frequency domain is discontinuous and the sampling points close to the discontinuity will result in a very large error. The maximum value of the error is Error max ¼ 10:6535 dB and this occurs where the discontinuities are located. The average of the error is Error average ¼ À38:7831 dB. With W r ¼ 90, N 1 ¼ 430, Table 5 shows the errors with respect to different value of N 1 and N 2 , from which Figs. 18 and 19 show the trend. From Fig. 18, increasing N 1 alone tends improve the average error. The maximum error does not change with N 1 , which is reasonable because of the discontinuity of the function in the frequency domain.
From Fig. 19, increasing N 2 leads to Error max and Error average first improving and then worsening. This is reasonable because when N 2 is less than the minimum requirement of 31 from sampling theorem, the test result is actually affected by both sampling point density (from the sampling theorem) and grid coverage (discussed in "Conclusion").  Increasing N 2 should give better results from the point of view of the sampling theorem but worse grid coverage. The result from the combined effects is dependent on the function properties. In the specific case of this function, when N 2 is bigger than 31, thereby implying that the angular sampling theorem has been satisfied-the results get worse with increasing N 2 .

Inverse transform
The error results for the 2D-IDFT of Four-term sinusoid & Sinc function with W r ¼ 90, N 1 ¼ 430 are shown in Fig. 20. The maximum value of the error is Error max ¼ À8:6734 dB. The average of the error is Error average ¼ À37:8119 dB. With W r ¼ 90, N 1 ¼ 430, Table 6 shows the errors with respect to different value of N 1 and N 2 , from which Figs. 21 and 22 show the trend.  From Fig. 21, it can be observed that the increasing N 1 alone improves the average error, as was expected. However, N 1 ¼ 380 gives an apparently worse average error than the other points. This could be caused by the discontinuity of the function in the frequency From graph of the original function and its 2D-DFT, it can be assumed that f ðr; uÞ is effectively space-limited with R ¼ 20, and Fðr; cÞ is effectively band-limited with W r ¼ 15, which gives j 0N 1 % 300. This results in N 1 ¼ 96: However, since the function explodes at the center area in both domains, relatively large values of R and W r should give a better approximation. Therefore, another case with R ¼ 40, W r ¼ 30 is tested. In this case, N 1 ¼ 383 is chosen.
In the angular direction, the highest frequency term is 12 sinð15uÞ. From the sampling theorem, the sampling frequency should be at least twice of the highest frequency of signal. Thus, N 2 ¼ 41 is chosen.

Forward transform
Here, the function is tested as a space limited function and Eqs. (14) and (15) are used to proceed with the forward and inverse transform in sequence. The error results with R ¼ 40; W r ¼ 30; N 1 ¼ 383 are shown in Fig. 24. The maximum value of the error is Error max ¼ À10:1535 dB and this occurs at the center area. The average of the error is Error average ¼ À32:7619 dB. This demonstrates that the discrete function approximates the sampled values of the continuous function quite well.

Inverse transform
The error results with R ¼ 40; W r ¼ 30; N 1 ¼ 383 are shown in Fig. 25.
The maximum value of the error is Error max ¼ 0:5579 dB and this occurs at the center. The average of the error isError average ¼ À68:7317 dB.
Performing 2D-DFT and 2D-IDFT results in e ¼ 1:421 Â e À12 , where e is calculated by Eq. (51). It can be observed that even for functions with the worst properties, the proposed transform can still be used to approximate the continuous FT with fairly small errors, as long as the function is sampled properly.

SUMMARY AND CONCLUSION Accuracy and precision of the transform
The proposed discrete 2D-FT in polar coordinates demonstrates an acceptable accuracy in providing discrete estimates to the continuous FT in polar coordinates. In Baddour & Chouinard (2015), Guizar-Sicairos & Gutiérrez-Vega (2004) and Higgins & Munson (1987), the one dimensional Hankel transform of a sinc function showed similar dynamic error, which could be used as a comparative measure. Since the DHT is one step of the proposed discrete 2D-FT, and the definition of the Hankel transform is based on Abbas, Sun & Foroosh (2017), a similar dynamic error should be expected.
The test cases showed that the transform introduced very small errors (e ¼ 1:4004 Â e À12 for worst case) by performing a forward transform and an inverse