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

Space-limited functions

Consider a function in the space domain $f(r,\mathrm{\theta })$ which is space limited to $r\in \left[0,R\right]$. This implies that the function is zero outside of the circle bounded by $r\in \left[0,R\right]$. An approximate relationship between sampled values of the continuous function and sampled values of its continuous forward 2D transform $F\left(\mathrm{\rho },\mathrm{\psi }\right)$ has been derived in the first part of the two-part paper as (4) $$F\left(\frac{j_{qm}}{R},\frac{2\mathrm{\pi }q}{N_2}\right)\approx 2\mathrm{\pi }{R}^2\sum _{k=1}^{N_1-1}\sum _{p=-M}^{M}f\left(\frac{j_{pk}R}{j_{p{N}_1}},\frac{2\mathrm{\pi }p}{N_2}\right)\frac{1}{N_2}\sum _{n=-M}^M\frac{2i^{-n}{J}_n\left(\frac{j_{nk}{j}_{nm}}{j_{n{N}_1}}\right)}{j_{n{N}_1}^2{J}_{n+1}^2\left(j_{nk}\right)}{e}^{-i\frac{2\mathrm{\pi }np}{N_2}}{e}^{+i\frac{2\mathrm{\pi }nq}{N_2}}$$

Similarly, an approximate relationship between sampled values of the continuous forward transform $F\left(\mathrm{\rho },\mathrm{\psi }\right)$ and sampled values of the continuous original function $f(r,\mathrm{\theta })$ was shown to be given by (5) $$f\left(\frac{j_{pk}R}{j_{p{N}_1}},\frac{2\mathrm{\pi }p}{N_2}\right)\approx \frac{1}{2\mathrm{\pi }{R}^2}\sum _{m=1}^{N_1-1}\sum _{q=-M}^{M}F\left(\frac{j_{qm}}{R},\frac{2\mathrm{\pi }q}{N_2}\right)\frac{1}{N_2}\sum _{n=-M}^M\frac{2i^n{J}_n\left(\frac{j_{nm}{j}_{nk}}{j_{n{N}_1}}\right)}{J_{n+1}^2\left(j_{nm}\right)}{e}^{+i\frac{2\mathrm{\pi }np}{N_2}}{e}^{-i\frac{2\mathrm{\pi }nq}{N_2}}$$

In Eqs. (4) and (5), f(r, θ) is the original function in 2D space and $F(\mathrm{\rho },\mathrm{\psi })$ 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(\mathrm{\rho },\mathrm{\psi })$, the process is as follows. For the forward transform, we start with the continuous f(r, θ), evaluate it at the sampling points and then assign this value to f_pk via (6) $$f_{pk}=f\left(\frac{j_{pk}R}{j_{pN1}},\frac{2\mathrm{\pi }p}{N_2}\right)$$

Then, $F_{qm}$ is calculated from the 2D DFT scaled by $2\mathrm{\pi }{R}^2$, Eq. (1), that is (7) $F_{qm}=2\mathrm{\pi }{R}^2\mathbb{F}\left(f_{pk}\right)=2\mathrm{\pi }{R}^2\sum _{k=1}^{N_1-1}\sum _{p=-M}^M{f}_{pk}{E}_{qm;pk}^{-}$

The factor of $2\mathrm{\pi }{R}^2$ is necessary so that the evaluation in Eq. (7) matches the expression in Eq. (4). To evaluate if the proposed 2D DFT can be used to approximate the continuous transform, the question becomes how well $F_{qm}$ calculated from the 2D DFT in Eq. (7) approximates $F\left(\frac{j_{qm}}{R},\frac{2\mathrm{\pi }q}{N_2}\right)$—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(\mathrm{\rho },\mathrm{\psi })$, evaluate it at the sampling points and assign this value to $F_{qm}$ via (8) $$F_{qm}=F\left(\frac{j_{qm}}{R},\frac{2\mathrm{\pi }q}{N_2}\right)$$

Now, $f_{pk}$ is calculated from a scaled version of the inverse 2D DFT, Eq. (2) that is (9) $$f_{pk}=\frac{1}{2\mathrm{\pi }{R}^2}{\mathbb{F}}^{-1}\left(F_{qm}\right)=\frac{1}{2\mathrm{\pi }{R}^2}\sum _{m=1}^{N_1-1}\sum _{q=-M}^M{F}_{qm}{E}_{qm;pk}^{+}$$

