Discrete Two-Dimensional Fourier Transform in Polar Coordinates Part I: Theory and Operational Rules

: 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 this paper, we propose and evaluate the theory of the 2D discrete Fourier transform (DFT) in polar coordinates. This discrete theory is shown to arise from discretization schemes that have been previously employed with the 1D DFT and the discrete Hankel transform (DHT). The proposed transform possesses orthogonality properties, which leads to invertibility of the transform. In the ﬁrst part of this two-part paper, the theory of the actual manipulated quantities is shown, including the standard set of shift, modulation, multiplication, and convolution rules. Parseval and modiﬁed Parseval relationships are shown, depending on which choice of kernel is used. Similar to its continuous counterpart, the 2D DFT in polar coordinates is shown to consist of a 1D DFT, DHT and 1D inverse DFT.


Introduction
The Fourier transform (FT) in continuous and discrete forms has seen much application in various disciplines [1]. It easily expands to multiple dimensions, with all the same rules of the one-dimensional (1D) case carrying into the multiple dimensions. Recent work has developed the complete toolkit for working with the continuous multidimensional Fourier transform in two-dimensional (2D) polar and three-dimensional (3D) spherical polar coordinates [2][3][4]. However, to date no discrete version of the 2D Fourier transform exists in polar coordinates. Hence, the aim of this paper is to develop the discrete version of the 2D Fourier transform in polar coordinates.
For the discrete version of the transform, the values of the transform will be available only at discrete points. To quote Bracewell [5], "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". This paper thus aims to develop the mathematical theory of the discrete two-dimensional Fourier transform in polar coordinates. Standard 'operational rules' associated with any Fourier transform (shift, modulation, multiplication, and convolution) will be developed. Parseval and modified Parseval relationships will also be shown, depending on the choice of kernel used.
To the best of the author's knowledge, there is no discrete version of the 2D Fourier transform in polar coordinates. It was shown in [2,4] that the 2D continuous Fourier transform in polar coordinates is actually a combination of a single dimensional Fourier transform, a Hankel transform, followed by an inverse Fourier transform. Of course, the discrete version of the 1D standard Fourier transform is very well known and the literature on this subject alone is vast. Recently, a discrete version of the Hankel transform has been proposed [6,7], yet this discrete transform is still in one dimension. We will show further on that the 2D Fourier transform in polar coordinates requires this transform.
Other researchers have defined the idea of a polar Fourier transform (polar FT), in which the original function in the spatial domain is in Cartesian coordinates but its FT is computed in polar coordinates, meaning discrete polar Fourier data and Cartesian spatial data [8][9][10]. Fast Fourier transforms (FFT) have also been developed for non-equispaced data, referred to as a unequally spaced FFT (USFFT) or non-uniform FFT (NUFFT) [11][12][13][14][15]. Using this approach, frequencies in a polar frequency domain can be considered to be unequally spaced and hence the problem of evaluating a polar FT can be considered as a special case of the USFFT. Averbuch et al. [8] compared the accuracy results of their proposed approach which used a pseudo-polar grid to those obtained by an USFFT approach and demonstrated that their approach show marked advantage over the USFFT. Fenn et al. [10] examined computing the FT on a polar, modified and pseudo-polar grid using the NUFFT, for both forward and backwards transforms. They demonstrated that the NUFFT was effective at this computation. Although the above 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 Fourier transform in polar coordinates that exists as a transform in its own right, with its own set of rules of the actual manipulated quantities.
The outline of the paper is as follows. Section 2 presents some of the necessary background material. Section 3 introduces an intuitive 'motivation' for the definition of the 2D Discrete Fourier Transform (DFT) in polar coordinates that will be introduced by considering space and band-limited functions. This leads to an intuitive discretization scheme and an intuitive kernel for the proposed 2D DFT, which is introduced in Section 4. Section 5 introduces the proposed transform while Section 6 derives the transform properties including modulation, shift, multiplication and convolution rules. Section 7 discusses Parseval relations while Section 8 demonstrates that the proposed transform can indeed be decomposed a sequence of DFT, Discrete Hankel Transform (DHT) and inverse DFT (IDFT), in keeping with the approach of the continuous version of the transform. Finally, Section 8 concludes the paper.

