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

View article
PeerJ Computer Science

Introduction

The Fourier Transform (FT) is a powerful analytical tool and has proved to be invaluable in many disciplines such as physics, mathematics and engineering. The development of the Fast Fourier Transform (FFT) algorithm (Cooley & Tukey, 1965), which computes the Discrete Fourier Transform (DFT) with a fast algorithm, firmly established the FT as a practical tool in diverse areas, most notably signal and image processing.

In two dimensions, the FFT can still be used to compute the DFT in Cartesian coordinates. However, in many applications such as photoacoustics (Xu, Feng & Wang, 2002) and tomography (Scott et al., 2012; Fahimian et al., 2013; Lee et al., 2008; Lewitt & Matej, 2003), it is often necessary to compute the FT in polar coordinates. Moreover, for functions that are naturally described in polar coordinates, a discrete version of the 2D FT in polar coordinates is needed. There have been some attempts to calculate the FT in polar coordinates, most notably through the Hankel transform, since the zeroth order Hankel transform is known to be a 2D FT in polar coordinates for rotationally symmetric functions. However, prior work has focused on numerically approximating the continuous transform. This stands in contrast to the FT, where the DFT can stand alone as an orthogonal transform, independent of the existence of its continuous counterpart.

The idea of a Polar FT has been previously investigated, where the spatial function is in Cartesian coordinates but its FT is computed in polar coordinates (Averbuch et al., 2006; Abbas, Sun & Foroosh, 2017; Fenn, Kunis & Potts, 2007). FTs have been proposed for non-equispaced data, referred to as Unequally Spaced FFT (USFFT) or Non-Uniform FFT (NUFFT) (Dutt & Rokhlin, 1993; Fourmont, 2003; Dutt & Rokhlin, 1995; Potts, Steidl & Tasche, 2001; Fessler & Sutton, 2003). A recent book gives a unified treatment of these topics (Plonka et al., 2018). Previous work has also considered the implications of using a polar grid (Stark, 1979; Stark & Wengrovitz, 1983). Although the above references demonstrate that the computation of a discrete 2D FT on a polar grid has previously been considered in the literature, there is to date no discrete 2D FT in polar coordinates that exists as a transform in its own right, with its own set of rules of the actual manipulated quantities.

In part I of this two part series, we proposed an independent discrete 2D FT in polar coordinates, which has been defined to be discrete from first principles (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, 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 fpkFqm is given by Fqm=F(fpk)=k=1N11p=MMfpkEqm;pk

where p,k,q,m,n, N1, and N2 are integers such that MnM, where 2M+1=N2 1m,k,N11 and Mp,qM. 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=F1(Fqm)=m=1N11q=MMFqmEqm;pk+

In Eqs. (1) and (2), Eqm;pk± are the kernels of the transformation. These can be chosen as the “non-symmetric” form given by Eqm;pk=1N2n=MMJn(jnkjnmjnN1)jnN12Jn+12(jnk)2inei2πnpN2e+i2πnqN2Eqm;pk+=1N2n=MMJn(jnmjnkjnN1)Jn+12(jnm)2i+ne+i2πnpN2ei2πnqN2

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×(N11) 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,θ) which is space limited to r[0,R]. This implies that the function is zero outside of the circle bounded by r[0,R]. An approximate relationship between sampled values of the continuous function and sampled values of its continuous forward 2D transform F(ρ,ψ) has been derived in the first part of the two-part paper as F(jqmR,2πqN2)2πR2k=1N11p=MMf(jpkRjpN1,2πpN2)1N2n=MM2inJn(jnkjnmjnN1)jnN12Jn+12(jnk)ei2πnpN2e+i2πnqN2

Similarly, an approximate relationship between sampled values of the continuous forward transform F(ρ,ψ) and sampled values of the continuous original function f(r,θ) was shown to be given by f(jpkRjpN1,2πpN2)12πR2m=1N11q=MMF(jqmR,2πqN2)1N2n=MM2inJn(jnmjnkjnN1)Jn+12(jnm)e+i2πnpN2ei2πnqN2

In Eqs. (4) and (5), f(r, θ) is the original function in 2D space and F(ρ,ψ) 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(ρ,ψ), 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(jpkRjpN1,2πpN2)

Then, Fqm is calculated from the 2D DFT scaled by 2πR2, Eq. (1), that is Fqm=2πR2F(fpk)=2πR2k=1N11p=MMfpkEqm;pk

The factor of 2πR2 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(jqmR,2πqN2)—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(ρ,ψ), evaluate it at the sampling points and assign this value to Fqm via Fqm=F(jqmR,2πqN2)