To evaluate if the proposed transform can approximate the continuous transform, the question becomes how well $f_{pk}$ calculated from Eq. (9) approximates $f\left(\frac{j_{pk}R}{j_{pN1}},\frac{2\mathrm{\pi }p}{N_2}\right)$—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\left(\mathrm{\rho },\mathrm{\psi }\right)$ with an effective band limit $\mathrm{\rho }\in \left[0,W_{\mathrm{\rho }}\right]$. That is, we suppose that the 2D FT $F\left(\mathrm{\rho },\mathrm{\psi }\right)$ of $f(r,\mathrm{\theta })$ is band-limited, meaning that $F\left(\mathrm{\rho },\mathrm{\psi }\right)$ is zero for $\mathrm{\rho }\ge W_{\mathrm{\rho }}=2\mathrm{\pi }W$. The variable $W_{\mathrm{\rho }}$ 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\mathrm{\pi }$ ensures that the final units are in $s^{-1}$ or $m^{-1}$. The approximate relationship between sampled values of the continuous 2D FT $F\left(\mathrm{\rho },\mathrm{\psi }\right)$ and sampled values of the original continuous function $f\left(r,\mathrm{\theta }\right)$ was derived in the first part of the paper and is given by (10) $$F\left(\frac{j_{qm}{W}_{\mathrm{\rho }}}{j_{q{N}_1}},\frac{2\mathrm{\pi }q}{N_2}\right)\approx \frac{2\mathrm{\pi }}{W_{\mathrm{\rho }}^2}\sum _{k=1}^{N_1-1}\sum _{p=-M}^{M}f\left(\frac{j_{pk}}{W_{\mathrm{\rho }}},\frac{2\mathrm{\pi }p}{N_2}\right)\frac{1}{N_2}\sum _{n=-M}^M\frac{2i^{-n}{J}_n\left(\frac{j_{nm}{j}_{nk}}{j_{n{N}_1}}\right)}{J_{n+1}^2\left(j_{nk}\right)}{e}^{-i\frac{2\mathrm{\pi }np}{N_2}}{e}^{+i\frac{2\mathrm{\pi }nq}{N_2}}$$

Similarly, the inverse relationship between sampled values of $F\left(\mathrm{\rho },\mathrm{\psi }\right)$ and sampled values of $f(r,\mathrm{\theta })$ was shown to be given by (11) $f\left(\frac{j_{pk}}{W_{\mathrm{\rho }}},\frac{2\mathrm{\pi }p}{N_2}\right)\approx \frac{W_{\mathrm{\rho }}^2}{2\mathrm{\pi }}\sum _{m=1}^{N_1-1}\sum _{q=-M}^{M}F\left(\frac{j_{qm}{W}_{\mathrm{\rho }}}{j_{q{N}_1}},\frac{2\mathrm{\pi }q}{N_2}\right)\frac{1}{N_2}\sum _{n=-M}^M\frac{2i^n{J}_n\left(\frac{j_{nk}{j}_{nm}}{j_{n{N}_1}}\right)}{j_{n{N}_1}^2{J}_{n+1}^2\left(j_{nm}\right)}{e}^{-i\frac{2\mathrm{\pi }nq}{N_2}}{e}^{+i\frac{2\mathrm{\pi }np}{N_2}}$

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\left(r,\mathrm{\theta }\right)$, evaluate it at the (bandlimited specific) sampling points and assign this value to $f_{pk}$ via (12) $$f_{pk}=f\left(\frac{j_{pk}}{W_{\mathrm{\rho }}},\frac{2\mathrm{\pi }p}{N_2}\right)$$

Then, $F_{qm}$ is calculated from the discrete transform scaled by $\frac{2\mathrm{\pi }}{W_{\mathrm{\rho }}^2}$, Eq. (1), that is (13) $$F_{qm}=\frac{2\mathrm{\pi }}{W_{\mathrm{\rho }}^2}\mathbb{F}\left(f_{pk}\right)=\frac{2\mathrm{\pi }}{W_{\mathrm{\rho }}^2}\sum _{k=1}^{N_1-1}\sum _{p=-M}^M{f}_{pk}{E}_{qm;pk}^{-}$$

To evaluate if the proposed 2D DFT can be used to approximate the continuous transform, the question is how well $F_{qm}$ calculated from Eq. (13) approximates $F\left(\frac{j_{qm}{W}_{\mathrm{\rho }}}{j_{q{N}_1}},\frac{2\mathrm{\pi }q}{N_2}\right)$, 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 $f_{pk}$ and $F_{qm}$ calculated via the proposed 2D DFT are in fact reasonable approximations to the sampled values of the continuous functions, as stated in the objectives of the paper.

Sampling points