Background: Continuous 2D Fourier Transforms in Polar Coordinates
The 2D Fourier transform of a function f ( → r ) = f (x, y) expressed in 2D Cartesian coordinates is defined as [4]: The inverse Fourier transform is given by: where the shorthand notation of → ω = ω x , ω y , → r = (x, y) has been used. For functions with cylindrical or circular symmetry, it is often more convenient to express both the original function f ( → r ) and its 2D Fourier transform F → ω in polar coordinates. If so, polar coordinates can be introduced as x = r cos θ, y = r sin θ and similarly in the spatial frequency domain as ω x = ρ cos ψ and ω y = ρ sin ψ, otherwise written as, r 2 = x 2 + y 2 , θ = arctan(y/x) and ρ 2 = ω 2 x + ω 2 y , ψ = arctan ω y /ω x . Given a function in polar coordinates f (r, θ), where θ is the angular variable and r is the radial variable, the function can be expanded into a Fourier series as: where the Fourier coefficients are given by: Similarly, the 2D Fourier transform of f (r, θ) is given by F(ρ, ψ). The function F(ρ, ψ), where ψ is the angular frequency variable and ρ is the radial frequency variable, can also be expanded into a Fourier series as: where: We note that F n (ρ) is NOT the Fourier transform of f n (r). The development details can be found in [4], where it is demonstrated that the relationship is given by: where H n {·} denotes an nth order Hankel transform [3], see Appendix A.1 [3]. The inverse relationship is given by: Thus, the nth term in the Fourier series of the original function will Hankel transform into the nth term of the Fourier series of the Fourier transform function via an nth order Hankel transform for the nth term. Therefore, the steps for finding the 2D Fourier transform F(ρ, ψ) of a function f (r, θ) are (i) finding its Fourier series coefficients in the angular variable f n (r), Equation (4), (ii) finding the Fourier series coefficient of the Fourier transform, F n (ρ) via F n (ρ) = 2π i −n H n f n (r) , then (iii) taking the inverse Fourier series transform (summing the series) with respect to the frequency angular variable, Equation (5).
The discrete equivalent to the relationships given by Equations (3) to (8) have not been developed and it is the goal of this paper to develop the discrete counterparts of these equations.

Space-Limited Functions
To motivate the discrete version of a 2D Fourier transform in polar coordinates, we follow the same path used to derive the classical discrete Fourier transform (DFT) and also the recently-proposed discrete Hankel transform (DHT) [6]. This approach starts with a space (or time for the traditional FT) limited function in one domain and then makes the assumption that the transform of the function is also limited in the corresponding frequency domain. While strictly speaking, functions cannot be limited in both space and spatial frequency domains, in practice, they can be made 'effectively' limited in the domain where they are not exactly limited by suitable truncation of an appropriate series. This is how the DFT and DHT were both motivated. The discrete transforms derived in this manner then have properties that exist in their own right, independent of their ability to approximate their continuous transform counterpart.
The same path is followed here. A function f (r, θ) in polar coordinates, where θ is the angular variable and r is the radial variable, is expanded into a Fourier series given by Equation (3) where f n (r) is given by Equation (4). It is now supposed that the function f (r, θ) is space-limited, meaning that f (r, θ) and, by virtue of Equation (4), all the Fourier coefficients f n (r) are zero for r ≥ R. Then, it follows that each of the Fourier coefficients f n (r) can be written in terms of a Fourier Bessel series (see [6] and Appendix A.2) as: where the order, n, of the Bessel function in (9) matches the order f n of the Fourier coefficient, C f nk denotes the kth coefficient of the Fourier-Bessel expansion of f n (r) and denotes the kth zero of the nth Bessel function. The C f nk can be found from [16]: Equation (7) gives the relationship between the Fourier coefficients of the function itself and its 2D Fourier transform. Using Equation (7) and making use of the space limited nature of f n (r), Equation (10) can be written as: Therefore, for r < R, Equation (9) becomes: Equation (12) with its infinite summation is exact. Now, evaluating Equation (12) at r = r nk = j nk R j nN 1 for any N 1 and where k < N 1 gives: For k < N 1 , then r nk = j nk R j nN 1 < R, and Equation (13), summing over infinite m, is still exact. For k ≥ N 1 , then r nk = j nk R j nN 1 ≥ R and by the assumption of the space-limited nature of the function, f (r nk ) = 0 for We now assume that the function is also effectively band limited, in addition to being space-limited. Now, a function cannot be finite in both space and spatial frequency (equivalently if using a standard Fourier transform it cannot be finite in both time and frequency). However, if a function is effectively band-limited, then there exists an integer N 1 for which F n j nm R ≈ 0 for m > N 1 . In other words, an interval can be found beyond which the Fourier transform coefficients F n (ρ) become very small. Since the convergence of the Fourier-Bessel series in (13) is known, then lim m→∞ F n j nm R = 0. In other words, for any arbitrarily small ρ, there exists an integer N 1 for which F n j nm R < ρ for m > N 1 . Hence, using notion of an effective band-limit as stated in the preceding paragraph, the series in Equation (13) can be terminated at a suitably chosen N 1 , thus giving an effective band limit. Termination of the series at m = N 1 is equivalent to assuming that F n (ρ) ≈ 0 for ρ > W ρ = j nN 1 R . It is noted that at m = N 1 , the last term in Equation (13) is J n j nN 1 j nk j nN 1 = J n ( j nk ) = 0, so that termination of the series at N 1 implies that Equation (13) becomes Equation (14) is the discrete equivalent of Equation (8) in that it demonstrates that the relationship between discrete samples of f n (r) and F n (r) is given by a discrete Hankel transform type of relationship, whereas the continuous relationship involved a continuous Hankel transform. The termination of the series at N 1 is equivalent to assuming an "effective" band-limit on the function. In other words, it states that for m > N 1 , the values of F n j nm R , which from Equation (11) are proportional to the Fourier-Bessel coefficients, are negligibly small. Of course, this is never exactly true, however, since the Fourier-Bessel series converges, it is always possible to choose N 1 so that the approximation introduced by truncating the series at N 1 is good [16].
The truncation of the series at N 1 also permits Equation (14) to be easily inverted. Multiplying both sides of (14) by 4J n j nk jnp j nN 1 n+1 ( j nk ) and summing over k gives: where we have used the discrete orthogonality of the Bessel functions as given in Appendix A.4. Hence, Equations (14) and (16) offer the basic structure on which to base the discrete transform formulation. Equation (16) is the basic structure to define the forward transform and Equation (14) offers the basic structure to define the inverse transform.
To proceed further, we need ways to compute f n j nk R j nN 1 and F n j nm R . Here, the theory of discrete Fourier transforms can be used. For n ∈ [−M, M] where N 2 = 2M + 1, it is shown in [17] that the Fourier coefficients f n (r) and F n (ρ) can be well approximated with expressions given by (see Appendix A.5): Hence, we will use Equation (17) to write: Equation (18) is a key assumption of the development. Note that in both cases, the function is sampled in the summation over p at the radial variable j pk,m R , that is, it is included in the summation index. However, the function on the left hand side of Equation (18) is sampled at j nm R We show in Appendix A.6 that this assumption is valid. This assumption is what also permits the invertibility of the discrete transforms, since without this assumption it would not be possible to propose an invertible, orthogonal discrete transform. Equation (18) will be used to derive the forward and inverse discrete transforms.