Now, fpk is calculated from a scaled version of the inverse 2D DFT, Eq. (2) that is fpk=12πR2F1(Fqm)=12πR2m=1N11q=MMFqmEqm;pk+

To evaluate if the proposed transform can approximate the continuous transform, the question becomes how well fpk calculated from Eq. (9) approximates f(jpkRjpN1,2πpN2)—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(ρ,ψ) with an effective band limit ρ[0,Wρ]. That is, we suppose that the 2D FT F(ρ,ψ) of f(r,θ) is band-limited, meaning that F(ρ,ψ) is zero for ρWρ=2πW. The variable Wρ 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 2π ensures that the final units are in s1 or m1. The approximate relationship between sampled values of the continuous 2D FT F(ρ,ψ) and sampled values of the original continuous function f(r,θ) was derived in the first part of the paper and is given by F(jqmWρjqN1,2πqN2)2πWρ2k=1N11p=MMf(jpkWρ,2πpN2)1N2n=MM2inJn(jnmjnkjnN1)Jn+12(jnk)ei2πnpN2e+i2πnqN2

Similarly, the inverse relationship between sampled values of F(ρ,ψ) and sampled values of f(r,θ) was shown to be given by f(jpkWρ,2πpN2)Wρ22πm=1N11q=MMF(jqmWρjqN1,2πqN2)1N2n=MM2inJn(jnkjnmjnN1)jnN12Jn+12(jnm)ei2πnqN2e+i2πnpN2

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,θ), evaluate it at the (bandlimited specific) sampling points and assign this value to fpk via fpk=f(jpkWρ,2πpN2)

Then, Fqm is calculated from the discrete transform scaled by 2πWρ2, Eq. (1), that is Fqm=2πWρ2F(fpk)=2πWρ2k=1N11p=MMfpkEqm;pk

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(jqmWρjqN1,2πqN2), 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 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,θ) space where 0rR and 0θ2π. The discrete sampling spaces used for radial and angular sampling points in regular r space (r,θ) and ω frequency (ρ,ψ) space are defined as rpk=jpkRjpN1θp=p2πN2

and ρqm=jqmRψq=q2πN2For a band limited function, the function is assumed band-limited to 0ρWρ, 0ψ2π. The sampling space used for radial and angular sampling points in regular ω frequency space (ρ,ψ) and r space (r,θ) for a bandlimited function is defined as rpk=jpkWρθp=p2πN2and ρqm=jqmWρjqN1ψq=q2πN2

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.

Figure 1: Spatial sampling grid for a space-limited function with R = 1, N1 = 16 and N2 = 15.

Figure 2: Frequency space sampling grid for a space-limited function with R = 1, N1 = 16 and N2 = 15.

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, the polar grids are plotted for the case R = 1, N1=96 and N2=95. These are shown in Figs. 3 and 4 respectively.

Figure 3: Spatial sampling grid for a space-limited function with R = 1, N1 = 96 and N2 = 95.

Figure 4: Frequency space sampling grid for a space-limited function with R = 1, N1 = 96 and N2 = 95.

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>2fmaxwhere 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) WρR=jnN1where Wρ 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 MpM, the order of the Bessel zero ranges from M to M, the required relationship becomes min(jpN1)WρR

The relationships jnN=jnN and j0N1<j±1N1<j±2N1<...<j±MN1 are valid (Lozier, 2003), hence Eq. (20) can be written as j0N1WρR

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 π and that the spacing becomes exactly π in the limit as z. 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)2πzcos[z(n+12)π2]

Therefore, an approximation to the Bessel zero, jnk is given by jnk(2k+n12)π2

Hence, Eq. (21) can be written to choose N1 approximately as N1πWρR=2πWRN12WRwhere 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. 14 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.

Space-limited function