For a space-limited function, we assume that the original function of interest is defined over continuous $\left(r,\mathrm{\theta }\right)$ space where $0\le r\le R$ and $0\le \mathrm{\theta }\le 2\mathrm{\pi }$. The discrete sampling spaces used for radial and angular sampling points in regular $\overrightarrow{r}$ space $\left(r,\mathrm{\theta }\right)$ and $\overrightarrow{\mathrm{\omega }}$ frequency $\left(\mathrm{\rho },\mathrm{\psi }\right)$ space are defined as (14) $$r_{pk}=\frac{j_{pk}R}{j_{p{N}_1}}{\mathrm{\theta }}_p=\frac{p2\mathrm{\pi }}{N_2}$$

and (15) $${\mathrm{\rho }}_{qm}=\frac{j_{qm}}{R}{\mathrm{\psi }}_q=\frac{q2\mathrm{\pi }}{N_2}$$For a band limited function, the function is assumed band-limited to $0\le \mathrm{\rho }\le W_{\mathrm{\rho }}$, $0\le \mathrm{\psi }\le 2\mathrm{\pi }$. The sampling space used for radial and angular sampling points in regular $\overrightarrow{\mathrm{\omega }}$ frequency space $\left(\mathrm{\rho },\mathrm{\psi }\right)$ and $\overrightarrow{r}$ space $\left(r,\mathrm{\theta }\right)$ for a bandlimited function is defined as (16) $$r_{pk}=\frac{j_{pk}}{W_{\mathrm{\rho }}}{\mathrm{\theta }}_p=\frac{p2\mathrm{\pi }}{N_2}$$and (17) $${\mathrm{\rho }}_{qm}=\frac{j_{qm}{W}_{\mathrm{\rho }}}{j_{q{N}_1}}{\mathrm{\psi }}_q=\frac{q2\mathrm{\pi }}{N_2}$$

Clearly, the density of the sampling points depends on the numbers of points chosen, that is on $N_1$ and $N_2$. Also clear is the fact that the grid is not equispaced in the radial variable. The sampling grid for a space-limited function are plotted below to enable visualization. In the first instance, the polar grids are plotted for the case $R=1$, $N_1=16$ and $N_2=15$. These are shown in space (r space) and frequency (ρ space) in Figs. 1 and 2 respectively. It should be noted that although we refer the grids in this article as polar grids, they are not true polar grids in the sense of equispaced sampling in the radial and angular coordinates.

Clearly, the grids in Figs. 1 and 2 are fairly sparse, but the low values of $N_2$ and $N_1$ have been chosen so that the structure of the sampling points can be easily seen. It can be observed that there is a hole at the center area in both domains which is caused by the special sampling points. For higher values of the $N_2$ and $N_1$, the grid becomes fairly dense, obtaining good coverage of both spaces, but details are harder to observe. To demonstrate, the polar grids are plotted for the case R = 1, $N_1=96$ and $N_2=95$. These are shown in Figs. 3 and 4 respectively.

From Figs. 3 and 4, by choosing higher values of $N_1$ and $N_2$, the sampling grid becomes denser, however there is still a gap in the center area. The sampling grids for band-limited functions are not plotted here since the sample grid for a band-limited function has the same shape as with space limited function but the domains are reversed.

Sample grid analysis

From part I of the paper, it was shown that the 2D-FT can be interpreted as a DFT in the angular direction, a DHT in the radial direction and then an IDFT in the angular direction. Hence, the sample size in the angular direction could have been decided by the Nyquist sampling theorem (Shannon, 1984), which states that (18) $$f_s>2f_{max}$$where $f_s$ is the sample frequency and $f_{max}$ is the highest frequency or band limit.

In the radial direction, the necessary relationship for the DHT is given by Baddour & Chouinard (2015) (19) $$W_{\mathrm{\rho }}R=j_{n{N}_1}$$where $W_{\mathrm{\rho }}$ is the effective band-limit, $R$ is the effective space limit and $j_{nN}$ is the Nth zero of $J_n\left(r\right)$. For the 2D FT, since $-M\le p\le M$, the order of the Bessel zero ranges from $-M$ to $M$, the required relationship becomes (20) $$\mathrm{m}\mathrm{i}\mathrm{n}(j_{p{N}_1})\ge W_{\mathrm{\rho }}R$$

The relationships $j_{nN}=j_{-nN}$ and $j_{0N_1}<j_{\pm 1N_1}<j_{\pm 2N_1}<...<j_{\pm M{N}_1}$ are valid (Lozier, 2003), hence Eq. (20) can be written as (21) $$j_{0N_1}\ge W_{\mathrm{\rho }}R$$