Forward Transform
For the forward transform, we can start with Equation (16), and use the key relationships given by Equation (18). Under these conditions, Equation (16) becomes: Equation (19) is the discrete equivalent of Equation (7). From Equation (19), multiply both sides by e Interchanging the order of summation on the left hand side of (20) and using the orthogonality relationship of the complex exponential (Appendix A.3) gives:

Inverse Transform
For the inverse transform, we start with the structure of Equation (14) and then use the key approximations given in Equation (18) to obtain: Multiplying both sides of Equation (22) by e +i 2πnp N 2 , summing from n = −M..M, interchanging the order of summation on the left hand side and using the orthogonality relationship of the discrete complex exponential gives:

Band-Limited Functions
The process in the previous section can be repeated by starting with the assumption that the function is band-limited. That is, we suppose that the 2D Fourier transform F(ρ, ψ) of f (r, θ) is band-limited, meaning that F(ρ, ψ) itself and therefore by virtue of the equivalent of Equation (9), all of its Fourier coefficients F n (ρ) are zero for ρ ≥ W ρ = 2πW. Typically, W would be given in units of Hz (cycles per second) if using temporal units, or cycles per meter if using spatial units. Hence, the definition of W ρ (with a multiplication by 2π) ensures that the final units are given in s −1 or m −1 . The details of this development follow the same steps as for the space-limited function but start with the assumption of a band-limited function and then impose a space-limit (i.e., truncation of the series). The results of this are summarized below.

Summary of Above Relationships
From the above, we summarize the derived relationships. In the case of a space-limited function, it is found that the forward transform is given by: and the inverse transform is given by: Similarly, starting from the assumption of a bandlimited function, the forward transform is given by: and the inverse transform is given by: It is noted that the forward-inverse transform pair defined by Equations (24) and (25) is similar to the transform pair defined by (26) and (27), with a few differences. First, the sampling points appear to be slightly different, depending on whether we started with the assumption of a space-limited function or a bandlimited function. The second observation is that the form of the transform itself might appear to be slightly different, depending on whether a space-limited or a band-limited function was assumed as a starting point. However, it was shown in [6] that for a nth order discrete Hankel transform, the required relationship between the band limit and space limit is given by W ρ R = j nN 1 . If the substitution W ρ R = j nN 1 is used in Equations (24) and (25), then it yields the same discrete transform as the transform pair defined by (26) and (27). Also, the relationship W ρ R = j nN 1 arose naturally in the development above when the truncation of the Fourier-Bessel series at N 1 was implemented, meaning that the truncation of the series at N 1 is the same as assuming W ρ R = j nN 1 .
Formally using the relationship W ρ R = j nN 1 , the expressions in Equations (24) and (25) can also be written using a symmetric forward/inverse transform pair, where the forward transform is given by: For the inverse transform, we can similarly write: The advantage of the formulation in Equations (28) and (29) shall be noted in the next section in that it suggests a symmetric form of the kernel for the 2D discrete transform in polar coordinates.
The above demonstrates that a natural, (N 1 − 1) × N 2 dimensional discretization scheme in finite space and finite frequency space is given by: and: where p, k, q, m, n, N 1 , and R can be used to formally switch from a finite frequency domain to a finite space domain. This is a 'formal' approach because in making this substitution, the index of the Bessel function is not fixed whereas W ρ and R are assumed fixed values. Nevertheless, it demonstrates the approach to switching from a space-limited based discretization scheme to a band-limited discretization scheme.