In this section, it is assumed that the function is a space limited function, defined in r[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=πR2πr¯2πR2100where 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 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=12(j01j0N1R+jM1jMN1R)

Hence, Eq. (25) for grid coverage can be written as Ar=[114(j01j0N1+jM1jMN1)2]100

Table 1 shows the different values of grid coverage Ar as the values of N1 and N2 are changed.

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%
DOI: 10.7717/peerj-cs.257/table-1

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 Aρ=πWρ2πρ¯2πWρ2100where ρ¯ denotes the average radius of the gap. Since ρ¯=(ρ01+ρM1)2=(j01+jM1)2R

Then, it follows that Eq. (28) for frequency domain grid coverage can be written as Aρ=[1(j01+jM1)24R2Wρ2]100%

From Eq. (30), it can be observed that the sample grid coverage in the frequency domain is affected by R, Wρ and M. Since N2=2M+1, in order to get a better grid coverage with a fixed Wρ, R and N2 can be adjusted. Table 2 shows the grid coverage Aρ 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 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 Wρ requires an increase in N1 at the same time.

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%
DOI: 10.7717/peerj-cs.257/table-2

Band-limited function

In this section, we suppose that the function is an effectively band limited function, defined on ρ[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=12(j01Wρ+jM1Wρ)

Therefore, the coverage of the grid in space domain can be written as Ar=[1(j01+jM1)24Wρ2R2]100

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 ρ¯=(ρ01+ρM1)2=12(j01j0N1Wρ+jM1jMN1Wρ)

The coverage of the grid in frequency domain can be written as Aρ=[114(j01j0N1+jM1jMN1)2]100

It can be observed that the grid coverage in the frequency domain of a band-limited function is the same as the grid coverage in the space domain of a space limited function.

Conclusion

Based on the discussion above, the following conclusions can be made:

  1. Increasing 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 Wρ 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 fpkFqm, 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=1N11 and n,p,q=MM, where 2M+1=N2. These steps can be summarized succinctly by rewriting Eq. (1) as Fqm=1N2n=MM[2πR2injnN1k=1N11Ym,knN1p=MMfpkein2πpN21DDFTcolumnwise]scaled1DDHTrowwisee+in2πqN2inverse1DDFTcolumnwisewhere the DHT is defined in Baddour & Chouinard (2015) via the transformation matrix Ym,knN1=2jnN1Jn+12(jnk)Jn(jnmjnkjnN1)1m,kN11

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=MM and k=1N11, Eq. (2) this can be expressed as fpk=1N2n=MM[jnN1i+n2πR2m=1N11Yk,mnN1q=MMFqmei2πnqN21DDFT(columnwise)]scaled1DDFT(rowwise)e+i2πnpN2inverse1DDFT(columnwise)

This parallels the steps taken for the continuous case, with each continuous operation (Fourier series, Hankel transform) replaced by its discrete counterpart (DFT, DHT).

Therefore, for both forward and inverse 2D-DFT, the sequence of operations is a DFT of each column of the starting matrix, followed by a DHT of each row, a term-by-term scaling, followed by an IDFT of each column. This is a significant computational improvement because by interpreting the transform this way, the Fast Fourier Transform (FFT) can be used, which reduces the computational time quite significantly in comparison with a direct implementation of the summation definitions in Eqs. (1) and (2).

Interpretation of the sampled forward transform in Matlab terms

To use the built-in Matlab function fft, a few operations are required. First, we define Matlab-friendly indices p=p+(M+1) and n=n+(M+1) so that p,n=MM become p,n=12M+1=1N2 (since 2M+1=N2). That is, the primed variables range from 12M+1 rather than MM. Hence, if the matrix f with entries fpk is defined, where p=1N2,k=1N11, then the first step in which is a column-wise DFT can be written as the Matlab-defined DFT as f¯nk=p=1N2fpke2πi(p1M)(n1M)N2

The overbar denotes a DFT. The definition of DFT in Matlab is actually given by the relationship f¯nk=p=1N2fpke2πi(p1)(n1)N2Since the relationship p=1N2fpke2πi(p1)(n1M)N2=p=1N2fpke2πi(p1M)(n1M)N2 is valid, we can sample the original function to obtain the discrete fpk values, put them in the matrix fpk then shift the matrix fpk 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 f¯nk=fft(circshift(fpk,M+1,1),N2,1)

In matrix operations, this is equivalent to stating that each column of fpk is DFT’ed to yield f¯nk.

The second step in Eq. (35) is a DHT of order n, transforming f¯nkf¯^nl so that the k subscript is Hankel transformed to the l subscript. The overhat denotes a DHT. In order to relate the order n to the index n, we need to shift f¯nk by (M+1) along column direction so that the order ranges from –M to M.

f¯^nl=k=1N112Jn(jnkjnljnN1)jnN1Jn+12(jnk)circshift(f¯nk,(M+1),1){forn=1N2,l=1N11wheren=n(M+1)

By using the Hankel transform matrix defined in Baddour & Chouinard (2015), Eq. (41) can be rewritten as f¯^nl=circshift(f¯nk,(M+1),1)(Yl,knN1)T{forn=1N2,l=1N11wheren=nM1

In matrix operations, this states that each row of f¯nk is DHT’ed to yield f¯^nl. These are now scaled to give the Fourier coefficients of the 2D DFT f¯^nlF¯nl. 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 F¯nl=circshift(2πR2jnN1inf¯^nl,M+1,1){forn=1N2,l=1N11wheren=n(M+1)

This last step is a 1D IDFT for each column of F¯nl to obtain Fql. Using 2M+1=N2, and q=q+1+M, this can be written as Fql=1N2n=1N2F¯nle+i(nM1)2π(q1M)N2forq=1N2,l=1..N11=1N2n=1N2F¯nle+i(n1)2π(q1M)N2=circshift(ifft(F¯nl,N2,1),(M+1),1)

Interpretation of the sampled inverse transform in Matlab terms

Similar to the forward transform, matlab-friendly indices q=q+(M+1) and n=n+(M+1) are also defined. Hence, if the matrix F with entries Fql is defined, where q=1..N2,l=1..N11, it then follows that the first 1D DFT step in Eq. (37) can be written as the Matlab-defined DFT as F¯nl=q=1N2Fqlei(nM1)2π(q1M)N2forn=1N2,l=1N11=q=1N2Fqlei(nM1)2π(q1)N2

If the original function can be sampled as Fql and then put into matrix Fql, then we need an circshift operation. So Eq. (45) can be written as F¯nl=fft(circshift(Fql,M+1,1),N2,1)

Subsequently, a DHT of order n is required, transforming F¯nlF¯^nl so that the l subscript is Hankel transformed to the k subscript. To achieve this, circshift is also needed here.

F¯^nk=circshift(F¯nl,(M+1),1)(Yk,lnN1)T{forn=1N2,l=1N11wheren=nM1

This is followed by a scaling operation to obtain F¯^nkf¯nk and then a circshift by (M+1) so that f¯nk=circshift(jnN12πR2i+nF¯^nk,(M+1),1){forn=1N2,k=1N11wheren=n(M+1)

This last step is a 1D IDFT for each column of f¯nk to get fpk. Using 2M+1=N2, and p=p+1, Eq. (37) can be written as fpk=1N2n=1N2f¯nke+i(nM1)2π(p1M)N2forp=1N2,k=1N11=1N2n=1N2f¯nke+i2π(n1)(p1M)N2=circshift(ifft(f¯nk,N2,1),(M+1),1)

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.

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)=20log10[|C(v)D(v)|max|D(v)|]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 ε=1(N11)N2n=1(N11)N2|ff|where f is the original function and f is the value obtained after sequentially performing a forward and then inverse transform. An ideal precision would result in the absolute error being zero.

Test functions

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

Gaussian

The first function chosen for evaluation is a circular symmetric function which is Gaussian in the radial direction. Specifically, the function in the space domain is given by f(r,θ)=ea2r2where 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(ρ,ψ)=πa2eρ24a2

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.

Figure 5: (A) Original function (Gaussian) and (B) its continuous 2D-DFT (which is also a Gaussian).

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 Wρ=10. From Eq. (21), j0N1RWρ=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 function in both the space and frequency domains is located in the center near the origin. Based on the discussion in “Conclusion”, relatively large values of R and Wρ are needed. The effective space limit R=40 and effective band-limit Wp=30 are thus chosen, which gives j0N1RWρ=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.

Figure 6: (A) sampled continuous transform and (B) discrete forward transform for a Gaussian function with R = 5 and N1 = 17.
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.

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.9115dB and this occurs at the center. The average error is Eavg.=30.4446dB.

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.3842dB and this occurs at the center. The average value of the error is Eavg.=63.8031dB. Clearly, the test with R=40, N1=383 gives a better approximation, which verifies the discussion in “Conclusion”.

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.

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.

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
DOI: 10.7717/peerj-cs.257/table-3
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.
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.

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 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. 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.

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.
Figure 12: Error between the sampled continuous inverse transform and discrete inverse transform for the Gaussian function for R = 5 and N1 = 17.

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.1954dB and this occurs at the center. The average of the error is Eavg.=25.7799dB.

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.2602dB and this occurs at the center. The average of the error is Eavg.=98.0316dB. 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.
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
DOI: 10.7717/peerj-cs.257/table-4
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.
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.

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 Wρ, 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.

Performing sequential 2D-DFT and 2D-IDFT results in ε=4.1656×e17 where ε 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,θ)=sin(ar)ar[3sin(θ)+sin(3θ)+4cos(10θ)+12sin(15θ)]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.

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.

The continuous 2D-FT can be calculated from Baddour (2011) F(ρ,ψ)=n=2πineinψ0fn(r)Jn(ρr)rdrwhere fn(r) is the Fourier series of f(r,θ) and can be written as fn(r)=12πππf(r,θ)einθdθ

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(ρ,ψ)={8πcos(10ψ)ρ10aa2ρ2(a+a2ρ2)10,ρ<a6πisin(ψ)aρρ2+a2+2πisin(3arcsin(aρ))sin(3ψ)ρ2+a28πsin(10arcsin(aρ))cos(10ψ)ρ2+a2+24πisin(15arcsin(aρ))sin(15ψ)ρ2+a2,ρ>a

In the angular direction, the highest frequency term in the function in the space domain is 12sin(15θ). 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,θ) is space-limited at R=15 and band-limited at Wρ=30. However, since most of the energy in the space domain is located at the origin, a relatively large band limit should be chosen based on the discussion in “Conclusion”. Therefore, Wρ=90, N1=430 are chosen.

Forward transform

The error results for the forward 2D-DFT of Four-term sinusoid & Sinc function with Wρ=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.6535dB and this occurs where the discontinuities are located. The average of the error is Erroraverage=38.7831dB.

Figure 17: Error results for the forward 2D Fourier transform of the Four-term sinusoid & Sinc function for Wp = 90 and N1 = 430.

With Wρ=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.

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
DOI: 10.7717/peerj-cs.257/table-5
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.
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.

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”). 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 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 Wρ=90, N1=430 are shown in Fig. 20. The maximum value of the error is Errormax=8.6734dB. The average of the error is Erroraverage=37.8119dB. With Wρ=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.
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
DOI: 10.7717/peerj-cs.257/table-6
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.
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.

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 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 ε=1.3117×e12 where ε is calculated by Eq. (51).

Four-term sinusoid and modified exponential

For the next test function, the function is given by f(r,θ)=earr[3sin(θ)+sin(3θ)+4cos(10θ)+12sin(15θ)]

Its continuous 2D-FT can be calculated as F(ρ,ψ)=6πisin(ψ)ρ2+a2aρρ2+a2+2πisin(3ψ)(ρ2+a2a)3ρ3ρ2+a28πcos(10ψ)(ρ2+a2a)10ρ10ρ2+a2+24πisin(15ψ)(ρ2+a2a)15ρ15ρ2+a2

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.

From graph of the original function and its 2D-DFT, it can be assumed that f(r,θ) is effectively space-limited with R=20, and F(ρ,ψ) is effectively band-limited with Wρ=15, which gives j0N1300. This results in N1=96. However, since the function explodes at the center area in both domains, relatively large values of R and Wρ should give a better approximation. Therefore, another case with R=40, Wρ=30 is tested. In this case, N1=383 is chosen.

In the angular direction, the highest frequency term is 12sin(15θ). 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,Wρ=30,N1=383 are shown in Fig. 24. The maximum value of the error is Errormax=10.1535dB and this occurs at the center area. The average of the error is Erroraverage=32.7619dB. This demonstrates that the discrete function approximates the sampled values of the continuous function quite well.

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.
Inverse transform

The error results with R=40,Wρ=30,N1=383 are shown in Fig. 25.

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.

The maximum value of the error is Errormax=0.5579dB and this occurs at the center. The average of the error isErroraverage=68.7317dB.

Performing 2D-DFT and 2D-IDFT results in ε=1.421×e12, where ε is calculated by Eq. (51).

It can be observed that even for functions with the worst properties, the proposed transform can still be used to approximate the continuous FT with fairly small errors, as long as the function is sampled properly.

Summary and Conclusion

Accuracy and precision of the transform

The proposed discrete 2D-FT in polar coordinates demonstrates an acceptable accuracy in providing discrete estimates to the continuous FT in polar coordinates. In Baddour & Chouinard (2015), Guizar-Sicairos & Gutiérrez-Vega (2004) and Higgins & Munson (1987), the one dimensional Hankel transform of a sinc function showed similar dynamic error, which could be used as a comparative measure. Since the DHT is one step of the proposed discrete 2D-FT, and the definition of the Hankel transform is based on Abbas, Sun & Foroosh (2017), a similar dynamic error should be expected.

The test cases showed that the transform introduced very small errors (ε=1.4004×e12 for worst case) by performing a forward transform and an inverse 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.

  3. Evaluate the transform as a DFT, DHT and IDFT in sequence with pre-calculating the Bessel zeros.

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
DOI: 10.7717/peerj-cs.257/table-7

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.

Supplemental Information

Sample Matlab Code for the discrete 2D Fourier transform in polar coordinates.

DOI: 10.7717/peerj-cs.257/supp-1
1 Citation   Views   Downloads