It is pointed out in Baddour (2019) and Guizar-Sicairos & Gutiérrez-Vega (2004) that the zeros of $J_n\left(z\right)$ are almost evenly spaced at intervals of $\mathrm{\pi }$ and that the spacing becomes exactly $\mathrm{\pi }$ in the limit as $z\to \mathrm{\infty }$. 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 (22) $$J_n\left(z\right)\approx \sqrt{\frac{2}{\mathrm{\pi }z}}\cos \left[z-\left(n+\frac{1}{2}\right)\frac{\mathrm{\pi }}{2}\right]$$

Therefore, an approximation to the Bessel zero, $j_{nk}$ is given by (23) $$j_{nk}\approx \left(2k+n-\frac{1}{2}\right)\frac{\mathrm{\pi }}{2}$$

Hence, Eq. (21) can be written to choose $N_1$ approximately as (24) $$\begin{array}{rl}{N}_1\mathrm{\pi }\ge & W_{\mathrm{\rho }}R=2\mathrm{\pi }WR\\ & \Rightarrow N_1\ge 2WR\end{array}$$where the reader is reminded that the units of $W$ is m⁻¹ (the space equivalent of Hz). $N_1/R$ is the spatial sampling frequency and we see that Eq. (24) effectively makes the same statement as Eq. (18), as it should.

Intuitively, more sample points lead to more information captured, which gives an expectation that increasing $N_1$ or $N_2$ individually will give a better sampling grid coverage. However, it can be seen from Figs. 1–4 that there is a gap in the center of the sample grid. From Eqs. (14) and (15), the area of the gap in the center is related to the ranges of $p$ and $k$, that is $N_2$ and $N_1$. In the sections below, it is assumed that the sampling theorems are already satisfied (that is, an appropriate space and band limit is selected) and the relationship between $N_2$, $N_1$ and the size of the gap will be discussed.

Conclusion

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

1. Increasing $N_2$ (angular direction) tends to decrease the sampling grid coverage in both domains. Increasing $N_1$ (radial direction) tends to increase the sampling coverage in the space domain for a space-limited function and in the frequency domain for a frequency-limited function. So, if a signal changes sharply in the angular direction such that large values of $N_2$ are needed, a large value of $N_1$ is also needed to compensate for the effect of increasing $N_2$ on the grid coverage.

2. For a space-limited function, if there is a lot of energy at the origin in the space domain, a larger value of $N_1$ will be required to ensure that the sampling grid gets as close to the origin as possible in the space domain. If the function has a lot of energy at the origin in the frequency domain, a large value for both $N_1$ and $R$ will be required to ensure adequate grid coverage.

3. For a band-limited function, if there is a lot of energy at the origin in the frequency domain, a large value of $N_1$ will be needed to ensure that the sample grid gets as close to the origin as possible in the frequency domain. If the function has a lot of energy at the origin in the space domain, large values for both $N_1$ and $W_{\mathrm{\rho }}$ are required.

Forward transform

The values $f_{pk}$ can be considered as the entries in a matrix. To transform $f_{pk}\to F_{qm}$, the operation is performed as a sequence of steps which are a 1D DFT (column-wise), followed by a scaled 1D DHT (row-wise), finally followed by a 1D IDFT (column-wise). The reader is reminded that the range of indices is given by $m,k=1\dots N_1-1$ and $n,p,q=-M\dots M$, where $2M+1=N_2$. These steps can be summarized succinctly by rewriting Eq. (1) as (35) $$F_{qm}=\underset{\mathrm{i}\mathrm{n}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{s}\mathrm{e}1\mathrm{D}\mathrm{D}\mathrm{F}\mathrm{T}\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{m}\mathrm{n}-\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}}{\underbrace{\frac{1}{N_2}\sum _{n=-M}^M\underset{\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{e}\mathrm{d}1\mathrm{D}\mathrm{D}\mathrm{H}\mathrm{T}\mathrm{r}\mathrm{o}\mathrm{w}-\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}}{\underbrace{\left[\frac{2\mathrm{\pi }{R}^2{i}^{-n}}{j_{n{N}_1}}\sum _{k=1}^{N_1-1}{Y}_{m,k}^{n{N}_1}\underset{1\mathrm{D}\mathrm{D}\mathrm{F}\mathrm{T}\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{m}\mathrm{n}-\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}}{\underbrace{\sum _{p=-M}^M{f}_{pk}{e}^{-in\frac{2\mathrm{\pi }p}{N_2}}}}\right]}}{e}^{+in\frac{2\mathrm{\pi }q}{N_2}}}}$$where the DHT is defined in Baddour & Chouinard (2015) via the transformation matrix (36) $$Y_{m,k}^{n{N}_1}=\frac{2}{j_{n{N}_1}{J}_{n+1}^2\left(j_{nk}\right)}{J}_n\left(\frac{j_{nm}{j}_{nk}}{j_{n{N}_1}}\right)1\le m,k\le N_1-1$$