Proposed Kernel for 2D Polar Discrete Fourier Transform
To work with the polar 2D DFT, a kernel for the transformation is required. Inspired by the formulations shown in Equations (24) and (25), we propose the following kernels: where p, k, q, m, n, N 1 , and It is noted that the proposed kernels in Equation (32) are almost complex conjugates of each other save for a factor of j 2 nN 1 in the denominator of E − qm;pk . The formulation in Equation (32) is proposed in order to emulate Equations (24) and (25). A symmetric formulation of the kernels, with one j nN 1 in the denominator of each of E − (qm; pk) and E + (qm; pk) would also be possible and would make E ± qm;pk complex conjugates of each other; however, such a kernel would be more of a departure from a discretization of the continuous transform. The integers N 1 , and N 2 denote the size of the working spaces, with N 2 giving the size in the angular direction and N 1 giving the size in the radial direction. Since N 2 = 2M + 1, it follows that N 2 must be an odd integer. The notation for E − (qm; pk) and E + (qm; pk) are chosen deliberately. The subscript (+ or -) indicate the sign on the i ± and on the exponent containing the p variable; the q variable exponent then takes the opposite sign.

Another Choice of Kernel
A second, more symmetric choice of kernel is also possible. We will see that this choice of kernel will allow for a more traditional version of Parseval's theorem. All the following expressions will hold with either form of kernel. Using as inspiration the forms written in Equations (28) and (29), then we suggest for a kernel the following expression: As before, p, k, q, m, n, N 1 , and qm;pk , as mentioned above.

Orthogonality of the Proposed Kernel
In what follows, we assume the ranges of the variables are such that p, k, q, m, n, N 1 , and N 2 are We state and prove that the following relationship is true: where δ pp is the Kronecker-delta function, defined as δ pp = 1 if p = p and δ pp = 0 otherwise. It is known that the continuous complex exponential expression can be written as: Hence, the form of the discrete kernel as proposed in (32) or (33) resembles discrete samples of the right hand side of Equation (35). It then follows that our proposed kernels in (32) or (33) can be considered to be the (discrete) corresponding form of the complex exponential kernel for the proposed discrete transform. The orthogonality relationship in (34) can then be considered to be the discrete version of: where the integration over the frequency vector → ω has been replaced with a discrete sum over the frequency vector indices (q,m). The proof of Equation (34) uses the orthogonality of the discrete complex exponential and the discrete Hankel transform and can be found in Appendix A.7.
It can be similarly shown that the following orthogonality relationship is also true: which is similarly to be considered to be the discrete version of: Once again, the integration over the vector → r has been replaced with a discrete sum over the → r vector indices (p,k). The proof of Equation (37) can also be found in Appendix A.7. The orthogonality expressions in Equations (34) and (37) still hold if E ± qm;pk is replaced with the symmetric E (s)± qm;pk since the only difference between the E ± qm;pk and E (s)± qm;pk is the attribution of a j nN 1 term in the denominator and this makes no difference when the two kernels are multiplied.

Proposed Transform
In this section, we propose a definition of the 2D discrete Fourier transform (DFT) in polar coordinates which is motivated by the results of the 2D Fourier transform applied to space-limited and band-limited functions and also by the proposed kernel. The 2D DFT in polar coordinates will be a transform that transforms a 2-subscript set of numbers (ie matrix) f pk to another set of values, matrix F qm where p, k, q, m, are integers such that 1 ≤ m, k, ≤ N 1 − 1 and −M ≤ p, q ≤ M where N 2 = 2M + 1 for integers N 1 , and N 2 .

Forward and Inverse Transform
The proposed forward transform, f pk → F qm is given by: where N 2 = 2M + 1 for some integer M. Similarly, for the inverse transform we propose: In the proposed transform, E (s)± qm;pk could easily be used in placed of E ± qm;pk and all the following expressions will still be valid.
Proof. Substituting Equation (39) into the right-hand side of (40), interchanging the order of summation and using the orthogonality relationships of the kernel given in Equation (34) gives: Similarly, substituting Equation (40) into the right-hand side of (39), interchanging the order of summation and using the orthogonality of the kernel given in Equation (37) gives: Hence, Equation (39) and (40) are inverses of each other. These expressions would also hold if E (s)± qm;pk were used instead of E ± qm;pk .

