Discrete two dimensional Fourier transform in polar coordinates part II: numerical computation and approximation of the continuous transform Discrete two dimensional Fourier transform in polar coordinates part II: numerical computation and approximation of the continuous transform Xueyang Yao1 and Natalie Baddour2 1 Department of Systems Design Engineering, University of Waterloo, Waterloo, ON, Canada 2 Department of Mechanical Engineering, University of Ottawa, Ottawa, ON, Canada ABSTRACT 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. Subjects Algorithms and Analysis of Algorithms, Scientific Computing and Simulation, Theory and Formal Methods Keywords Fourier theory, DFT in polar coordinates, Polar coordinates, Multidimensional DFT, Discrete Hankel Transform, Discrete Fourier Transform, Orthogonality 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 How to cite this article Yao X, Baddour N. 2020. Discrete two dimensional Fourier transform in polar coordinates part II: numerical computation and approximation of the continuous transform. PeerJ Comput. Sci. 6:e257 DOI 10.7717/peerj-cs.257 Submitted 12 July 2019 Accepted 19 January 2020 Published 2 March 2020 Corresponding author Natalie Baddour, nbaddour@uottawa.ca Academic editor Yilun Shang Additional Information and Declarations can be found on page 36 DOI 10.7717/peerj-cs.257 Copyright 2020 Yao and Baddour Distributed under Creative Commons CC-BY 4.0 http://dx.doi.org/10.7717/peerj-cs.257 mailto:nbaddour@�uottawa.�ca https://peerj.com/academic-boards/editors/ https://peerj.com/academic-boards/editors/ http://dx.doi.org/10.7717/peerj-cs.257 http://www.creativecommons.org/licenses/by/4.0/ http://www.creativecommons.org/licenses/by/4.0/ https://peerj.com/computer-science/ 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 (Baddour, 2019). 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 (Baddour, 2019), 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, Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 2/38 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ 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) fpk to the matrix (double-subscripted series) Fql such that fpk ! Fqm is given by Fqm ¼ F fpk � � ¼ XN1�1 k¼1 XM p¼�M fpkE � qm;pk (1) where p; k; q; m; n, N1, and N2 are integers such that �M � n � M, where 2M þ 1 ¼ N2 1 � m; k; � N1 � 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, N1, and N2 are within these stated ranges. Similarly, for the inverse transform we propose fpk ¼ F�1 Fqm � � ¼ XN1�1 m¼1 XM q¼�M FqmE þ qm;pk (2) Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 3/38 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ In Eqs. (1) and (2), E�qm;pk are the kernels of the transformation. These can be chosen as the “non-symmetric” form given by E�qm;pk ¼ 1 N2 XM n¼�M Jn jnkjnm jnN1 � � j2nN1J 2 nþ1 jnkð Þ 2i�ne �i 2pnpN2 e þi 2pnqN2 Eþqm;pk ¼ 1 N2 XM n¼�M Jn jnmjnk jnN1 � � J2nþ1 jnmð Þ 2iþne þi 2pnpN2 e �i 2pnqN2 (3) Here, Jn zð Þ is the nth order Bessel function of the first kind and jnk denotes the kth zero of the nth Bessel function. The subscript (+ or −) indicated the sign on the i� and on the exponent containing the p variable; the q variable exponent then takes the opposite sign. From a matrix point of view, both fpk and Fql are N2 � N1 � 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 (Baddour, 2019). 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 F jqm R ; 2pq N2 � � �2pR2 XN1�1 k¼1 XM p¼�M f jpkR jpN1 ; 2pp N2 � � 1 N2 XM n¼�M 2i�nJn jnkjnm jnN1 � � j2nN1J 2 nþ1 jnkð Þ e �i 2pnpN2 e þi 2pnqN2 (4) 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 f jpkR jpN1 ; 2pp N2 � � � 1 2pR2 XN1�1 m¼1 XM q¼�M F jqm R ; 2pq N2 � � 1 N2 XM n¼�M 2inJn jnmjnk jnN1 � � J2nþ1 jnmð Þ e þi 2pnpN2 e �i 2pnqN2 (5) Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 4/38 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ 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 fpk via fpk ¼ f jpkR jpN1 ; 2pp N2 � � (6) Then, Fqm is calculated from the 2D DFT scaled by 2pR 2, Eq. (1), that is Fqm ¼ 2pR2F fpk � � ¼ 2pR2 XN1�1 k¼1 XM p¼�M fpkE � qm;pk (7) The factor of 2pR2 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 Fqm calculated from the 2D DFT in Eq. (7) approximates F jqm R ; 2pq N2 � � —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 Fqm via Fqm ¼ F jqm R ; 2pq N2 � � (8) Now, fpk is calculated from a scaled version of the inverse 2D DFT, Eq. (2) that is fpk ¼ 1 2pR2 F�1 Fqm � � ¼ 1 2pR2 XN1�1 m¼1 XM q¼�M FqmE þ qm;pk (9) To evaluate if the proposed transform can approximate the continuous transform, the question becomes how well fpk calculated from Eq. (9) approximates f jpkR jpN1 ; 2pp N2 � � —the values of the continuous function evaluated on the sampling grid. 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; Wr � . 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 � Wr ¼ 2pW. The variable Wr 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 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 5/38 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ 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 F jqmWr jqN1 ; 2pq N2 � � � 2p W2r XN1�1 k¼1 XM p¼�M f jpk Wr ; 2pp N2 � � 1 N2 XM n¼�M 2i�nJn jnmjnk jnN1 � � J2nþ1 jnkð Þ e �i 2pnpN2 e þi 2pnqN2 (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 f jpk Wr ; 2pp N2 � � � W2r 2p XN1�1 m¼1 XM q¼�M F jqmWr jqN1 ; 2pq N2 � � 1 N2 XM n¼�M 2inJn jnkjnm jnN1 � � j2nN1J 2 nþ1 jnmð Þ e �i 2pnqN2 e þi 2pnpN2 (11) 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 fpk via fpk ¼ f jpk Wr ; 2pp N2 � � (12) Then, Fqm is calculated from the discrete transform scaled by 2p W2r , Eq. (1), that is Fqm ¼ 2p W2r F fpk � � ¼ 2p W2r XN1�1 k¼1 XM p¼�M fpkE � qm;pk (13) To evaluate if the proposed 2D DFT can be used to approximate the continuous transform, the question is how well Fqm calculated from Eq. (13) approximates F jqmWr jqN1 ; 2pq N2 � � , which are the values of the continuous 2D FT, evaluated on the sampling grid. The evaluation of the inverse transform for the band-limited function proceeds similarly by comparing values obtained from the inverse 2D DFT to the values obtained by sampling the continuous function directly. The relationships given by Eqs. (4), (5), (10) and (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 fpk and Fqm 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 fpk to yield its forward transform Fqm, 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 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 6/38 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ 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 regular~r space r; uð Þ and ~v frequency r; cð Þ space are defined as rpk ¼ jpkR jpN1 up ¼ p2p N2 (14) and rqm ¼ jqm R cq ¼ q2p N2 (15) For a band limited function, the function is assumed band-limited to 0 � r � Wr, 0 � c � 2p. The sampling space used for radial and angular sampling points in regular ~v frequency space r; cð Þ and~r space r; uð Þ for a bandlimited function is defined as rpk ¼ jpk Wr up ¼ p2p N2 (16) and rqm ¼ jqmWr jqN1 cq ¼ q2p N2 (17) Clearly, the density of the sampling points depends on the numbers of points chosen, that is on N1 and N2. 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, N1 ¼ 16 and N2 ¼ 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 N2 and N1 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 N2 and N1, the grid becomes fairly dense, obtaining good coverage of both spaces, but details are harder to observe. To demonstrate, Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 7/38 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ the polar grids are plotted for the case R = 1, N1 ¼ 96 and N2 ¼ 95. These are shown in Figs. 3 and 4 respectively. From Figs. 3 and 4, by choosing higher values of N1 and N2, 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 fs > 2fmax (18) where fs is the sample frequency and fmax is the highest frequency or band limit. In the radial direction, the necessary relationship for the DHT is given by Baddour & Chouinard (2015) WrR ¼ jnN1 (19) Figure 1 Spatial sampling grid for a space-limited function with R = 1, N1 = 16 and N2 = 15. Full-size DOI: 10.7717/peerj-cs.257/fig-1 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 8/38 http://dx.doi.org/10.7717/peerj-cs.257/fig-1 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ where Wr is the effective band-limit, R is the effective space limit and jnN is the Nth zero of Jn rð Þ. For the 2D FT, since �M � p � M, the order of the Bessel zero ranges from �M to M, the required relationship becomes minðjpN1Þ � WrR (20) The relationships jnN ¼ j�nN and j0N1 < j�1N1 < j�2N1 < … < j�MN1 are valid (Lozier, 2003), hence Eq. (20) can be written as j0N1 � WrR (21) It is pointed out in Baddour (2019) and Guizar-Sicairos & Gutiérrez-Vega (2004) that the zeros of Jn 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 Jn zð Þ � ffiffiffiffiffiffi 2 pz r cos z � n þ 1 2 � � p 2 � � (22) Figure 2 Frequency space sampling grid for a space-limited function with R = 1, N1 = 16 and N2 = 15. Full-size DOI: 10.7717/peerj-cs.257/fig-2 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 9/38 http://dx.doi.org/10.7717/peerj-cs.257/fig-2 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ Therefore, an approximation to the Bessel zero, jnk is given by jnk � 2k þ n � 1 2 � � p 2 (23) Hence, Eq. (21) can be written to choose N1 approximately as N1p � WrR ¼ 2pWR ) N1 � 2WR (24) where the reader is reminded that the units of W is m−1 (the space equivalent of Hz). N1=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 N1 or N2 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 N2 and N1. 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 N2, N1 and the size of the gap will be discussed. Figure 3 Spatial sampling grid for a space-limited function with R = 1, N1 = 96 and N2 = 95. Full-size DOI: 10.7717/peerj-cs.257/fig-3 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 10/38 http://dx.doi.org/10.7717/peerj-cs.257/fig-3 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ 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 N2, N1 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 N2 and N1 affect the coverage of the grid in space domain, consider the following definition of ‘grid coverage’ Ar ¼ pR2 � p�r2 pR2 100 (25) where �r denotes the average radius of the gap (the hole in the middle of the grid). Ar 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 Ar would be 100%, that is, complete coverage. Based on the observation of Figs. 1 and 3, the relationship Figure 4 Frequency space sampling grid for a space-limited function with R = 1, N1 = 96 and N2 = 95. Full-size DOI: 10.7717/peerj-cs.257/fig-4 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 11/38 http://dx.doi.org/10.7717/peerj-cs.257/fig-4 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ r01 < r�11 < r�21 < r�M1 is valid. Therefore, from Eq. (14), the average radius of the gap is given by �r ¼ ðr01 þ rM1Þ 2 ¼ 1 2 j01 j0N1 R þ jM1 jMN1 R � � (26) Hence, Eq. (25) for grid coverage can be written as Ar ¼ 1 � 1 4 j01 j0N1 þ jM1 jMN1 � �2" # 100 (27) Table 1 shows the different values of grid coverage Ar as the values of N1 and N2 are changed. From Table 1, it can be seen that increasing N1 (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 N1 actually increases the effective band limit. However, increasing N2 (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 Ar ¼ pW2r � p�r2 pW2r 100 (28) where �r denotes the average radius of the gap. Since �r ¼ ðr01 þ rM1Þ 2 ¼ ðj01 þ jM1Þ 2R (29) Then, it follows that Eq. (28) for frequency domain grid coverage can be written as Ar ¼ 1 � ðj01 þ jM1Þ2 4R2W2r " # 100% (30) From Eq. (30), it can be observed that the sample grid coverage in the frequency domain is affected by R, Wr and M. Since N2 ¼ 2M þ 1, in order to get a better grid coverage with a fixed Wr, R and N2 can be adjusted. Table 2 shows the grid coverage Ar for different values of R and N2. 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 Table 1 Spatial grid coverage, Ar, with respect to different values of N1 and N2 (R is fixed). N2 N1 15 75 150 300 15 Ar = 98.48% Ar = 99.92% Ar = 99.98% Ar = 99.99% 75 Ar = 93.78% Ar = 99.36% Ar = 99.81% Ar = 99.95% 151 Ar = 90.14% Ar = 98.42% Ar = 99.46% Ar = 99.84% 301 Ar = 86.17% Ar = 96.58% Ar = 98.59% Ar = 99.51% Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 12/38 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ frequency domain, while increasing N2 (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 Wr requires an increase in N1 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; Wp�. 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 N2, N1 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 �r ¼ ðr01 þ rM1Þ 2 ¼ 1 2 j01 Wr þ jM1 Wr � � (31) Therefore, the coverage of the grid in space domain can be written as Ar¼ 1 � ðj01 þ jM1Þ2 4W2rR 2 " # 100 (32) 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 �r ¼ ðr01 þ rM1Þ 2 ¼ 1 2 j01 j0N1 Wr þ jM1 jMN1 Wr � � (33) The coverage of the grid in frequency domain can be written as Ar ¼ 1 � 1 4 j01 j0N1 þ jM1 jMN1 � �2" # 100 (34) 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. Table 2 Frequency grid coverage, Aρ, with respect to different values of R and N2 (Wρ is fixed). N2 R 15 75 150 300 15 Aρ = 99.80% Aρ = 99.99% Aρ = 100.00% Aρ = 100.00% 75 Aρ = 97.66% Aρ = 99.91% Aρ = 99.98% Aρ = 99.99% 151 Aρ = 91.88% Aρ = 99.68% Aρ = 99.92% Aρ = 99.98% 301 Aρ = 70.67% Aρ = 98.83% Aρ = 99.71% Aρ = 99.93% Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 13/38 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ Conclusion Based on the discussion above, the following conclusions can be made: 1. Increasing N2 (angular direction) tends to decrease the sampling grid coverage in both domains. Increasing N1 (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 N2 are needed, a large value of N1 is also needed to compensate for the effect of increasing N2 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 N1 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 N1 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 N1 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 N1 and Wr 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 fpk can be considered as the entries in a matrix. To transform fpk ! Fqm, 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 . . . N1 � 1 and n; p; q ¼ �M . . . M, where 2M þ 1 ¼ N2. These steps can be summarized succinctly by rewriting Eq. (1) as Fqm ¼ 1 N2 XM n¼�M 2pR2i�n jnN1 XN1�1 k¼1 YnN1m;k XM p¼�M fpke �in 2ppN2 |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} 1D DFT column‐wise 2 66664 3 77775 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} scaled 1D DHT row‐wise e þin 2pqN2 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} inverse 1D DFT column‐wise (35) Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 14/38 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ where the DHT is defined in Baddour & Chouinard (2015) via the transformation matrix YnN1m;k ¼ 2 jnN1J 2 nþ1 jnkð Þ Jn jnmjnk jnN1 � � 1 � m; k � N1 � 1 (36) 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 . . . N1 � 1, Eq. (2) this can be expressed as fpk ¼ 1 N2 XM n¼�M jnN1i þn 2pR2 XN1�1 m¼1 YnN1k;m XM q¼M Fqme �i 2pnqN2 |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} 1D DFT ðcolumn‐wiseÞ 2 66664 3 77775 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} scaled 1D DFT ðrow‐wiseÞ e þi 2pnpN2 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} inverse 1D DFT ðcolumn‐wiseÞ (37) 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 p0 ¼ p þ ðM þ 1Þ and n0 ¼ n þ ðM þ 1Þ so that p; n ¼ �M . . . M become p0; n0 ¼ 1 . . . 2M þ 1 ¼ 1 . . . N2 (since 2M þ 1 ¼ N2). That is, the primed variables range from 1 . . . 2M þ 1 rather than �M . . . M. Hence, if the matrix f with entries fp0k is defined, where p 0 ¼ 1 . . . N2; k ¼ 1 . . . N1 � 1, then the first step in which is a column-wise DFT can be written as the Matlab-defined DFT as �fn0k ¼ XN2 p0¼1 fpke �2piðp0�1�MÞðn0�1�MÞ N2 (38) The overbar denotes a DFT. The definition of DFT in Matlab is actually given by the relationship �fn0k ¼ XN2 p0¼1 fp0ke �2piðp0�1Þðn0�1Þ N2 (39) Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 15/38 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ Since the relationship PN2 p0¼1 fp0ke �2piðp0�1Þðn0�1�MÞ N2 ¼ PN2p0¼1 fpke�2piðp 0�1�MÞðn0�1�MÞ N2 is valid, we can sample the original function to obtain the discrete fpk values, put them in the matrix fp0k then shift the matrix fp0k 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 �fn0k ¼ fft circshift fp0k; M þ 1; 1 � � ; N2; 1 � � (40) In matrix operations, this is equivalent to stating that each column of fp0k is DFT’ed to yield �fn0k. The second step in Eq. (35) is a DHT of order n, transforming �fn0k ! �̂f n0l 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 n0, we need to shift �fn0k by �ðM þ 1Þ along column direction so that the order ranges from –M to M. �̂f n0l ¼ XN1�1 k¼1 2Jn jnkjnl jnN1 � � jnN1J 2 nþ1 jnkð Þ circshift �f n0k;�ðM þ1Þ;1 � � for n0 ¼ 1...N2; l ¼ 1...N1 �1 where n ¼ n0 � M þ1ð Þ � (41) By using the Hankel transform matrix defined in Baddour & Chouinard (2015), Eq. (41) can be rewritten as �̂f n0l ¼ circshift �f n0k;�ðM þ 1Þ; 1 � � YnN1l;k � �T for n0 ¼ 1 . . . N2; l ¼ 1 . . .N1 � 1 where n ¼ n0 � M � 1 � (42) In matrix operations, this states that each row of �fn0k is DHT’ed to yield �̂f n0l. These are now scaled to give the Fourier coefficients of the 2D DFT �̂f n0l ! �Fn0l. 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 �Fn0l ¼ circshift 2pR2 jnN1 i�n�̂f n0l; M þ 1; 1 � � for n0 ¼ 1 . . . N2; l ¼ 1 . . . N1 � 1 where n ¼ n0 � M þ 1ð Þ � (43) This last step is a 1D IDFT for each column of �Fn0l to obtain Fql. Using 2M þ 1 ¼ N2, and q0 ¼ q þ 1 þ M, this can be written as Fq0l ¼ 1 N2 XN2 n0¼1 �Fnle þi n0�M�1ð Þ 2p q 0�1�Mð Þ N2 for q0 ¼ 1 :: N2; l ¼ 1 :: N1 � 1 ¼ 1 N2 XN2 n0¼1 �Fn0le þiðn0�1Þ 2p q 0�1�Mð Þ N2 ¼ circshift ifft �Fn0l; N2; 1ð Þ; � M þ 1ð Þ; 1ð Þ (44) Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 16/38 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ Interpretation of the sampled inverse transform in Matlab terms Similar to the forward transform, matlab-friendly indices q0 ¼ q þ ðM þ 1Þ and n0 ¼ n þ ðM þ 1Þ are also defined. Hence, if the matrix F with entries Fq0l is defined, where q0 ¼ 1 . . . N2; l ¼ 1 . . . N1 � 1, it then follows that the first 1D DFT step in Eq. (37) can be written as the Matlab-defined DFT as �Fn0 l ¼ XN2 q0¼1 Fqle �iðn0�M�1Þ 2pðq 0�1�MÞ N2 for n0 ¼ 1 . . . N2; l ¼ 1 . . . N1 � 1 ¼ XN2 q0¼1 Fq0le �iðn0�M�1Þ 2pðq 0�1Þ N2 (45) If the original function can be sampled as Fql and then put into matrix Fq0l, then we need an circshift operation. So Eq. (45) can be written as �Fn0l ¼ fft circshiftðFq0l; M þ 1; 1Þ; N2; 1 � � (46) Subsequently, a DHT of order n is required, transforming �Fn0l ! �̂Fn0l so that the l subscript is Hankel transformed to the k subscript. To achieve this, circshift is also needed here. �̂Fn0k ¼ circshift �Fn0l; �ðM þ 1Þ; 1ð Þ YnN1k;l � �T for n0 ¼ 1 . . . N2; l ¼ 1 . . . N1 � 1 where n ¼ n0 � M � 1 � (47) This is followed by a scaling operation to obtain �̂Fn0k ! �fn0k and then a circshift by ðM þ 1Þ so that �fn0k ¼ circshift jnN1 2pR2 iþn�̂Fn0k; ðM þ 1Þ; 1 � � for n0 ¼ 1 . . . N2; k ¼ 1 . . . N1 � 1 where n ¼ n0 � M þ 1ð Þ � (48) This last step is a 1D IDFT for each column of �fn0k to get fp0k. Using 2M þ 1 ¼ N2, and p0 ¼ p þ 1, Eq. (37) can be written as fp0k ¼ 1 N2 XN2 n0¼1 �f nke þi n0�M�1ð Þ 2p p 0�1�Mð Þ N2 for p0 ¼ 1 . . . N2; k ¼ 1 . . . N1 � 1 ¼ 1 N2 XN2 n0¼1 �f n0ke þi 2p n 0�1ð Þ p0�1�Mð Þ N2 ¼ circshift ifft �f n0k; N2; 1 � � ; �ðM þ 1Þ; 1 � � (49) 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. Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 17/38 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ Method for testing the algorithm 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 log10 CðvÞ � DðvÞj j max DðvÞj j � � (50) 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 e ¼ 1 N1 � 1ð Þ N2 XN1�1ð Þ N2 n¼1 f � f j j (51) 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 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 18/38 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ 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 f ðr; uÞ ¼ e�a2r2 (52) 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 Fðr; cÞ ¼ p a2 e �r2 4a2 (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: N2, R and N1. In the angular direction, since the function in the spatial domain is circularly symmetric, N2 can be chosen to be small. Thus, N2 ¼ 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 Wr ¼ 10. From Eq. (21), j0N1 � R Wr ¼ 50. Therefore, N1 ¼ 17 is chosen (we could also have obtained a rough estimate of N1 from Eq. (24)). However, most of the energy of the Figure 5 (A) Original function (Gaussian) and (B) its continuous 2D-DFT (which is also a Gaussian). Full-size DOI: 10.7717/peerj-cs.257/fig-5 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 19/38 http://dx.doi.org/10.7717/peerj-cs.257/fig-5 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ 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 Wr are needed. The effective space limit R ¼ 40 and effective band-limit Wp ¼ 30 are thus chosen, which gives j0N1 � R Wr ¼ 1200. Therefore N1 ¼ 383 is chosen in order to satisfy this constraint. Both cases discussed here (N1 ¼ 17 and N1 ¼ 383) are tested in following. Forward transform Test results with R ¼ 5, N1 ¼ 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 Emax ¼ �0:9115 dB and this occurs at the center. The average error is Eavg: ¼ �30:4446 dB. Error test results with R ¼ 40, N1 ¼ 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 Emax ¼ �8:3842 dB and this occurs at the center. The average value of the error is Eavg: ¼ �63:8031 dB. Clearly, the test with R ¼ 40, N1 ¼ 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 N1 and N2. The trends as functions of N1 and N2 are shown as plots in Figs. 9 and 10. From Fig. 9, it can be seen that when N1 individually (N2 is fixed at N2 ¼ 15) is less than the minimum of 383 obtained from the sampling theorem, increasing N1 will lead to smaller errors, as expected. When N1 is bigger than the sampling-theorem threshold Figure 6 (A) sampled continuous transform and (B) discrete forward transform for a Gaussian function with R = 5 and N1 = 17. Full-size DOI: 10.7717/peerj-cs.257/fig-6 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 20/38 http://dx.doi.org/10.7717/peerj-cs.257/fig-6 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ of 383, increasing N1 still decreases the error which verifies the discussion about sample grid coverage in “Conclusion”. Increasing N1 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 N2 alone (i.e., without a corresponding increase in N1) leads to larger errors, both Errormax and Erroraverage. Although at first counterintuitive, this result is actually reasonable because the function is radially symmetric which implies that N2 ¼ 1 should be sufficient based on the sampling theorem for the angular direction. Figure 7 Error between the sampled values of the continuous transform and the discretely calculated values for a Gaussian function with R = 5 and N1 = 17. Full-size DOI: 10.7717/peerj-cs.257/fig-7 Figure 8 Error between the sampled values of the continuous transform and the discretely calculated values for a Gaussian function with R = 40 and N1 = 383.Full-size DOI: 10.7717/peerj-cs.257/fig-8 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 21/38 http://dx.doi.org/10.7717/peerj-cs.257/fig-7 http://dx.doi.org/10.7717/peerj-cs.257/fig-8 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ Therefore, increasing N2 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 N2 gets bigger because more information from the center is lost. This problem can be solved by increasing N1 at the same time, but it could be computationally time consuming. Therefore, choosing N2 properly is very important from the standpoint of accuracy and computational efficiency. Figure 9 Error trend between the sampled values of the continuous transform and the discretely calculated values for a Gaussian function, as a function of N1. Full-size DOI: 10.7717/peerj-cs.257/fig-9 Table 3 Error (dB) of forward transform of Gaussian function with R = 40, different value of N1 and N2. N2 N1 283 333 383 433 483 3 Emax. = −21.6 Emax. = −23.0 Emax. = −24.3 Emax. = −25.4 Emax. = −26.3 Eavg. = −71.3 Eavg. = −76.9 Eavg. = −81.8 Eavg. = −86.0 Eavg. = −89.8 7 Emax. = −12.9 Emax. = −14.4 Emax. = −15.7 Emax. = −16.9 Emax. = −17.8 Eavg. = −62.6 Eavg. = −68.3 Eavg. = −73.2 Eavg. = −77.5 Eavg. = −81.4 15 Emax. = −5.4 Emax. = −7.0 Emax. = −8.4 Emax. = −9.6 Emax. = −10.6 Eavg. = −53.1 Eavg. = −58.9 Eavg. = −63.8 Eavg. = −68.1 Eavg. = −72.0 31 Emax. = 2.3 Emax. = 0.5 Emax. = −1.0 Emax. = −2.3 Emax. = −3.4 Eavg. = −42.0 Eavg. = −47.6 Eavg. = −52.5 Eavg. = −56.9 Eavg. = −60.7 61 Emax. = 9.7 Emax. = 7.9 Emax. = 6.4 Emax. = 5.0 Emax. = 3.8 Eavg. = −32.5 Eavg. = −37.5 Eavg. = −42.0 Eavg. = −46.1 Eavg. = −49.8 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 22/38 http://dx.doi.org/10.7717/peerj-cs.257/fig-9 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ Inverse transform Test results for the inverse transform with R ¼ 5, N1 ¼ 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. Figure 11 (A) sampled continuous inverse transform and (B) discrete inverse transform for the Gaussian function for R = 5 and N1 = 17. Full-size DOI: 10.7717/peerj-cs.257/fig-11 Figure 10 Error trend between the sampled values of the continuous transform and the discretely calculated values for a Gaussian function, as a function of N2. Full-size DOI: 10.7717/peerj-cs.257/fig-10 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 23/38 http://dx.doi.org/10.7717/peerj-cs.257/fig-11 http://dx.doi.org/10.7717/peerj-cs.257/fig-10 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ 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 Emax ¼ 3:1954 dB and this occurs at the center. The average of the error is Eavg: ¼ �25:7799 dB. Error test results for the inverse transform with R ¼ 40, N1 ¼ 383 are shown in Fig. 13. In this case, the maximum value of the error is Emax ¼ �12:2602 dB and this occurs at the center. The average of the error is Eavg: ¼ �98:0316 dB. Table 4 shows the errors with respect to different value of N1 and N2, from which Figs. 14 and 15 demonstrate the trend. Figure 13 Error between the sampled continuous inverse transform and discrete inverse transform for the Gaussian function for R = 40 and N1 = 383. Full-size DOI: 10.7717/peerj-cs.257/fig-13 Figure 12 Error between the sampled continuous inverse transform and discrete inverse transform for the Gaussian function for R = 5 and N1 = 17. Full-size DOI: 10.7717/peerj-cs.257/fig-12 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 24/38 http://dx.doi.org/10.7717/peerj-cs.257/fig-13 http://dx.doi.org/10.7717/peerj-cs.257/fig-12 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ From Fig. 15 it can be observed that increasing N1 tends to improve the result but not significantly. This could be explained by the discussion for R ¼ 40, N1 ¼ 383 that with fixed R and Wr, increasing N1 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 N2 (with fixed N1 ¼ 383) leads to a worse approximation which verifies the discussion for R ¼ 40, N1 ¼ 383. Figure 14 Error trend between the sampled values of the continuous inverse transform and the discretely calculated values for a Gaussian function, as a function of N1. Full-size DOI: 10.7717/peerj-cs.257/fig-14 Table 4 Error (dB) of inverse transform of Gaussian function with R = 40, different value of N1 and N2. N2 N1 283 333 383 433 483 3 Emax. = −25.9 Emax. = −27.5 Emax. = −28.9 Emax. = −30.2 Emax. = −31.3 Eavg. = −115.3 Eavg. = −115.4 Eavg. = −115.4 Eavg. = −115.5 Eavg. = −115.5 7 Emax. = −16.5 Emax. = −18.1 Emax. = −19.4 Emax. = −20.5 Emax. = −21.6 Eavg. = −107.0 Eavg. = −107.1 Eavg. = −107.2 Eavg. = −107.2 Eavg. = −107.2 15 Emax. = −9.7 Emax. = −11.0 Emax. = −12.3 Emax. = −13.4 Emax. = −14.4 Eavg. = −97.9 Eavg. = −98.0 Eavg. = −98.0 Eavg. = −98.1 Eavg. = −98.1 34 Emax. = −4.4 Emax. = −5.5 Emax. = −6.5 Emax. = −7.5 Emax. = −8.3 Eavg. = −86.9 Eavg. = −86.9 Eavg. = −87.0 Eavg. = −87.0 Eavg. = −87.0 61 Emax. = −1.1 Emax. = −1.7 Emax. = −2.4 Emax. = −3.0 Emax. = −3.7 Eavg. = −75.6 Eavg. = −75.6 Eavg. = −75.6 Eavg. = −75.6 Eavg. = −75.7 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 25/38 http://dx.doi.org/10.7717/peerj-cs.257/fig-14 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ 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 band- limited function. Therefore Eqs. (16) and (17) were used to implement the forward and inverse transform. The continuous 2D-FT can be calculated from Baddour (2011) Fðr; cÞ ¼ X1 n¼�1 2pi�neinc Z1 0 fnðrÞJnðrrÞrdr (55) where fnðrÞ is the Fourier series of f ðr; uÞ and can be written as fnðrÞ ¼ 1 2p Zp �p f ðr; uÞe�inudu (56) Figure 15 Error trend between the sampled values of the continuous inverse transform and the discretely calculated values for a Gaussian function, as a function of N2. Full-size DOI: 10.7717/peerj-cs.257/fig-15 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 26/38 http://dx.doi.org/10.7717/peerj-cs.257/fig-15 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ From the sampling theorem for the angular direction, the highest angular frequency in Eq. (54) results in N2 being at least 31 required to reconstruct the signal. Therefore, at least 31 terms are required to calculate the continuous 2D-FT, which can be written as Fðr;cÞ¼ 8pcosð10cÞr10 a ffiffiffiffiffiffiffiffiffiffiffiffiffi a2�r2 p ðaþ ffiffiffiffiffiffiffiffiffiffiffiffiffi a2�r2 p Þ10 ; r, a � 6pisinðcÞ ar ffiffiffiffiffiffiffiffiffiffiffiffiffi r2þa2 p þ2pisin 3arcsin a r � �� � sinð3cÞffiffiffiffiffiffiffiffiffiffiffiffiffi r2þa2 p �8psin 10arcsin a r � �� � cosð10cÞffiffiffiffiffiffiffiffiffiffiffiffiffi r2þa2 p þ 24pisin 15arcsin a r � �� � sinð15cÞffiffiffiffiffiffiffiffiffiffiffiffiffi r2þa2 p ; r.a 8>>>>>>>>>>>>>< >>>>>>>>>>>>>: (57) 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, N2 ¼ 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 Wr ¼ 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, Wr ¼ 90, N1 ¼ 430 are chosen. Forward transform The error results for the forward 2D-DFT of Four-term sinusoid & Sinc function with Wr ¼ 90, N1 ¼ 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 Errormax ¼ 10:6535 dB Figure 16 Plots of the (A) original function (four-term sinusoid and sinc) and (B) the magnitude of its continuous forward 2D Fourier transform with a = 5. Full-size DOI: 10.7717/peerj-cs.257/fig-16 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 27/38 http://dx.doi.org/10.7717/peerj-cs.257/fig-16 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ and this occurs where the discontinuities are located. The average of the error is Erroraverage ¼ �38:7831 dB. With Wr ¼ 90, N1 ¼ 430, Table 5 shows the errors with respect to different value of N1 and N2, from which Figs. 18 and 19 show the trend. From Fig. 18, increasing N1 alone tends improve the average error. The maximum error does not change with N1, which is reasonable because of the discontinuity of the function in the frequency domain. From Fig. 19, increasing N2 leads to Errormax and Erroraverage first improving and then worsening. This is reasonable because when N2 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”). Figure 17 Error results for the forward 2D Fourier transform of the Four-term sinusoid & Sinc function for Wp = 90 and N1 = 430. Full-size DOI: 10.7717/peerj-cs.257/fig-17 Table 5 Error (dB) of the forward transform of ‘four-term sinusoid & Sinc’ function with different value of N1 and N2 of forward transform. N2 N1 330 380 430 480 530 11 Emax. = 4.6 Emax. = 7.1 Emax. = 3.4 Emax. = 9.0 Emax. = 2.8 Eavg. = −33.6 Eavg. = −33.4 Eavg. = −33.5 Eavg. = −35.1 Eavg. = −35.5 21 Emax. = 6.7 Emax. = 10.5 Emax. = 3.2 Emax. = 6.9 Emax. = 3.6 Eavg. = −33.9 Eavg. = −34.6 Eavg. = −37.2 Eavg. = −38.0 Eavg. = −38.1 41 Emax. = 8.5 Emax. = 35.1 Emax. = 10.7 Emax. = 14.6 Emax. = 11.1 Eavg. = −38.7 Eavg. = −38.9 Eavg. = −38.8 Eavg. = −39.8 Eavg. = −41.3 81 Emax. = 9.7 Emax. = 32.7 Emax. = 14.8 Emax. = 22.6 Emax. = 14.5 Eavg. = −34.3 Eavg. = 35.5 Eavg. = −36.2 Eavg. = −37.3 Eavg. = −37.5 161 Emax. = 19.9 Emax. = 30.2 Emax. = 22.5 Emax. = 22.5 Emax. = 16.1 Eavg. = −29.4 Eavg. = −30.7 Eavg. = −31.1 Eavg. = −32.2 Eavg. = −32.8 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 28/38 http://dx.doi.org/10.7717/peerj-cs.257/fig-17 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ Increasing N2 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 N2 is bigger than 31, thereby Figure 18 Error trend between the sampled values of the continuous forward transform and the discretely calculated values for a four-term sinusoid and sinc as a function of N1. Full-size DOI: 10.7717/peerj-cs.257/fig-18 Figure 19 Error trend between the sampled values of the continuous forward transform and the discretely calculated values for a four-term sinusoid and sinc as a function of N2. Full-size DOI: 10.7717/peerj-cs.257/fig-19 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 29/38 http://dx.doi.org/10.7717/peerj-cs.257/fig-18 http://dx.doi.org/10.7717/peerj-cs.257/fig-19 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ implying that the angular sampling theorem has been satisfied—the results get worse with increasing N2. Inverse transform The error results for the 2D-IDFT of Four-term sinusoid & Sinc function with Wr ¼ 90, N1 ¼ 430 are shown in Fig. 20. The maximum value of the error is Errormax ¼ �8:6734 dB. The average of the error is Erroraverage ¼ �37:8119 dB. With Wr ¼ 90, N1 ¼ 430, Table 6 shows the errors with respect to different value of N1 and N2, from which Figs. 21 and 22 show the trend. Figure 20 Error results for the 2D inverse discrete Fourier transform of the four-term sinusoid and sinc function for Wp = 90 and N1 = 430. Full-size DOI: 10.7717/peerj-cs.257/fig-20 Table 6 Error (dB) of inverse transform of ‘four-term sinusoid & Sinc’ function with different value of N1 and N2. N2 N1 330 380 430 480 530 11 Emax. = 0.1 Emax. = 0.1 Emax. = 0.1 Emax. = 0.1 Emax. = 0.1 Eavg. = −43.7 Eavg. = −43.7 Eavg. = −46.6 Eavg. = −45.6 Eavg. = −48.1 21 Emax. = 0.7 Emax. = 0.7 Emax. = 0.6 Emax. = 0.6 Emax. = 0.7 Eavg. = −38.3 Eavg. = −38.0 Eavg. = −40.4 Eavg. = −40.6 Eavg. = −42.2 41 Emax. = −9.0 Emax. = −8.5 Emax. = −8.7 Emax. = −8.8 Emax. = −8.6 Eavg. = −35.9 Eavg. = −24.7 Eavg. = −37.8 Eavg. = −38.2 Eavg. = −39.0 81 Emax. = −4.5 Emax. = −4.7 Emax. = −4.5 Emax. = −4.6 Emax. = −4.5 Eavg. = −35.7 Eavg. = −26.5 Eavg. = −37.5 Eavg. = −36.2 Eavg. = −39.0 161 Emax. = 0.8 Emax. = 0.7 Emax. = 0.7 Emax. = 0.7 Emax. = 0.7 Eavg. = −35.6 Eavg. = −32.5 Eavg. = −36.6 Eavg. = −37.2 Eavg. = −39.2 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 30/38 http://dx.doi.org/10.7717/peerj-cs.257/fig-20 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ From Fig. 21, it can be observed that the increasing N1 alone improves the average error, as was expected. However, N1 ¼ 380 gives an apparently worse average error than the other points. This could be caused by the discontinuity of the function in the frequency Figure 21 Error trend between the sampled values of the continuous inverse transform and the discretely calculated values for a four-term sinusoid and sinc function, as a function of N1. Full-size DOI: 10.7717/peerj-cs.257/fig-21 Figure 22 Error trend between the sampled values of the continuous inverse transform and the discretely calculated values for a four-term sinusoid and sinc function, as a function of N2. Full-size DOI: 10.7717/peerj-cs.257/fig-22 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 31/38 http://dx.doi.org/10.7717/peerj-cs.257/fig-21 http://dx.doi.org/10.7717/peerj-cs.257/fig-22 http://dx.doi.org/10.7717/peerj-cs.257 https://peerj.com/computer-science/ domain. Changing to N1 ¼ 381, the average error becomes −37.0 dB which proves that the large error is caused by the discontinuity. From Fig. 22, increasing N2 does not lead to worse results, which is different from previous cases. However, from Fig. 16 it can be seen that the function in the frequency domain does not have much information in the center area. So, even though increasing N2 causes a larger hole in the center as discussed in “Conclusion”, it does not lead to losing much energy which explains why Fig. 22 shows a different trend from the previous cases. Performing 2D-DFT and 2D-IDFT sequentially results in e ¼ 1:3117 � e�12 where e is calculated by Eq. (51). Four-term sinusoid and modified exponential For the next test function, the function is given by f ðr; uÞ ¼ e �ar r ½3 sinðuÞ þ sinð3uÞ þ 4 cosð10uÞ þ 12 sinð15uÞ� (58) Its continuous 2D-FT can be calculated as Fðr; cÞ ¼ �6pi sinðcÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r2 þ a2 p � a r ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r2 þ a2 p þ 2pi sinð3cÞ ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r2 þ a2 p � aÞ3 r3 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r2 þ a2 p � 8p cosð10cÞ ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r2 þ a2 p � aÞ10 r10 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r2 þ a2 p þ 24pi sinð15cÞ ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r2 þ a2 p � aÞ15 r15 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r2 þ a2 p (59) The graphs for the original function and the magnitude of its continuous 2D-FT with a = 0.1 are shown in Fig. 23. From Fig. 23, it can be observed that the function has a spike in both domains, which is a more difficult scenario based on the discussion in “Conclusion”. In this case, the function can be assumed as space-limited or band-limited. This function will be tested as being space-limited. Figure 23 Plots for (A) the original function and (B) the magnitude of its continuous 2D discrete Fourier transform with a = 0.1 for a four-term sinusoid and modified exponential function. Full-size DOI: 10.7717/peerj-cs.257/fig-23 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 32/38 http://dx.doi.org/10.7717/peerj-cs.257/fig-23 https://peerj.com/computer-science/ http://dx.doi.org/10.7717/peerj-cs.257 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 Wr ¼ 15, which gives j0N1 � 300. This results in N1 ¼ 96: However, since the function explodes at the center area in both domains, relatively large values of R and Wr should give a better approximation. Therefore, another case with R ¼ 40, Wr ¼ 30 is tested. In this case, N1 ¼ 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, N2 ¼ 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; Wr ¼ 30; N1 ¼ 383 are shown in Fig. 24. The maximum value of the error is Errormax ¼ �10:1535 dB and this occurs at the center area. The average of the error is Erroraverage ¼ �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; Wr ¼ 30; N1 ¼ 383 are shown in Fig. 25. The maximum value of the error is Errormax ¼ 0:5579 dB and this occurs at the center. The average of the error isErroraverage ¼ �68:7317 dB. Performing 2D-DFT and 2D-IDFT results in e ¼ 1:421 � e�12, where e is calculated by Eq. (51). Figure 24 Error between the sampled values of the continuous forward transform and the discretely calculated values for the four-term sinusoid and modified exponential function with R = 40, Wp = 30 and N1 = 383. Full-size DOI: 10.7717/peerj-cs.257/fig-24 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 33/38 http://dx.doi.org/10.7717/peerj-cs.257/fig-24 https://peerj.com/computer-science/ http://dx.doi.org/10.7717/peerj-cs.257 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 Figure 25 Error between the sampled values of the continuous inverse transform and the discretely calculated values for the four-term sinusoid and modified exponential function with R = 40, Wp = 30 and N1 = 383. Full-size DOI: 10.7717/peerj-cs.257/fig-25 Table 7 Computing time of three cases: Case1: Run the transform as matrixes in matrix without pre-calculating the Bessel zeros; Case2: Run the transform as DFT, DHT and IDFT in sequence without pre-calculating the Bessel zeros; Case3: Run the transform as DFT, DHT. Test cases Total computing time (s) Case 1 3,346.5 Case 2 321.1 Case 3 14.3 Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 34/38 http://dx.doi.org/10.7717/peerj-cs.257/fig-25 https://peerj.com/computer-science/ http://dx.doi.org/10.7717/peerj-cs.257 transform sequentially, which demonstrates that the discrete transform shows good precision. Guidelines of choosing sample size As discussed in “Conclusion” and proved by test cases, the sample size N1 (sample size in the radial direction) and N2 (sample size in the angular direction) do not have to be of the same order. For functions with different properties, sample size in the different directions could be very different. To approximate the continuous FT properly, sample size should be chosen based on the discussion in “Conclusion”. Interpretation of the transform By interpreting the transform as a 1DDFT, 1D DHT and 1D IDFT, the computing time of the transform is improved to a useful level since the FFT can be used to compute the DFT. APPENDIX: IMPROVING THE COMPUTING TIME OF THE TRANSFORM One of the advantages of the traditional FT is that the computing speed is fast by using the now well-established fft algorithm. To reduce the computing time of the 2D DFT in polar coordinates, the following steps are recommended: 1. Interpreting the transform as three sequential operations (DFT, DHT, IDFT) instead of a single four-dimensional matrix. 2. Pre-calculating and saving the Bessel zeros. Reducing computing time by interpreting the transform as three operations in sequence As explained above, the essence of the transform is that the matrix fpk is transformed into the matrix Fql. The intuitive way to interpret the transform kernel is to think of it as a four-dimensional matrix or matrices in a matrix. However, interpreting the transform as a 1D-DFT of each column, a 1D-DHT of each row and then a 1D-IDFT of each column makes it possible to use the Matlab built in functions fft and ifft, which significantly reduced the computational time. Reduce computing time by pre-calculating the Bessel Zeros After defining the transform as three operations in sequence and using the matrix for the DHT defined in Lozier (2003), it was found that a lot of computational time was used to calculate the Bessel zeros for every different test case, even though the Bessel zeros are the same in every case. Pre-calculating the Bessel zeros and storing the results for large numbers of N1 and N2 saves a lot of time. Table 7 shows the computing time of a forward transform on the same computer (Processor: Intel(R) Core(TM) i7-4710HQ CPU, RAM:12GB) for three cases: 1. Evaluate the transform as matrices in a matrix without pre-calculating the Bessel zeros. 2. Evaluate the transform as a DFT, DHT and IDFT in sequence without pre-calculating the Bessel zeros. Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 35/38 https://peerj.com/computer-science/ http://dx.doi.org/10.7717/peerj-cs.257 3. Evaluate the transform as a DFT, DHT and IDFT in sequence with pre-calculating the Bessel zeros. The Gaussian function was used as the test function therefore N1 ¼ 383 and N2 ¼ 15. From Table 7, it can be clearly observed that the computing time by running the transform as matrices in a matrix costs 3,346.5 s, which is not acceptable for the transform to be useful. Testing the transform as three operations turns 3,346.5 s into 321.1 s. This is much better. Finally, pre-calculating the Bessel Zeros makes the transform much faster and applicable. ADDITIONAL INFORMATION AND DECLARATIONS Funding This work was financially supported by the Natural Sciences and Engineering Research Council of Canada, grant number RGPIN-2016-04190. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Grant Disclosures The following grant information was disclosed by the authors: Natural Sciences and Engineering Research Council of Canada: RGPIN-2016-04190. Competing Interests The authors declare that they have no competing interests. Author Contributions � Xueyang Yao performed the experiments, analyzed the data, performed the computation work, prepared figures and/or tables, authored or reviewed drafts of the paper, and approved the final draft. � Natalie Baddour conceived and designed the experiments, analyzed the data, authored or reviewed drafts of the paper, and approved the final draft. Data Availability The following information was supplied regarding data availability: Sample Matlab code is available as a Supplemental File. REFERENCES Abbas SA, Sun Q, Foroosh H. 2017. An exact and fast computation of discrete Fourier transform for polar and spherical grid. IEEE Transactions on Signal Processing 65(8):2033–2048 DOI 10.1109/TSP.2016.2645510. Averbuch A, Coifman RR, Donoho DL, Elad M, Israeli M. 2006. Fast and accurate Polar Fourier transform. Applied and Computational Harmonic Analysis 21(2):145–167 DOI 10.1016/j.acha.2005.11.003. Baddour N. 2011. Two-dimensional Fourier transforms in polar coordinates. Advances in Imaging and Electron Physics 165:1–45. Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 36/38 http://dx.doi.org/10.7717/peerj-cs.257#supplemental-information http://dx.doi.org/10.1109/TSP.2016.2645510 http://dx.doi.org/10.1016/j.acha.2005.11.003 https://peerj.com/computer-science/ http://dx.doi.org/10.7717/peerj-cs.257 Baddour N. 2019. The discrete Hankel transform. In: Nikolic G, ed. Fourier Transforms: Century of Digitalization and Increasing Expectations. London: IntechOpen. Baddour N. 2019. Discrete two-dimensional Fourier transform in polar coordinates part I: theory and operational rules. Mathematics 7(8):698 DOI 10.3390/math7080698. Baddour N, Chouinard U. 2015. Theory and operational rules for the discrete Hankel transform. Journal of the Optical Society of America A 32(4):611–622 DOI 10.1364/JOSAA.32.000611. Baddour N, Chouinard U. 2017. Matlab code for the discrete Hankel transform. Journal of Open Research Software 5(1):4 DOI 10.5334/jors.82. Bracewell R. 1999. The Fourier transform and its applications. New York: McGraw-Hill. Cooley JW, Tukey JW. 1965. An algorithm for the machine calculation of complex Fourier series. Mathematics of Computation 19(90):297–301 DOI 10.1090/S0025-5718-1965-0178586-1. Dutt A, Rokhlin V. 1993. Fast Fourier transforms for nonequispaced data. SIAM Journal on Scientific Computing 14(6):1368–1393 DOI 10.1137/0914081. Dutt A, Rokhlin V. 1995. Fast Fourier transforms for nonequispaced data, II. Applied and Computational Harmonic Analysis 2(1):85–100 DOI 10.1006/acha.1995.1007. Fahimian BP, Zhao Y, Huang Z, Fung R, Mao Y, Zhu C, Khatonabadi M, DeMarco JJ, Osher SJ, McNitt-Gray MF, Miao J. 2013. Radiation dose reduction in medical X-ray CT via Fourier- based iterative reconstruction. Medical Physics 40(3):031914 DOI 10.1118/1.4791644. Fenn M, Kunis S, Potts D. 2007. On the computation of the polar FFT. Applied and Computational Harmonic Analysis 22(2):257–263 DOI 10.1016/j.acha.2006.05.009. Fessler JA, Sutton BP. 2003. Nonuniform fast Fourier transforms using min–max interpolation. IEEE Transactions on Signal Processing 51(2):560–574 DOI 10.1109/TSP.2002.807005. Fourmont K. 2003. Non-equispaced fast Fourier transforms with applications to tomography. Journal of Fourier Analysis and Applications 9(5):431–450 DOI 10.1007/s00041-003-0021-1. Guizar-Sicairos M, Gutiérrez-Vega JC. 2004. Computation of quasi-discrete Hankel transforms of integer order for propagating optical wave fields. Journal of the Optical Society of America A 21(1):53–58 DOI 10.1364/JOSAA.21.000053. Higgins W, Munson D Jr. 1987. An algorithm for computing general integer-order Hankel transforms. IEEE Transactions on Acoustics Speech and Signal Processing 35(1):86–97. Lee E, Fahimian BP, Iancu CV, Suloway C, Murphy GE, Wright ER, Castaño-Díez D, Jensen GJ, Miao J. 2008. Radiation dose reduction and image enhancement in biological imaging through equally-sloped tomography. Journal of Structural Biology 164(2):221–227 DOI 10.1016/j.jsb.2008.07.011. Lewitt RM, Matej S. 2003. Overview of methods for image reconstruction from projections in emission computed tomography. Proceedings of the IEEE 91(10):1588–1611 DOI 10.1109/JPROC.2003.817882. Lozier DW. 2003. NIST digital library of mathematical functions. Annals of Mathematics and Artificial Intelligence 38(1–3):105–119 DOI 10.1023/A:1022915830921. Plonka G, Potts D, Steidl G, Tasche M. 2018. Numerical Fourier analysis. Basel: Birkhäuser. Potts D, Steidl G, Tasche M. 2001. Fast Fourier transforms for nonequispaced data: a tutorial. In: Benedetto JJ, Ferreira PJSG, eds. Modern Sampling Theory: Mathematics and Applications. Boston: Birkhäuser Boston, 247–270. Poularikas AD. 2010. Transforms and applications handbook. Third Edition. Boca Raton: CRC Press. Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 37/38 http://dx.doi.org/10.3390/math7080698 http://dx.doi.org/10.1364/JOSAA.32.000611 http://dx.doi.org/10.5334/jors.82 http://dx.doi.org/10.1090/S0025-5718-1965-0178586-1 http://dx.doi.org/10.1137/0914081 http://dx.doi.org/10.1006/acha.1995.1007 http://dx.doi.org/10.1118/1.4791644 http://dx.doi.org/10.1016/j.acha.2006.05.009 http://dx.doi.org/10.1109/TSP.2002.807005 http://dx.doi.org/10.1007/s00041-003-0021-1 http://dx.doi.org/10.1364/JOSAA.21.000053 http://dx.doi.org/10.1016/j.jsb.2008.07.011 http://dx.doi.org/10.1109/JPROC.2003.817882 http://dx.doi.org/10.1023/A:1022915830921 https://peerj.com/computer-science/ http://dx.doi.org/10.7717/peerj-cs.257 Scott MC, Chen C-C, Mecklenburg M, Zhu C, Xu R, Ercius P, Dahmen U, Regan BC, Miao J. 2012. Electron tomography at 2.4-ångström resolution. Nature 483(7390):444 DOI 10.1038/nature10934. Shannon CE. 1984. Communication in the presence of noise. Proceedings of the IEEE 72(9):1192–1201 DOI 10.1109/PROC.1984.12998. Stark H. 1979. Sampling theorems in polar coordinates. Journal of the Optical Society of America 69(11):1519–1525 DOI 10.1364/JOSA.69.001519. Stark H, Wengrovitz M. 1983. Comments and corrections on the use of polar sampling theorems in CT. IEEE Transactions on Acoustics, Speech, and Signal Processing 31(5):1329–1331 DOI 10.1109/TASSP.1983.1164195. Xu Y, Feng D, Wang LV. 2002. Exact frequency-domain reconstruction for thermoacoustic tomography. I. Planar geometry. IEEE Transactions on Medical Imaging 21(7):823–828 DOI 10.1109/TMI.2002.801172. Yao and Baddour (2020), PeerJ Comput. Sci., DOI 10.7717/peerj-cs.257 38/38 http://dx.doi.org/10.1038/nature10934 http://dx.doi.org/10.1109/PROC.1984.12998 http://dx.doi.org/10.1364/JOSA.69.001519 http://dx.doi.org/10.1109/TASSP.1983.1164195 http://dx.doi.org/10.1109/TMI.2002.801172 https://peerj.com/computer-science/ http://dx.doi.org/10.7717/peerj-cs.257 Discrete two dimensional Fourier transform in polar coordinates part II: numerical computation and approximation of the continuous transform ... Introduction Definition of the discrete 2d ft in polar coordinates Discrete transform to approximate the continuous transform Discretization points and sampling grid Numerical computation of the transform Numerical evaluation of the 2d dft in polar coordinates to approximate the continuous ft Summary and conclusion Appendix: improving the computing time of the transform References << /ASCII85EncodePages false /AllowTransparency false /AutoPositionEPSFiles true /AutoRotatePages /None /Binding /Left /CalGrayProfile (Dot Gain 20%) /CalRGBProfile (sRGB IEC61966-2.1) /CalCMYKProfile (U.S. Web Coated \050SWOP\051 v2) /sRGBProfile (sRGB IEC61966-2.1) /CannotEmbedFontPolicy /Warning /CompatibilityLevel 1.4 /CompressObjects /Off /CompressPages true /ConvertImagesToIndexed true /PassThroughJPEGImages true /CreateJobTicket false /DefaultRenderingIntent /Default /DetectBlends true /DetectCurves 0.0000 /ColorConversionStrategy /LeaveColorUnchanged /DoThumbnails false /EmbedAllFonts true /EmbedOpenType false /ParseICCProfilesInComments true /EmbedJobOptions true /DSCReportingLevel 0 /EmitDSCWarnings false /EndPage -1 /ImageMemory 1048576 /LockDistillerParams false /MaxSubsetPct 100 /Optimize true /OPM 1 /ParseDSCComments true /ParseDSCCommentsForDocInfo true /PreserveCopyPage true /PreserveDICMYKValues true /PreserveEPSInfo true /PreserveFlatness true /PreserveHalftoneInfo false /PreserveOPIComments false /PreserveOverprintSettings true /StartPage 1 /SubsetFonts true /TransferFunctionInfo /Apply /UCRandBGInfo /Preserve /UsePrologue false /ColorSettingsFile (None) /AlwaysEmbed [ true ] /NeverEmbed [ true ] /AntiAliasColorImages false /CropColorImages true /ColorImageMinResolution 300 /ColorImageMinResolutionPolicy /OK /DownsampleColorImages false /ColorImageDownsampleType /Average /ColorImageResolution 300 /ColorImageDepth 8 /ColorImageMinDownsampleDepth 1 /ColorImageDownsampleThreshold 1.50000 /EncodeColorImages true /ColorImageFilter /FlateEncode /AutoFilterColorImages false /ColorImageAutoFilterStrategy /JPEG /ColorACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /ColorImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000ColorACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000ColorImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasGrayImages false /CropGrayImages true /GrayImageMinResolution 300 /GrayImageMinResolutionPolicy /OK /DownsampleGrayImages false /GrayImageDownsampleType /Average /GrayImageResolution 300 /GrayImageDepth 8 /GrayImageMinDownsampleDepth 2 /GrayImageDownsampleThreshold 1.50000 /EncodeGrayImages true /GrayImageFilter /FlateEncode /AutoFilterGrayImages false /GrayImageAutoFilterStrategy /JPEG /GrayACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /GrayImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000GrayACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000GrayImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasMonoImages false /CropMonoImages true /MonoImageMinResolution 1200 /MonoImageMinResolutionPolicy /OK /DownsampleMonoImages false /MonoImageDownsampleType /Average /MonoImageResolution 1200 /MonoImageDepth -1 /MonoImageDownsampleThreshold 1.50000 /EncodeMonoImages true /MonoImageFilter /CCITTFaxEncode /MonoImageDict << /K -1 >> /AllowPSXObjects false /CheckCompliance [ /None ] /PDFX1aCheck false /PDFX3Check false /PDFXCompliantPDFOnly false /PDFXNoTrimBoxError true /PDFXTrimBoxToMediaBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXSetBleedBoxToMediaBox true /PDFXBleedBoxToTrimBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXOutputIntentProfile (None) /PDFXOutputConditionIdentifier () /PDFXOutputCondition () /PDFXRegistryName () /PDFXTrapped /False /CreateJDFFile false /Description << /CHS /CHT /DAN /DEU /ESP /FRA /ITA /JPN /KOR /NLD (Gebruik deze instellingen om Adobe PDF-documenten te maken voor kwaliteitsafdrukken op desktopprinters en proofers. De gemaakte PDF-documenten kunnen worden geopend met Acrobat en Adobe Reader 5.0 en hoger.) /NOR /PTB /SUO /SVE /ENU (Use these settings to create Adobe PDF documents for quality printing on desktop printers and proofers. Created PDF documents can be opened with Acrobat and Adobe Reader 5.0 and later.) >> /Namespace [ (Adobe) (Common) (1.0) ] /OtherNamespaces [ << /AsReaderSpreads false /CropImagesToFrames true /ErrorControl /WarnAndContinue /FlattenerIgnoreSpreadOverrides false /IncludeGuidesGrids false /IncludeNonPrinting false /IncludeSlug false /Namespace [ (Adobe) (InDesign) (4.0) ] /OmitPlacedBitmaps false /OmitPlacedEPS false /OmitPlacedPDF false /SimulateOverprint /Legacy >> << /AddBleedMarks false /AddColorBars false /AddCropMarks false /AddPageInfo false /AddRegMarks false /ConvertColors /NoConversion /DestinationProfileName () /DestinationProfileSelector /NA /Downsample16BitImages true /FlattenerPreset << /PresetSelector /MediumResolution >> /FormElements false /GenerateStructure true /IncludeBookmarks false /IncludeHyperlinks false /IncludeInteractive false /IncludeLayers false /IncludeProfiles true /MultimediaHandling /UseObjectSettings /Namespace [ (Adobe) (CreativeSuite) (2.0) ] /PDFXOutputIntentProfileSelector /NA /PreserveEditing true /UntaggedCMYKHandling /LeaveUntagged /UntaggedRGBHandling /LeaveUntagged /UseDocumentBleed false >> ] >> setdistillerparams << /HWResolution [2400 2400] /PageSize [612.000 792.000] >> setpagedevice