Matlab code for the DHT is described in Baddour & Chouinard (2017). The inverse 2D DFT can be similarly interpreted, as shown in “Inverse Transform”.

Inverse transform

The steps of the inverse 2D DFT are the reverse of the steps outlined above for the forward 2D DFT. For $p=-M\dots M$ and $k=1\dots N_1-1$, Eq. (2) this can be expressed as (37) $$f_{pk}=\underset{\mathrm{i}\mathrm{n}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{s}\mathrm{e}1\mathrm{D}\mathrm{D}\mathrm{F}\mathrm{T}(\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{m}\mathrm{n}-\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e})}{\underbrace{\frac{1}{N_2}\sum _{n=-M}^M\underset{\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{e}\mathrm{d}1\mathrm{D}\mathrm{D}\mathrm{F}\mathrm{T}(\mathrm{r}\mathrm{o}\mathrm{w}-\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e})}{\underbrace{\left[\frac{j_{n{N}_1}{i}^{+n}}{2\mathrm{\pi }{R}^2}\sum _{m=1}^{N_1-1}{Y}_{k,m}^{n{N}_1}\underset{1\mathrm{D}\mathrm{D}\mathrm{F}\mathrm{T}(\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{m}\mathrm{n}-\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e})}{\underbrace{\sum _{q=M}^M{F}_{qm}{e}^{-i\frac{2\mathrm{\pi }nq}{N_2}}}}\right]}}{e}^{+i\frac{2\mathrm{\pi }np}{N_2}}}}$$

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^{\prime }=p+(M+1)$ and $n^{\prime }=n+(M+1)$ so that $p,n=-M\dots M$ become $p^{\prime },n^{\prime }=1\dots 2M+1=1\dots N_2$ (since $2M+1=N_2$). That is, the primed variables range from $1\dots 2M+1$ rather than $-M\dots M$. Hence, if the matrix $\mathbf{f}$ with entries $f_{p^{\prime }k}$ is defined, where $p^{\prime }=1\dots N_2,k=1\dots N_1-1$, then the first step in which is a column-wise DFT can be written as the Matlab-defined DFT as (38) $${\overline{f}}_{n^{\prime }k}=\sum _{p^{\prime }=1}^{N_2}{f}_{pk}{e}^{\frac{-2\mathrm{\pi }i(p^{\prime }-1-M)(n^{\prime }-1-M)}{N_2}}$$

The overbar denotes a DFT. The definition of DFT in Matlab is actually given by the relationship (39) $${\overline{f}}_{n^{\prime }k}=\sum _{p^{\prime }=1}^{N_2}{f}_{p^{\prime }k}{e}^{\frac{-2\mathrm{\pi }i(p^{\prime }-1)(n^{\prime }-1)}{N_2}}$$Since the relationship $\sum _{p^{\prime }=1}^{N_2}{f}_{p^{\prime }k}{e}^{\frac{-2\mathrm{\pi }i(p^{\prime }-1)(n{}^{\prime }-1-M)}{N_2}}=\sum _{p^{\prime }=1}^{N_2}{f}_{pk}{e}^{\frac{-2\mathrm{\pi }i(p^{\prime }-1-M)(n^{\prime }-1-M)}{N_2}}$ is valid, we can sample the original function to obtain the discrete $f_{pk}$ values, put them in the matrix $f_{p^{\prime }k}$ then shift the matrix $f_{p^{\prime }k}$ by $M+1$ along the column direction. In Matlab, the function $circshift\left(A,K,\mathrm{d}\mathrm{i}\mathrm{m}\right)$ can be used, which circularly shifts the values in array $A$ by $K$ positions along dimension dim. Inputs $K$ and $\mathrm{d}\mathrm{i}\mathrm{m}$ must be scalars. Specifically, $\mathrm{d}\mathrm{i}\mathrm{m}=1$ indicates the columns of matrix $A$ and $\mathrm{d}\mathrm{i}\mathrm{m}=2$ indicates the rows of matrix A. Hence, Eq. (38) can be written as (40) $${\overline{f}}_{n^{\prime }k}=fft\left(circshift\left(f_{p^{\prime }k},M+1,1\right),N_2,1\right)$$

In matrix operations, this is equivalent to stating that each column of $f_{p^{\prime }k}$ is DFT’ed to yield ${\overline{f}}_{n^{\prime }k}$.