The Complex Exponential
For the discrete case, the functions E − qm;pk and E + qm;pk as introduced above are the complex exponentials for this space, satisfying the required orthogonality condition and functioning as the kernel for the 2D-DFT in polar coordinates. These kernels are not e ±i → ω· → r evaluated at particular points because the evaluation of the discrete radial variables in regular and frequency space varies with the order of the Bessel function. Nevertheless, these functions are the 'effective' complex exponentials for the space under consideration. From the orthogonality condition of the 2D polar DFT kernel, it can be shown that the expected Fourier rule of a complex exponential transforming to a delta function applies. Specifically, the 2D DFT of f pk = E + q 0 m 0 ;pk for some fixed, given values (q 0 , m 0 ) is given by: Hence, f pk = E + q 0 m 0 ;pk transforms to δ qq 0 δ mm 0 or in compact notation, E + q 0 m 0 ;pk ⇔ δ qq 0 δ mm 0 . This is the discrete version of the transform of exp → ω 0 · → r .

The Delta Function
Clearly, the discrete equivalent of the Dirac-delta function is the Kronecker-delta function and in 2D, this needs to be a 2-subscript function. Thus, the discrete function whose 2D DFT is sought is given by f pk = δ pp 0 δ kk 0 , which defines a matrix indexed by (p, k) where all the entries are zero except for the index where p = p 0 and k = k 0 . The dimensions of this matrix are in keeping with all the dimensions assumed for the space which are p, k, q, m, n, N 1 , and N 2 are integers such that −M ≤ n ≤ M, 1 ≤ m, k, ≤ N 1 − 1 and −M ≤ p, q ≤ M and where 2M + 1 = N 2 . Finding the 2D DFT of this function gives: Hence, as in the continuous case, the delta function transforms to the complex exponential (with a negative sign in the exponent). Hence we have another the Fourier pair δ pp 0 δ kk 0 ⇔ E − qm;p 0 k 0 .

The Generalized Shift Operator
For a one dimensional Fourier transform, the shift rule is one of the known transform rules. This rule says that a shift in time is equivalent to a modulation in frequency. Mathematically, this is stated as: Using this result as motivation, a generalized-shift operator is defined by finding the inverse DFT of the DFT of the function multiplied by the DFT kernel (modulation). A generalized shift operator was first proposed by Levitan [18], and our definition is a discretized version of this definition. Levitan suggested the complex conjugate of the Fourier operator as a generalized shift operator, which for Fourier transforms is the inverse transform operator. This approach to a generalized shift operator has previously been used with the Hankel transform itself [6,19]. Thus, we define the definition of a generalized-shifted function f p 0 k 0 pk as the inverse Fourier transform of the function multiplied by the inverse transform operator. That is, it is defined as: Here, f pk is the original (unshifted) function with 2D DFT F qm such that f pk → F qm . f p 0 k 0 pk is the shifted function where p 0 k 0 denotes the amount of the shift (the equivalent of a in Equation (45)).
The shifted function f p 0 k 0 pk can also be expressed in terms of the unshifted function f pk by writing F qm in terms of f pk such as: By interchanging the order of summation, this can be rewritten as: Equation (48) permits the definition of a shift operator so that the shift operator in the spatial domain is defined as: This triple-product shift operator resembles previous definitions of shift operators for multidimensional Fourier transforms [2,3], generalized Hankel convolutions [20][21][22] and also discrete Hankel transforms [6].

Forward Transform of the Generalized Shift
We now consider the forward 2D Fourier transform of the generalized shifted function f p 0 k 0 pk . From the definition of the shifted function given in Equation (46), it is obvious that the forward transform of the shifted function is given by: The above can also be verified directly. The 2D Fourier transform of the shifted function can be found from: where the definition in (46) was used. Interchanging the order of summation and using the orthogonality result in (37) gives: This gives another transform pair and also defines the shift-modulation rule. This rule is in analogy with the shift-modulation rule for regular Fourier transforms that states that a shift in the spatial/time domain is equivalent to modulation in the frequency domain:

Modulation
We suppose that the forward 2D-DHT of a function g pk is 'modulated' in the space domain so that the function whose transform we seek is f pk = E + q 0 m 0 ;pk g pk . This is the discrete equivalent of a function g(t) modulated as e iat g(t) Here, the interpretation of f pk = E + q 0 m 0 ;pk g pk is as follows: f pk = E + q 0 m 0 ;pk g pk f pk = g pk Again, we implement the definition of the forward transform on the modulated function f pk = E + q 0 m 0 ;pk g pk so that: and write g pk in terms of its inverse transform: So that Equation (56) becomes: Interchanging the order of summation gives: By comparing Equation (59) with Equation (49), we recognize the shift operator as shown in (59). This follows from a shift over the (q,m) variables and defines a shift operator in the frequency domain as: Hence, Equation (59) can be written as: The shift operator in the frequency domain over the (q,m) variables as given by Equation (60) can be compared to the shift operator over the (p,k) variables in the space domain as shown in (49). We note that operations in the spatial domains are operations that involve the (p,k) variables or the second group of variables in E ± qm;pk . Similarly, operations in the frequency domain involve operations over the (q,m) variables or the first set of variables in E ± qm;pk . Hence, the above development shows the derivation of a modulation-shift rule, where the forward 2D-DHT of a modulated function is equivalent to a generalized shift in the frequency domain. This gives the following transform pair: Otherwise stated, Equation (62) shows that modulation in the space domain is equivalent to shift in the frequency domain, in keeping with expectations for a (generalized) Fourier transform.