The second step in Eq. (35) is a DHT of order $n$, transforming ${\overline{f}}_{n\mathrm{\prime }k}\to {\hat{\overline{f}}}_{n\mathrm{\prime }l}$ so that the $k$ subscript is Hankel transformed to the $l$ subscript. The overhat denotes a DHT. In order to relate the order $n$ to the index $n^{\prime }$, we need to shift ${\overline{f}}_{n^{\prime }k}$ by $-(M+1)$ along column direction so that the order ranges from –M to M.

By using the Hankel transform matrix defined in Baddour & Chouinard (2015), Eq. (41) can be rewritten as (42) $${\hat{\overline{f}}}_{n{}^{\mathrm{\prime }}l}=circshift({\overline{f}}_{n^{\mathrm{\prime }}k},-(M+1),1){\left(Y_{l,k}^{n{N}_1}\right)}^T\{\begin{array}{c}\mathrm{f}\mathrm{o}\mathrm{r}{n}^{\mathrm{\prime }}=1\dots N_2,l=1\dots N_1-1\\ \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}n=n^{\mathrm{\prime }}-M-1\end{array}$$

In matrix operations, this states that each row of ${\overline{f}}_{n^{\prime }k}$ is DHT’ed to yield ${\hat{\overline{f}}}_{n^{\mathrm{\prime }}l}$. These are now scaled to give the Fourier coefficients of the 2D DFT ${\hat{\overline{f}}}_{n^{\mathrm{\prime }}l}\to {\overline{F}}_{n^{\mathrm{\prime }}l}$. In order to proceed to an IDFT in the next step, it is necessary to shift the matrix by $M+1$ along the column direction after scaling (43) $${\overline{F}}_{n^{\mathrm{\prime }}l}=\mathrm{c}\mathrm{i}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}\left(\frac{2\mathrm{\pi }{R}^2}{j_{n{N}_1}}{i}^{-n}{\hat{\overline{f}}}_{n^{\mathrm{\prime }}l},M+1,1\right)\{\begin{array}{c}\mathrm{f}\mathrm{o}\mathrm{r}{n}^{\mathrm{\prime }}=1\dots N_2,l=1\dots N_1-1\\ \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}n=n^{\mathrm{\prime }}-\left(M+1\right)\end{array}$$

This last step is a 1D IDFT for each column of ${\overline{F}}_{n{}^{\prime }l}$ to obtain $F_{ql}$. Using $2M+1=N_2$, and $q^{\prime }=q+1+M$, this can be written as (44) $$\begin{array}{rl}& F_{q^{\prime }l}=\frac{1}{N_2}\sum _{n^{\prime }=1}^{N_2}{\overline{F}}_{nl}{e}^{+i\left(n^{\prime }-M-1\right)\frac{2\mathrm{\pi }\left(q^{\prime }-1-M\right)}{N_2}}\mathrm{f}\mathrm{o}\mathrm{r}{q}^{\prime }=1\dots N_2,l=\mathrm{1..}{N}_1-1\\ & =\frac{1}{N_2}\sum _{n^{\prime }=1}^{N_2}{\overline{F}}_{n^{\prime }l}{e}^{+i(n^{\prime }-1)\frac{2\mathrm{\pi }\left(q^{\prime }-1-M\right)}{N_2}}\\ & =\mathrm{c}\mathrm{i}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}\left(ifft\left({\overline{F}}_{n^{\prime }l},N_2,1\right),-\left(M+1\right),1\right)\end{array}$$

Interpretation of the sampled inverse transform in Matlab terms

Similar to the forward transform, matlab-friendly indices $q^{\prime }=q+(M+1)$ and $n^{\prime }=n+(M+1)$ are also defined. Hence, if the matrix $F$ with entries $F_{q^{\prime }l}$ is defined, where $q^{\prime }=\mathrm{1..}{N}_{2,}l=\mathrm{1..}{N}_1-1$, it then follows that the first 1D DFT step in Eq. (37) can be written as the Matlab-defined DFT as (45) $$\begin{array}{rl}& {\overline{F}}_{n^{{}^{\prime }}l}=\sum _{q^{\prime }=1}^{N_2}{F}_{ql}{e}^{-i(n^{\prime }-M-1)\frac{2\mathrm{\pi }(q^{\prime }-1-M)}{N_2}}\mathrm{f}\mathrm{o}\mathrm{r}{n}^{\prime }=1\dots N_2,l=1\dots N_1-1\\ & =\sum _{q^{\prime }=1}^{N_2}{F}_{q^{\prime }l}{e}^{-i(n^{\prime }-M-1)\frac{2\mathrm{\pi }(q^{\prime }-1)}{N_2}}\end{array}$$

If the original function can be sampled as $F_{ql}$ and then put into matrix $F_{q^{\prime }l}$, then we need an $circshift$ operation. So Eq. (45) can be written as (46) $${\overline{F}}_{n^{\prime }l}=fft\left(circshift(F_{q^{\prime }l},M+1,1),N_2,1\right)$$

Subsequently, a DHT of order $n$ is required, transforming ${\overline{F}}_{n^{\prime }l}\to {\hat{\overline{F}}}_{n^{\prime }l}$ so that the $l$ subscript is Hankel transformed to the $k$ subscript. To achieve this, $circshift$ is also needed here.

This is followed by a scaling operation to obtain ${\hat{\overline{F}}}_{n^{\prime }k}\to {\overline{f}}_{n^{\prime }k}$ and then a $circshift$ by $(M+1)$ so that (48) $${\overline{f}}_{n^{\prime }k}=circshift\left(\frac{j_{n{N}_1}}{2\mathrm{\pi }{R}^2}{i}^{+n}{\hat{\overline{F}}}_{n^{\prime }k},(M+1),1\right)\{\begin{array}{c}\mathrm{f}\mathrm{o}\mathrm{r}{n}^{\prime }=1\dots N_2,k=1\dots N_1-1\\ \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}n=n^{\prime }-\left(M+1\right)\end{array}$$

This last step is a 1D IDFT for each column of ${\overline{f}}_{n^{\prime }k}$ to get $f_{p^{\prime }k}$. Using $2M+1=N_2$, and $p^{\prime }=p+1$, Eq. (37) can be written as (49) $$\begin{array}{rl}& f_{p^{\prime }k}=\frac{1}{N_2}\sum _{n^{\prime }=1}^{N_2}{\overline{f}}_{nk}{e}^{+i\left(n^{\prime }-M-1\right)\frac{2\mathrm{\pi }\left(p^{\prime }-1-M\right)}{N_2}}\mathrm{f}\mathrm{o}\mathrm{r}{p}^{\prime }=1\dots N_2,k=1\dots N_1-1\\ & =\frac{1}{N_2}\sum _{n^{\prime }=1}^{N_2}{\overline{f}}_{n^{\prime }k}{e}^{+i\frac{2\mathrm{\pi }\left(n^{\prime }-1\right)\left(p^{\prime }-1-M\right)}{N_2}}\\ & =circshift\left(ifft\left({\overline{f}}_{n^{\prime }k},N_2,1\right),-(M+1),1\right)\end{array}$$

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.

Method for testing the algorithm

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 (52) $$f(r,\mathrm{\theta })=e^{-a^2{r}^2}$$where a is some real constant. Since the function is circularly symmetric, the 2D-DFT is a zeroth-order Hankel Transform (Poularikas, 2010) and is given by (53) $$F(\mathrm{\rho },\mathrm{\psi })=\frac{\mathrm{\pi }}{a^2}{e}^{\frac{-{\mathrm{\rho }}^2}{4a^2}}$$

The graphs for the original function and its continuous 2D-DFT (which is also a Gaussian) are plotted with $a=1$ and shown in Fig. 5. From Fig. 5, the function is circular symmetric and fairly smooth in the radial direction. Moreover, the function can be considered as either an effectively space-limited function or an effectively band-limited function. For the purposes of testing it, it shall be considered as a space-limited function and Eqs. (14) and (15) will be used to proceed with the forward and inverse transform in sequence.

To perform the transform, the following variables need to be chosen: $N_2$, $R$ and $N_1$. In the angular direction, since the function in the spatial domain is circularly symmetric, $N_2$ can be chosen to be small. Thus, $N_2=15$ is chosen.

In the radial direction, from plotting the function, it can be seen that the effective space limit can be taken to be $R=5$ and the effective band limit can be taken to be $W_{\mathrm{\rho }}=10$. From Eq. (21), $j_{0N_1}\ge R\cdot W_{\mathrm{\rho }}=50$. Therefore, $N_1=17$ is chosen (we could also have obtained a rough estimate of $N_1$ from Eq. (24)). However, most of the energy of the function in both the space and frequency domains is located in the center near the origin. Based on the discussion in “Conclusion”, relatively large values of $R$ and $W_{\mathrm{\rho }}$ are needed. The effective space limit $R=40$ and effective band-limit $W_p=30$ are thus chosen, which gives $j_{0N_1}\ge R\cdot W_{\mathrm{\rho }}=1200$. Therefore $N_1=383$ is chosen in order to satisfy this constraint. Both cases discussed here ($N_1=17$ and $N_1=383$) are tested in following.

Forward transform

Test results with $R=5$, $N_1=17$ are shown in Figs. 6 and 7. Figure 6 shows the sampled continuous forward transform and the discrete forward transform. Figure 7 shows the error between the sampled values of the continuous transform and the discretely calculated values.