Convolution-Multiplication
For a 2D convolution/multiplication rule, we consider a 2D convolution in the space domain. The convolution is defined in the traditional manner as the product of a shifted function with another unshifted function, and then the summation over all possible shifts. Specifically, we write it as: summation over all possible shifts h p 0 k 0 pk shifted function where h p 0 k 0 pk is the h pk shifted by p 0 k 0 given by: The summation in Equation (63) is then over all the possible shifts. Taking the forward transform of f pk as defined in (63) gives: Interchanging the order of summation so that the summation over p,k is performed first and using the orthogonality of the kernel gives: Now summing over p 0 , k 0 and again using the orthogonality of the kernel gives: In other words, we have the result that: Equation (68) is, of course, the expected convolution-multiplication rule where convolution in the space domain is equivalent to multiplication in the frequency domain.

Multiplication-Convolution Rule
We now consider the forward 2D FT of a term-by-term product in the space domain so that f pk = h pk g pk . Then, the forward transform of the term-by-term product is given by: Using the definitions of the inverse 2D FT to write h pk and g pk then: In Equation (70), we have used the modulation rule E + q 0 m 0 ;pk g pk ⇔ G q 0 m 0 qm . In other words, Equation (70) states that: Hence, h pk g pk ⇔ H qm * * G qm which is the multiplication-convolution rule where multiplication in the space domain is equivalent to convolution in the frequency domain.

Rotation
It is generally known that rotating a function in 2D space also rotates its 2D Fourier transform. We demonstrate that this is still true with our definition of the discrete 2D DFT in polar coordinates. To see this, we consider a shift of the function in frequency space, meaning consider F (q−q 0 )m where a shift by q 0 in the angular coordinate has been implemented. In this case, since the circular direction is circularly periodic, we interpret q − q 0 in the sense of modulo N 2 . So consider the inverse discrete 2D DFT of F (q−q 0 )m , that is from the definition in Equation (40) Now suppose that q = q − q 0 so that q = q + q 0 and q = −M implies q = −q 0 − M and also q = +M implies q = −q 0 + M. Hence, Equation (72) becomes: But because of the circular (N 2 ) periodicity of the function, then: Hence, Equation (72) becomes: As above, f (p−q 0 )k is to be interpreted in the sense of modulo N 2 . However, what this clearly demonstrates is that rotating the Fourier transform by q 0 is equivalent to rotating the original function by q 0 , as is expected of a 2D Fourier transform.

Generalized Parseval Theorem
Under the proposed transform, inner products are preserved and, therefore, energies are preserved with the symmetric version of the transform. With the non-symmetric version of the transform, a modified version of Parseval's theorem is possible. This will be demonstrated in the following subsections.

Parseval's Theorem with the Symmetric Kernel
Consider the total energy of the term-by-term product (Hadamard product) of two matrices in the spatial domain f pk = h pk g pk . We use the overbar notation to denote the complex conjugate, so that g pk denotes the complex conjugate of g pk . We recall that in the case of the symmetric kernel, the complex conjugate of E (s)+ qm;pk is E (s)− qm;pk , which is what will enable the Parseval relationship to exist in its expected form, as will be shown. More specifically, it is noted that: However, Hence, Equation (76) becomes: For the special case that g pk = h pk then Equation (78) yields: Equations (78) and (79) are the expected for of the Parseval relationship, which essentially states that the energy computed in one domain is equivalent to the energy computed in the other domain. The reader is reminded that the symmetric kernel was used for the derivation in (79).

Parseval's Theorem with the Non-Symmetric Kernel
For the non-symmetric kernel, some modifications to the above Parseval relationship are necessary. Again, we consider the total energy of a Hadamard product of two matrices in the spatial domain. However, now we need to define a more 'general' version of a complex conjugate expression in order for the Parseval relationship to exist. We denote this more general version as g pk * (over bar and star) and define this expression as: g pk * := We note in Equation (80) that g pk * uses E − q m ;pk instead of E + q m ;pk (where the latter would normally be used for the complex conjugate). The reason for this is that with the non-symmetric kernel, using E + q m ;pk will not lead to the required orthogonality condition. However, with our 'modified' version of the complex conjugate as denoted by the g pk * and defined in Equation (80), it then follows that: Using the orthogonality of the kernel, Equation (81) becomes: Similarly, we can consider the special product in the frequency domain F qm = H qm G qm * where again the special expression G qm * needs to be defined as follows: Consider: Interchanging the order of summation and summing over the (q,m) variables first gives: Using the orthogonality of the kernel, the last line can be rewritten as: In summary, Equation (82) shows how to interpret H q m G q m and Equation (86) shows how to interpret h p k g p k and also shows that they are not quite equivalent as was the case for the symmetric kernel. In summary, In the special case that h = g, then Equation (87) becomes:

Discussion: Interpretation of the Transform
In the previous sections, we demonstrated that the 2D DFT in polar coordinates is most conveniently defined in terms of the kernels E ± qm;pk or E (s)± qm;pk , and indeed this definition allows many of the proofs of the DFT properties to assume a straightforward form that exploits the properties of the kernel. In this section, we demonstrate that the proposed forms of the 2D DFT can be interpreted in terms of a sequence of 1D DFT, DHT and IDFT discrete transforms, thereby demonstrating that the proposed transform follows the same path as the continuous 2D transform in that it can be decomposed into a sequence of Fourier, Hankel and inverse Fourier transforms [2].

Interpretation of the 2D Forward DFT in Polar Coordinates
Let us reconsider the definition of the forward 2D DFT, Equation (39), and rewrite it as: We can consider these as a sequence of 1D discrete Fourier transforms, with a discrete Hankel transform, as explained in the following. The first step is a forward 1D DFT transforming f pk → f nk where the p subscript is transformed to the n subscript as: The tilde is used to indicate a standard 1D DFT. In matrix terms, this says that each column of f pk is DFT'ed to yield f nk . The second step of Equation (89) is a discrete Hankel transform of order n that transforms f nk →ˆ f nm , where the k subscript is Hankel transformed to the m subscript via: The overhat denotes a DHT, as defined in [6]. Using the transformation matrix notation defined in [6], we define: Hence Equation (91) can be written aŝ In matrix terms, this shows that each row of f nk is nth-order DHT'ed to yieldˆ f nm . The nth row is nth order DHT'ed (with some loose interpretation of row counters since in this case the index n takes on negative values). A scaling operation then gives the Fourier coefficients of the 2D DFTˆ f nm → F nm such that: It is noted that Equation (94) exactly parallels the equivalent step of the continuous form of the transform where F n (ρ) = 2π i −n H n f n (r) , see [4,6]. If the symmetric form of the kernel is used, that is, The final step to compute the forward 2D DFT in polar coordinates is then a standard inverse 1D DFT. Here, each column of F nm → F qm is transformed so that the n subscript is (inverse) transformed to the q subscript via This last step is a 1D IDFT for each column of F nm to obtain F qm . It was shown in [2,4] that a continuous 2D Fourier transform in polar coordinates is a sequence of operations consisting of (i) a Fourier series transform (transforming the continuous function to its discrete set of Fourier coefficients), (ii) a Hankel transform for each Fourier coefficient (an nth order transform for the nth coefficient), and (iii) an inverse Fourier series transform (a set of Fourier coefficients is transformed back to a continuous function via the infinite Fourier series summation). Hence, we have shown here that the proposed 2D DFT in polar coordinates consists of the same sequence of transforms: a forward DFT, a forward DHT and then an inverse DFT.

Interpretation of the 2D Inverse DFT in Polar Coordinates
Similarly, we can decompose the inverse 2D DFT in polar coordinates, from Equation (40) written as: The steps of the inverse 2D DFT are the reverse of those outlined for the forward 2DDFT. First F qm → F nm via a forward 1D DFT: This is followed by a discrete Hankel transform to obtain F nm →ˆ F nk The next step is a scaling operation to obtainˆ F nk → f nk via: The step in Equation (99) follows the pattern of the continuous form transform, where f n (r) = i n 2π H n F n (ρ) , see [2,4,6]. As before, if the symmetric form of the kernel is used (Equation (33)), then Equation (99) is replaced with f nk = i +nˆ F nk . Finally, an inverse 1D DFT is used to obtain f nk → f pk via: As previously mentioned, this parallels the steps taken for the continuous case, with each continuous operation (Fourier series, Hankel transform) replaced by its discrete counterpart (DFT, DHT).
For both forward and inverse 2D DFT, the same sequence of steps are followed. The operations are a 1D DFT of each column of the given matrix, then a DHT of each row, then a term-by-term scaling, and finally an IDFT of each column.

Conclusions
In this paper, a discrete 2D Fourier transform in polar coordinates was motivated and proposed by applying a discretization and truncation approach to the continuous 2D Fourier transform in polar coordinates. This new transform stands in its own right and, unlike previous approaches to a polar FT, is not an evaluation of the Cartesian form of the transform on a polar grid. This approach yields two possible kernels for the discrete 2D transform in polar coordinates. One of these two kernels is closer to the continuous version of the transform and the second kernel is symmetric, in that the kernel for the forward transform is the complex conjugate of the kernel for the inverse transform. Both versions of the kernel yield a 2D transform that transform a 2-subscripted entity (matrix) to another one. The standard set of shift, modulation, multiplication and convolution rules were derived for both kernels and are the same for either form of the kernel. However, only the symmetric kernel yields the expected Parseval relationship. It was also shown that the 2D discrete transform can be interpreted as a 1D discrete Fourier transform (DFT), followed by a 1D discrete Hankel transform (DHT), followed by a 1D inverse DFT. This DFT-DHT-IDFT pattern mimics the manner in which the continuous 2D Fourier transform in polar coordinates is evaluated. In conclusion, part I of the paper proposes the form of the 2D DFT in polar coordinates, and demonstrates the expected operational rules for this transform. Part II of the paper will examine how the proposed 2D DFT in polar coordinates can be used to approximate the continuous FT at certain discrete points.