From Fig. 7, it can be observed that the error gets bigger at the center, which is as expected because the sampling grid shows that the sampling points can never attain the origin. The maximum value of the error is $E_{max}=-0.9115\mathrm{d}\mathrm{B}$ and this occurs at the center. The average error is $E_{\mathrm{a}\mathrm{v}\mathrm{g}.}=-30.4446\mathrm{d}\mathrm{B}$.

Error test results with $R=40$, $N_1=383$ are shown in Fig. 8. Similar to the previous case, the error gets larger at the center, as expected. However, the maximum value of the error is $E_{max}=-8.3842\mathrm{d}\mathrm{B}$ and this occurs at the center. The average value of the error is $E_{\mathrm{a}\mathrm{v}\mathrm{g}.}=-63.8031\mathrm{d}\mathrm{B}$. Clearly, the test with $R=40$, $N_1=383$ gives a better approximation, which verifies the discussion in “Conclusion”.

With $R=40$, Table 3 shows the errors (max and average error) with respect to different value of $N_1$ and $N_2$. The trends as functions of $N_1$ and $N_2$ are shown as plots in Figs. 9 and 10.

Table 3:

Error (dB) of forward transform of Gaussian function with R = 40, different value of N₁ and N₂.

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

From Fig. 9, it can be seen that when $N_1$ individually ($N_2$ is fixed at $N_2=15$) is less than the minimum of 383 obtained from the sampling theorem, increasing $N_1$ will lead to smaller errors, as expected. When $N_1$ is bigger than the sampling-theorem threshold of 383, increasing $N_1$ still decreases the error which verifies the discussion about sample grid coverage in “Conclusion”. Increasing $N_1$ tends to increase the sample grid coverage and capture more information at the center area and thus leads to smaller errors.

From Fig. 10, increasing $N_2$ alone (i.e., without a corresponding increase in $N_1$) leads to larger errors, both ${\mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}}_{max}$ and ${\mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}}_{\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{g}\mathrm{e}}$. Although at first counterintuitive, this result is actually reasonable because the function is radially symmetric which implies that $N_2=1$ should be sufficient based on the sampling theorem for the angular direction. Therefore, increasing $N_2$ will not lead to a better approximation. Moreover, from the discussion of the sample grid coverage in “Conclusion”, the sampling grid coverage in both domains gets worse when $N_2$ gets bigger because more information from the center is lost. This problem can be solved by increasing $N_1$ at the same time, but it could be computationally time consuming. Therefore, choosing $N_2$ properly is very important from the standpoint of accuracy and computational efficiency.

Inverse transform

Test results for the inverse transform with $R=5$, $N_1=17$ are shown in Figs. 11 and 12. Figure 11 shows the sampled continuous inverse transform and discrete inverse transform and Fig. 12 shows the error between the sampled continuous and discretely calculated values.

Similar to the case for the forward transform, the error gets larger at the center, which is as expected because the sampling grid shows that the sampling points never attain the center. The maximum value of the error is $E_{max}=3.1954\mathrm{d}\mathrm{B}$ and this occurs at the center. The average of the error is $E_{\mathrm{a}\mathrm{v}\mathrm{g}.}=-25.7799\mathrm{d}\mathrm{B}$.

Error test results for the inverse transform with $R=40$, $N_1=383$ are shown in Fig. 13. In this case, the maximum value of the error is $E_{max}=-12.2602\mathrm{d}\mathrm{B}$ and this occurs at the center. The average of the error is $E_{\mathrm{a}\mathrm{v}\mathrm{g}.}=-98.0316\mathrm{d}\mathrm{B}$. Table 4 shows the errors with respect to different value of $N_1$ and $N_2$, from which Figs. 14 and 15 demonstrate the trend.

Table 4:

Error (dB) of inverse transform of Gaussian function with R = 40, different value of N₁ and N₂.

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

From Fig. 15 it can be observed that increasing $N_1$ tends to improve the result but not significantly. This could be explained by the discussion for $R=40$, $N_1=383$ that with fixed $R$ and $W_{\mathrm{\rho }}$, increasing $N_1$ will not allow the sampling grid in the frequency domain to get any closer to the origin to capture more information. From Fig. 15, increasing $N_2$ (with fixed $N_1=383$) leads to a worse approximation which verifies the discussion for $R=40$, $N_1=383$.

Performing sequential 2D-DFT and 2D-IDFT results in $\mathrm{\epsilon }=4.1656\times e^{-17}$ where $\mathrm{\epsilon }$ is calculated with Eq. (51). Therefore, performing sequential forward and inverse transforms does not add much error.