Conflicts of Interest:
The author declares no conflict of interest.

Appendix A
The following appendices give important definitions for Hankel transforms (Appendix A.1), Fourier Bessel series (Appendix A.2), finite Fourier transforms (Appendix A.5) and also contain statements of the orthogonality of the discrete complex exponential (Appendix A.3) and the discrete Bessel functions (Appendix A.4). Section Appendix A.6 contains a discussion on the sampling points and how they affect the proposed evaluation of the discrete Fourier coefficients at the chosen sampling points. Proofs of the orthogonality of the proposed kernel can be found in Appendix A.7.
where the Fourier coefficients are given by: A principal application of the finite Fourier transform (FFT) is to approximately compute samples of the Fourier transform of a function. We define the FFT partial sum of the samples f 2πp N 2 of the continuous function f (θ) as: where N 2 is an integer such that N 2 = 2M + 1 for some other integer M. The over square-hat notation f indicates the taking of a (finite) Fourier transform. Clearly, Equation (A13) is a Riemann sum for the integral in (A12). It is generally asserted in the signal processing literature that f n ≈ f n , and it is specifically shown in [17] that f n provides a uniformly good estimate for f n for n ∈ [−M, M].
It is also shown in [17] that the finite Fourier transform partial sum given by: is almost as good an approximation to f (θ) as the usual partial sum: Appendix A.6. Sampling Points In this section, the difference between including the radial sampling points in the index of summation for the discrete Fourier transform is discussed. We noted above in Equation (18) that the radial sampling point is included in the index of summation of the discrete Fourier transform. In other words, we wrote for n ∈ [−M, M] that: However, strictly speaking, the radial sampling points should be fixed to the value of the radial sampling point on the left hand side, that is the expected discrete definition of F n j nl R should be given by: Note that in both Equations (A16) and (A17), the index of summation is p, and the radial sampling point is j pl in (A16) but j nl in (A17).
Which of the definitions for F n j nl R is correct? Definition (i), as given in Equation (A16), or definition (ii) as given in Equation (A17)? Traditionally, (ii) of Equation (A17) would be expected but taking this form does not allow the 2D discrete transform that ensues to be invertible. We showed above in the main text of the manuscript, that version (i) with Equation (A16) leads to an invertible, discrete 2D transform. We show in this section that if we confine ourselves to the chosen sampling points, then both versions are equivalent.
Considering the reconstruction formula based on (A14) which says: Then, sampling at θ = 2πm N 2 gives: So now consider the right hand side of Equation (A19) under the two different sampling assumptions implied by (i) or (ii). That is: Therefore, the (ii) version also works the way it is expected to work. Therefore, both (i) and (ii) work properly.
However, if we try evaluating at different values of angular position, say θ = 2πr N 2 where now the sampling index on the angle and the Bessel function do not match, in other words: In this case, (ii) does not yield the expected result, but (i) does. So the question of (i) vs (ii) becomes a question of where on the theta (angular position) the total function needs to be evaluated-not only a question of evaluating on a discrete radial position. However, if the fixed set of sampling points that have been proposed for the discrete 2D transform are used, where the indices on radial and angular position match, then the results are as expected.

Appendix A.7. Proofs of Orthogonality of the Proposed Kernel
In what follows, we assume the ranges of the variables are such that 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 0 ≤ p, q ≤ N 2 − 1.
Appendix A.7.1. Proof of Orthogonality of the Kernel over the Frequency Indices We state and prove that the following relationship is true: The proof is as follows. We start by substituting the definition of the kernel into the expression: Summing over the index q and using the orthogonality of the discrete complex exponential (Appendix A.3) returns a N 2 δ nn so that n = n and Equation (A26) becomes: where the orthogonality relationship of the discrete complex exponential has been used again.
Appendix A.7.2. Proof of Orthogonality of the Kernel over the Spatial Indices It can be similarly shown that the following orthogonality relationship is also true: The proof is as follows. We start by substituting the definition of the kernel into the expression: Now summation over k gives the right hand side of (A32) and using the discrete orthogonality of the Bessel functions (Appendix A.4) gives for the right hand side: Then, finally summation over n and using the orthogonality of the discrete complex exponential (Appendix A.3) finally gives: qm;pk E + q m ;pk = δ mm δ qq (A34) as required.