1.

Introduction

The Wigner distribution function (WDF) is a classical time-frequency analysis tool that has been applied in the areas of signal processing and optics widely. It allows us to interpret various systems geometrically in a space-frequency domain such as Fourier, fractional, Fresnel, and linear canonical transforms,¹ optical propagation problems including paraxial² and non-paraxial propagation,³ digital holography,^4–6 and phase retrieval.⁷ In addition, linking the description of wave optics and the ray concept in geometrical optics by connecting the position on the source plane and the light propagation angle has been demonstrated.⁸ Its cross-terms can also characterize the coherence information of partially coherent optical fields.^9,10 Furthermore, the linear canonical transform (LCT) is a type of linear integral transform and the generalization of the Fourier transform, fractional Fourier transform, and Fresnel transform. In Fourier optics, it is a significant tool used to describe first-order optical systems, and it has also been applied in various areas including quantitative phase reconstruction¹¹ and filter design.¹²

In recent years, many novel distributions combining the WDF and LCT have been proposed. In a recent survey of such distributions, we identified the Bai distribution function (BDF)¹³ as the most independently interesting one among them.¹⁴ As a generalization of the WDF, it has been demonstrated that the BDF can extend the application of the WDF in linear frequency-modulated signal parameter estimation¹³ and partial coherence¹⁴ areas. It may also link the description of geometrical and wave optics not only by connecting the positions and propagation angles but also by connecting the positions on the source plane and the positions on the observation plane at any propagation distances¹⁴ or after arbitrary first-order optical systems. It can interpolate between the instantaneous auto-correlation function (IAF) and the WDF of an optical field, which creates a continuum of representations of the signal that may generalize the function of either extremum in a variety of situations. Furthermore, its cross-terms can capture the coherence property or the interference between two optical fields at different propagation distances during the entire propagation process. Therefore, the BDF can be potentially applied in various areas of optics, not least those that incorporate partial coherence. However, the BDF is only defined for continuous signals currently. This is insufficient for simulation purposes, as the signals must necessarily be discrete. It is also insufficient for experimental data, which is captured using digital cameras. This necessitates the development of a BDF analysis workflow for discrete signals. In this paper, we construct a discrete BDF formulation. We also analyze the sampling requirements for the proposed discrete BDF and examine their accuracy by comparing the analytical continuous BDF of a continuous test signal and its discrete numerical approximation. Our results will provide an efficient and accurate numerical calculation method for a discrete BDF, which can benefit various areas and problems in optics mentioned above.

This paper is organized as follows. We start by introducing the definition of the WDF and BDF in Sec. 2. Then, we review several discrete WDF definitions and discuss their strengths, weaknesses, and significance in Sec. 3. The derivation of a discrete BDF is presented mathematically and geometrically in Sec. 4. While arrived at by different methods, we show that one of the most important definitions of a discrete WDF is a special case of the proposed discrete BDF. In Sec. 4.1, we first derive a mathematical expression for a discrete-space BDF. The sampling requirement is analyzed using the phase-space diagram (PSD) in Sec. 4.2. A PSD analysis tracks the support of the WDF of a 1D signal to track its width and bandwidth as a sequence of operations is carried out on the signal. We consider the BDF to consist of two parts: It is the LCT of the IAF and determines the PSD of the IAF of a signal of a given width and bandwidth. Then, we track that PSD as it is transformed by the direct method algorithm for calculating the discrete-space LCT. In Sec. 4.3, based on the derived discrete-space BDF, a discrete BDF calculation process is presented. In addition, some simulation results of the discrete BDF are provided, and a comparison with its analytical continuous version is carried out in Sec. 5. Finally, we present our conclusions.

2.

Definition of Continuous Wigner and Bai Distribution

The WDF of a function $f(x)$ is defined as:

Assuming the Fourier transform of $f(x)$ is $F(k)$, then the WDF can be defined in the Fourier frequency domain alternatively by applying the Parseval formula as shown below:

In Fourier optics, any optical field can be decomposed into a set of plane wave components with different propagation directions, and the Fourier spatial spectrum of this field is closely related to the amplitude distribution of these components. By applying the Fourier transform to the IAF of a signal, the WDF can provide an intermediate representation between the space and spatial-frequency domain. This representation can characterize the energy of the rays or the plane wave components emitted from specific positions on the source plane and propagate with specific angles. It indicates that the WDF can link the wave optics with the ray concept in geometrical optics in free space.^8,14 In addition, the cross-terms of the WDF can characterize the interference or the coherence properties of the optical fields so that the WDF can also be applied in the partial coherence area.

The LCT of a function $f(x)$ for the LCT parameters $M=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)$ is defined as

Then, the BDF of a function $f(x)$ for the parameter $M=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)$ is defined by replacing the Fourier kernel $e^{-ik\tau }$ in Eq. (1) with the LCT kernel in Eq. (5), which is given by

For $M=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)$, Eq. (6) reduces to the WDF in Eq. (1). Furthermore, an alternative BDF definition in the Fourier frequency domain by applying Parseval’s theorem is presented below when $a=0$:

The BDF, such as the WDF, is separable in higher dimensions.

Various mathematical properties for the BDF have been determined, including the shift theorem, multiplication theorem, support, and nature of the cross-terms. The known properties are reviewed and extended in Ref. 14.

Similar to the WDF, by applying the LCT to a signal’s IAF, the BDF can provide an intermediate expression between the space $x$ and its linear canonical frequency domain $u$. It can interpolate between the IAF and the WDF of an optical field, illustrate how the IAF evolves through the propagation in free space, and eventually result in a WDF when the propagation distance is large. It can also reveal how the cross-terms interfere with each other during the propagation process and result in the fluctuated pattern in the WDF where the cross-terms can capture coherence information of an optical field.¹⁴ In addition, if the frequency marginal property of the BDF can be explained appropriately, the BDF may also extend the application of the WDF in linking the wave optics and geometrical optics. It may connect the positions on the source plane with not only the rays’ propagation angles but also the positions on the observation plane at any propagation distances in free-space propagation or after arbitrary first-order optical systems.

3.

Discrete Wigner Distribution Function

The BDF was proposed as a generalization of the WDF. To cite our discrete BDF in context, we briefly review the literature on discrete WDFs. Claasen and Mecklenbrauker¹⁵ first introduced a discrete-space WDF, $W_{CM,f}(n,k)$, of a discrete signal $f(n)$, which is given by

Chan¹⁷ proposed an alternative discrete-space WDF definition. Assume $L$ is an operator on the discrete signal $f(n)$ defined by

Then, the discrete-space WDF proposed by Chan, $W_{\text{Chan},f}(n,k)$, is given by

Although Chan claimed that this discrete-space WDF is free of aliasing, it was disproved soon afterward by Claasen and Mecklenbrauker.¹⁸ However, in this expression, the aliased replicas alternate between positive and negative values on even- and odd-numbered samples, so proper averaging processes can suppress them. This also leads to the fact that $W_{\text{Chan},f}(n,k)$ can satisfy more mathematical properties, especially the integral-type properties such as frequency marginal, compared with $W_{CM,f}(n,k)$. The zero-padding process in this method means that the calculation efficiency of this method is much worse than $W_{CM,f}(n,k)$. A discrete frequency version of $W_{\text{Chan},f}(n,k)$ can be obtained by sampling its frequency domain as well.¹⁶

Claasen and Mecklenbrauker extended Chan’s work to propose three new expressions of the discrete-space WDF to address the aliasing problem,¹⁸ which are given by

Another discrete-space WDF definition was proposed by Peyrin and Prost,¹⁹ $W_{PP,f}(n,k)$, which is given by

By sampling its frequency domain, the discrete WDF is expressed as

This definition is closely related to the expression proposed by Chan, $W_{\text{Chan},f}(n,k)$.¹⁶ Therefore, they have very similar characteristics, and the aliasing components in $W_{PP,f}(n,j)$ can be suppressed using averaging processes as well. Furthermore, Peyrin and Prost also proposed that instead of oversampling the signal or taking the averaging processes, an analytic signal sampled at normal Nyquist frequency can also be applied to solve the aliasing problem, which can be considered as an alternative method for oversampling a real signal.

Beiker et al.²⁰ and Chassante-Mottin et al.²¹ applied a similar averaging method proposed by Claasen and Mecklenbrauker¹⁸ to address the aliasing issue. One of their discrete WDF definitions is given by

Richman et al.²³ applied a group theoretic definition to propose a discrete WDF, and O’Neill et al.²⁴ presented the same definition using a different derivation method. Their discrete WDF is given by

This expression assumes $N$ is odd. Although it satisfies most desired mathematical properties, it still suffers from the aliasing problem.¹⁶

The proposed discrete or discrete-time WDF definitions reviewed above have different strengths and weaknesses. They are the candidates of an ideal discrete WDF, where the ideal discrete WDF should closely represent its continuous counterpart, satisfy all the mathematical properties, be alias-free under Nyquist sampling criteria, and be easy to evaluate. These aspects should also be considered when defining a discrete BDF. In 2009, O’Toole¹⁶ found that almost all these definitions in previous literature can be reformulated as the linear combinations of two types of discrete WDFs. These WDFs are based on the first and second sampling methods mentioned above proposed by Claasen and Mecklenbrauker,¹⁵ and Chan,¹⁷ respectively, $W_{CM,f}(n,j)$ and $W_{\text{Chan},f}(n,j)$. In addition, he further proposed two related discrete WDF definitions to reduce the computational efficiency problem that occurred in $W_{\text{Chan},f}(n,j)$. According to his results, although $W_{CM,f}(n,j)$ and $W_{\text{Chan},f}(n,j)$ are not the best or optimal definitions in this area, they have become the most significant parts of establishing the optimal discrete WDF. Therefore, to discretize the BDF, it is reasonable to adapt one of these sampling methods for making contributions to formulate an optimal discrete BDF.

4.

Discrete Bai Distribution Function

The following sections present the analysis of a discrete BDF mathematically and geometrically, which includes two steps. The first step is the discretization of the continuous BDF’s space domain $x$. In Sec. 4.1, we start by sampling the test signal $f(x)$ with a uniform sampling period $T$, then generate a discrete IAF based on Claasen and Mecklenbrauker’s sampling method.¹⁵ After that, by taking the discrete-space LCT of the sampled $R_f(x,\tau )$ with respect to $\tau$, a mathematical expression for a discrete-space and continuous-frequency BDF (discrete-space BDF) is established. This discrete-space BDF is periodic in its LCT frequency domain $u$. In Sec. 4.2, we first calculate the support of the BDF in the $u$ direction. Then, based on the periodicity and the support, the requirement for the sampling period $T$ is analyzed to avoid overlaps of the replicas created by sampling the signal. The second step is the discretization of the LCT frequency domain $u$ of the discrete-space BDF, which is presented in Sec. 4.3. This is done by making the sampled $R_f(x,\tau )$ chirp-periodic in the $\tau$ direction with an appropriate period. Then, a discrete-space and discrete-frequency BDF (discrete BDF) can be obtained by taking the discrete LCT of the sampled and chirp-periodic $R_f(x,\tau )$ with respect to $\tau$. Finally, we will show that the derived discrete BDF is the generalization of the discrete WDF proposed by Claasen and Mecklenbrauker.

4.1.

Discrete-Space Bai Distribution Function

In this section, we determine an expression for a discrete-space BDF.

Assume a function $f(x)$ has finite support as shown in Fig. 1(a). Its IAF $R_f(x,\tau )$ has support as shown in Fig. 1(b). It can be observed that the signal length in the $\tau$ direction is twice as large as in the $x$ direction because of the scaling effect on $\tau$.

Fig. 1

(a) Support of the continuous signal $f(x)$ on $x$. (b) Support of the continuous IAF of $f(x)$ on $x$ and $\tau$.

To obtain an expression of a discrete-space BDF, the $R_f(x,\tau )$ needs to be sampled first. In this paper, we first consider a naive and uniform sampling method; the sampling frequency in both $x$ and $\tau$ directions is the same, which is $\frac{1}{T}$, where $T$ is the sampling period as shown in Fig. 2(a). By substituting $x=nT$ and $\tau =mT$, the discrete version of the $R_f(x,\tau )$ using the naive sampling method, we obtain

Eq. (18)

$$R_f(nT,mT)=f(nT+\frac{m}{2}T)f^{*}(nT-\frac{m}{2}T),$$

where $m$ and $n$ are integers. However, this method will only be feasible when we can directly access and sample a continuous $R_f(x,\tau )$. Practically, the discrete IAF should be obtained from a discrete test signal $f(nT)$, which is impossible without interpolation techniques because we can only access the samples at $nT$ and cannot access the samples at non-integer multiples of $T$, i.e., for odd $m$. To address this issue, an alternative sampling approach proposed by Claasen and Mecklenbrauker, which we will refer to as the CM approach, is shown in Fig. 2(b) by simply ignoring the samples when $m$ is odd in the naive method.¹⁵ In this sampling method, the sampling period in the space domain $x$ remains $T$ but will be doubled in the $\tau$ domain, which is $2T$. The discrete version of the IAF using the CM sampling method is given by

Eq. (19)

$$R_f(nT,2mT)=f(nT+mT)f^{*}(nT-mT).$$

Fig. 2

(a) Naive sampling method. (b) CM sampling method.

Unlike the naive sampling approach $R_f(nT,mT)$, the CM approach $R_f(nT,2mT)$ can be generated directly from a discrete signal $f(nT)$.

After that, to obtain the discrete-space BDF, we simply need to take the discrete-space LCT of $R_f(nT,2mT)$. One algorithm to calculate this comprises four steps:²⁵

• 1. First chirp multiplication

• 2. Discrete-space Fourier transform

• 3. Scaling operation

• 4. Second chirp multiplication

• Step 1: Assume $g(nT,2mT)$ equals the product of $R_f(nT,2mT)$ and a sampled chirp signal (using the same sampling method as $R_f(nT,2mT)$),

  Eq. (20)

  $$g(nT,2mT)=R_f(nT,2mT)e^{\frac{i\pi a}{b}(2mT{)}^2}.$$

• Step 2: Take the discrete-space Fourier transform of $g(nT,2mT)$ with respect to $m$. The resultant $G(nT,k)$ will be discrete in $n$, and continuous and periodic with a period of $\frac{1}{2T}$ in $k$, i.e.,

  Eq. (21)

  $$G(nT,k)=\sum _p{G}_c(nT,k-\frac{1}{2T}p),$$
  where $G_c(nT,k)$ is the Fourier transform of $g(nT,\tau )$ on $\tau$. Note that we only sample $g(x,\tau )$ in space domain $x$ with a sampling period $T$ and $\tau$ remains continuous.

• Step 3: Scale the variable $k$ in Eq. (21) by a factor $b$ and substitute $u=kb$,

  Eq. (22)

  $$G(nT,u)=\sum _p{G}_c(nT,u-\frac{b}{2T}p).$$

• Step 4: Finally, the discrete-space LCT of $R_f(nT,2mT)$ with respect to $m$ can be obtained by a multiplication between the $G(nT,u)$ in Eq. (22) and a chirp signal $e^{\frac{i\pi d}{2b}{u}^2}$, which is defined as the discrete-space BDF ${\tilde{B}}_f(nT,u)$,

  Eq. (23)

  $${\tilde{B}}_f(nT,u)=e^{\frac{i\pi d}{b}{u}^2}\sum _p{G}_c(nT,u-\frac{b}{2T}p).$$

According to Eq. (23), the discrete-space BDF is chirp-periodic. The replicas are placed by $\frac{|b|}{2T}$ in its linear canonical frequency domain $u$. It can be observed that the period mainly depends on the value of $T$. Therefore, the requirement for the sampling period $T$ can be analyzed by preventing the replicas from overlapping with each other.

Then, we will show that the derived discrete-space BDF is a generalization of the discrete-space WDF proposed by Claasen and Mecklenbrauker. Recall that when we set the LCT matrix as $M_F=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)$, as mentioned in Sec. 2, the BDF reduces to the WDF. In such a situation, if we substitute this matrix into the discrete-space LCT calculation process shown above, the chirps in Steps 1 and 4, $e^{\frac{i\pi a}{b}(2mT{)}^2}$ and $e^{\frac{i\pi d}{b}{u}^2}$, simply become 1. In addition, the scaling factor in Step 3 becomes $b=1$ as well. Therefore, the resultant discrete-space BDF simply reduces to the discrete-space Fourier transform of the sampled IAF, $R_f(nT,2mT)=f(nT+mT)f^{*}(nT-mT)$, with respect to $m$, which is given by

Eq. (24)

$${\tilde{B}}_{f,M_F}(nT,k)=\sum _{m}f(nT+mT)f^{*}(nT-mT)e^{-i4\pi kmT}.$$

It can be observed that Eq. (24) has the same form as 7, which means the discrete-space WDF proposed by Claasen and Mecklenbrauker is a special case of the derived BDF.

4.2.

Sampling Requirement Analysis for the Discrete-Space BDF

In Eq. (23), we defined a discrete-space BDF, i.e., a BDF of the discrete signal $f(nT)$. Methods to determine an appropriate value of $T$ based on, e.g., the support of the Fourier transform of $f(x)$, are well known. However, it may be necessary to decrease $T$ so that the discrete-space BDF of $f(nT)$ well approximates the BDF of $f(x)$. To do that, we will use PSDs. The basic idea goes back at least as far as Adolf Lohmann in the 1960s.²⁶ Lohmann sketched figures marking the limits of a signal in phase space (space-spatial frequency or time-frequency, as appropriate) and used them to analyze problems such as the storage capacity of different kinds of holograms on a given piece of film. There are three foundational ideas here. One is that we can usefully bind a signal in phase space. In fact, this is a mild fiction. A signal that has compact support in space (for example) must by theorem have infinite bandwidth, see, e.g., Ref. 27. The PSD of such a signal would be an infinitely tall rectangle and require by Nyquist’s theorem an infinitely large sampling rate. For perfect reconstruction, this is true, but normally, we compromise with something that is arbitrarily good enough. This arbitrary sampling rate implies a pseudo-bandwidth in the Fourier domain, and this is the width in frequency depicted in a PSD. An example is shown in Fig. 3(b). The second key idea is that the effects of various optical transformations on a PSD are well known and in many cases are simple linear coordinate transformations, which lead to PSDs being used to analyze the sampling requirements of LCT algorithms based on the decomposition of the ABCD matrix, e.g., Refs. 1 and 28. This also leads to PSDs that assume (approximate) compact support in domains other than Fourier.²⁷ An example is shown in Fig. 3(a). Depending on the specific distribution of the signal, one or another representation may be more efficient (essentially produce a PSD with a smaller area), but this latter kind of PSD lends itself to certain analyses and can physically model, e.g., truncation by a camera aperture in a non-imaging plane. For example, consider Ref. 28. Whatever the assumptions we make, the point of a PSD is to rise above considering a specific signal and instead address the cohort of all signals with the same (approximate) compact support in space and some other plane. A third key idea here is that periodic replicas arise in frequency when we sample, which is clearly visible in a derivation of Nyquist’s theorem. These replicas also appear in WDFs and can be included in the PSD analysis of numerical algorithms, e.g., Ref. 28. There is a simplifying assumption here: cross-terms appear between all of the infinite replicas in the WDF or BDF of a discrete signal and are not accounted for in this analysis. We justify this based on three things: First, a study of the impact of replicas on the numerical approximation of the WDF shows that their effects are not large.⁹ Second, the results in this paper, shown in Fig. 14, show a typical case in which the effect of the replicas is much stronger than the effect of the cross-terms. This is shown by the substantial reduction in the slope of the error curve immediately after the sampling rate identified in this paper is met. Third, we note that our analysis, while performed using different tools, reduces to that of Claasen and Mechlenbraucher in the special case. We note finally that the goals of our analysis are twofold: (1) the replicas must not overlap, and (2) the replicas should tile without gaps so as not to waste samples representing no or negligible information.

Fig. 3

(a) PSD of the continuous test signal $f(\tau )$ based on its LCT width and (b) based on its Fourier bandwidth.

Our next step is to determine the width or support of the BDF in the $u$ direction. From this, we can determine the sampling requirement for the discrete-space BDF.

We use a PSD analysis as described above. We interpret the BDF here as the LCT of the $R_f(x,\tau )$ with respect to $\tau$. This means that we must analyze how the PSD of $R_f(x,\tau )$ for a specific $x$ value will evolve through the LCT and the sampling process to investigate the sampling requirement. Finally, we identify how the $x$ value will impact the PSD and the sampling requirement. We further decompose the problem into the components of the $R_f(x,\tau )$. This can be a little challenging to follow, so let us review the steps here: Given $f(x)$ and some assumptions about it, we will determine the following sequence of PSDs:

• 1. $f(x-\tau /2)$.

• 2. $f(x+\tau /2)$.

• 3. Hence, their product, giving $R_f(x,\tau )$ for fixed $x$.

• 4. Hence, the discrete-space LCT of $R_f(x,\tau )$, which is the discrete-space BDF.

We will do this in two different ways, differing in the assumptions about what information is known about the signal. Initially, we will assume that we know the (Fourier) bandwidth, and finally, we will assume that we know the width of the signal’s LCT with the same ABCD parameters as the BDF. The Fourier assumption is where most people start, but we show that it leads to an unsatisfying and inefficient discretization of the distribution.

Figure 3 shows the PSD of the two initial assumptions, with Fig. 3(a) being where we know the LCT’s extent (sometimes called, e.g., LCT-bandwidth, or ABCD-bandwidth) and Fig. 3(b) where we know the Fourier transform’s extent (i.e., the signal’s bandwidth). Both of them are PSDs of a continuous function $f(\tau )$, but the assumptions are different, and as discussed previously, both describe a finite region of phase space in which most of the signal’s energy is contained. For Fig. 3(a), the two sides that are parallel to the $k$ axis define the width of the signal in the $\tau$ domain, $W_{\tau }$. The other two sides with the same slope of $-\frac{4a}{b}$ are separated by a distance $|\frac{1}{b}|W_u$ in the $k$ direction, where $W_u$ represents the width of $f(\tau )$’s LCT when the parameter matrix is $M=\left(\begin{array}{cc}-4a& b\\ c& d\end{array}\right)$. For Fig. 3(b), the width in the $\tau$ direction is also $W_{\tau }$, and the two sides separated by a distance $W_k$ define the width of the signal in its Fourier frequency domain. Note that the dot in the corner is an aid to tracking the orientation of the PSD.

Then, by scaling the $f(\tau )$ by a factor of 2, the $f(\frac{\tau }{2})$ can be obtained. According to the scaling theorem of the WDF, the two PSDs of $f(\tau )$ in Fig. 3 will be extended in $\tau$ and compressed in $k$ directions as shown in Figs. 4(a) and 4(b). After that, we flip the PSD in Fig. 4(b) along the $\tau$ direction and take its conjugate to obtain the PSD of $f^{*}(-\frac{\tau }{2})$ as shown in Fig. 4(c). Note that the complex conjugate will flip a function’s PSD upside down. However, due to its rectangular shape in our case, this operation will only impact its orientation (the position of the dot).

Fig. 4

(a) PSD of $f(\frac{\tau }{2})$ based on $f(\tau )$’s LCT width. (b) PSD of $f(\frac{\tau }{2})$ based on $f(\tau )$’s Fourier bandwidth. (c) PSD of $f^{*}(-\frac{\tau }{2})$ based on $f(\tau )$’s Fourier bandwidth.

Next, by shifting $f(\frac{\tau }{2})$ [Fig. 4(a)] and $f^{*}(-\frac{\tau }{2})$ [Fig. 4(c)] by a distance $2x$ to different directions in $\tau$, the functions $f(\frac{\tau +2x}{2})=f(x+\frac{\tau }{2})$ and $f^{*}(-\frac{\tau -2x}{2})=f^{*}(x-\frac{\tau }{2})$ can be obtained. Their PSDs are shown in Figs. 5(a) and 5(b), respectively.

Fig. 5

(a) PSD of $f(x+\frac{\tau }{2})$ based on $f(\tau )$’s LCT width and (b) PSD of $f^{*}(x-\frac{\tau }{2})$ based on $f(\tau )$’s Fourier bandwidth.

To obtain the PSD of the IAF, $R_f(x,\tau )=f(x+\frac{\tau }{2})f^{*}(x-\frac{\tau }{2})$, we first assume that the WDFs of $R_f(x,\tau )$, $f(x+\frac{\tau }{2})$, and $f^{*}(x-\frac{\tau }{2})$ with respect to $\tau$ are $W_{R_f}(\tau ,k)$, $W_{f_1}(\tau ,k)$, and $W_{f_2}(\tau ,k)$. Then, we apply the multiplication theorem of the WDF to find their relationship, which is given by

Eq. (25)

$$W_{R_f}(\tau ,k)={\int }_{-\infty }^{\infty }{W}_{f_1}(\tau ,p)W_{f_2}(\tau ,k-p)\mathrm{d}p.$$

Hence, the WDF of $R_f(x,\tau )$ is the convolution between $W_{f_1}(\tau ,k)$ and $W_{f_2}(\tau ,k)$ with respect to $k$, and the PSD of $R_f(x,\tau )$ is the support of $W_{R_f}(\tau ,k)$. To calculate this convolution, $W_{f_2}(\tau ,k)$ is flipped and shifted along the $k$ direction, multiplied by $W_{f_1}(\tau ,p)$, and integrated. Figures 6(a) and 6(b) show two extreme conditions for this shift process. The width of the overlapping area between the supports of $W_{f_1}(\tau ,p)$ and $W_{f_2}(\tau ,k-p)$ in the $\tau$ direction is ${\widehat{W}}_{\tau }=2(W_{\tau }-2x)$, which defines the support of $W_{R_f}(\tau ,k)$ in $\tau$. The resultant PSD will be a parallelogram, and the slopes of the upper and lower sides are the same as in Fig. 5(a), which is $-\frac{a}{b}$.

Fig. 6

(a) and (b) Two extreme conditions of the shift process in the convolution.

Then, the two conditions shown in Figs. 7(a) and 7(b) define the width of the two vertical lines of $R_f(x,\tau )$’s PSD, ${\widehat{W}}_u$, where $d_1=\frac{1}{2}(|\frac{1}{2b}|W_u+|\frac{2a}{b}|W_{\tau })-|\frac{2a}{b}|W_{\tau }+\frac{1}{4}{W}_k$ and $d_2=\frac{1}{2}(|\frac{1}{2b}|W_u+|\frac{2a}{b}|W_{\tau })+\frac{1}{4}{W}_k$. Figure 7(a) shows the situation when the lower side of the rectangular just reaches the upper-right corner of the parallelogram. Therefore, the IAF’s PSD is shown in Fig. 7(c), and the widths of the lines parallel to the $k$ direction are ${\widehat{W}}_u=d_1+d_2=\frac{1}{2}{W}_k+|\frac{1}{2b}|W_u$. Although this PSD is not centralized at the coordinate origin, this would not impact the following sampling analysis for the discrete BDF.

Fig. 7

Conditions (a) and (b) define the vertical lines’ width of $R_f(x,\tau )$’s PSD. (c) The PSD of the IAF, $R_f(x,\tau )$.

Now that we have obtained the PSD of the continuous $R_f(x,\tau )$, we can analyze the discretization process introduced in Sec. 4.1. For the CM sampling method, the sampling periods in the $x$ and $\tau$ directions are $T$ and $2T$, respectively. The sampling of the variable $\tau$ will result in a periodic replication in the $k$ direction of its PSD with a period $\frac{1}{2T}$ as shown in Fig. 8(a) and a substitution of $\tau =2mT$. Then, by applying the four-step discrete-space LCT calculation method to $R_f(nT,2mT)$ with respect to $m$, we can get the discrete-space BDF. Figures 8(b)–8(e) present how the PSD evolves through those four steps.

• Step 1: The first chirp multiplication will introduce a skew in the $k$ direction as shown in Fig. 8(b). The signal width in the $m$ direction remains unchanged, but the bandwidth in the frequency domain $k$ becomes ${\widehat{W}}_u+|\frac{a}{b}|{\widehat{W}}_{\tau }-|\frac{a}{b}|{\widehat{W}}_{\tau }={\widehat{W}}_u$, where $a$ and $b$ are LCT parameters.

• Step 2: The Fourier transform will introduce a 90 deg rotation as shown in Fig. 8(c). Both ${\widehat{W}}_u$ and ${\widehat{W}}_{\tau }$ remain unchanged.

• Step 3: For the scaling operation, we substitute $u=kb$, then the horizontal axis becomes the LCT frequency axis $u$ in the BDF, and this will impact both widths in the two directions and the periodicity. In the $u$ direction, the period becomes $\frac{|b|}{2T}$, and ${\widehat{W}}_u$ changes to $|b|{\widehat{W}}_u$ as shown in Fig. 8(d). Note that the sign of the parameter of $b$ can only impact the orientation of the PSD in the $u$ direction but not the width.

• Step 4: Finally, the second chirp multiplication will introduce another skew in the $m$ direction, which will not impact the width in the $u$ direction, and the PSD of the discrete-space BDF for specific $n$ values is obtained as shown in Fig. 8(e). Furthermore, the period in the $u$ direction is exactly the same as in Eq. (23), which is $\frac{|b|}{2T}$.

Fig. 8

(a) PSD of $R_f(nT,2mT)$. PSD of $R_f(nT,2mT)$ (b) after the first chirp multiplication, (c) after the discrete-space Fourier transform, (d) after the scaling operation, and (e) after the second chirp multiplication.

In Fig. 8(e), it can be observed that all the terms are orthogonal with each other, and the central term can be extracted easily. However, with the increase of the sampling period, the terms will overlap with each other and lead to aliasing. Therefore, to avoid the overlapping problem, the period $\frac{|b|}{2T}$ should always be larger than the width or support of the BDF in the $u$ direction $|b|{\widehat{W}}_u$. Then, we have $\frac{|b|}{2T}\ge |b|{\widehat{W}}_u$, where $|b|{\widehat{W}}_u=|b|(|\frac{1}{2b}|W_u+\frac{1}{2}{W}_k)$. Substitute ${\widehat{W}}_u$ into the inequality, $\frac{|b|}{2T}\ge |b|(|\frac{1}{2b}|W_u+\frac{1}{2}{W}_k)$, and the sampling requirement is given by

Eq. (26)

$$f_s\ge W_k+|\frac{1}{b}|W_u,$$

where $f_s=\frac{1}{T}$ is the sampling frequency.

4.3.

Discrete Bai Distribution Function

In previous sections, we have already derived the expression for the discrete-space BDF and analyzed its sampling requirement based on the PSD method. In this section, we will show the derivation process of a discrete BDF by applying the same method. We start by constructing a chirp-periodic version of the IAF, $R_f(x,\tau )$, to tile the terms in the discrete BDF without gaps. The first step is multiplying $R_f(x,\tau )$ by a chirp signal, $e^{i\frac{\pi a}{b}{\tau }^2}$, to cancel the skew in the $k$ direction of its PSD. The resultant PSD is shown in Fig. 9(a), and the resultant expression is given by

Eq. (27)

$$R_f(x,\tau )e^{i\frac{\pi a}{b}{\tau }^2}.$$

Fig. 9

PSD of the IAF (a) after multiplying by a chirp and (b) after periodizing.

Then, by making Eq. (27) periodic in the $\tau$ direction with a period $\frac{|b|}{\mathrm{\Lambda }}$, we obtain

Eq. (28)

$$\sum _p{R}_f(x,\tau -\frac{|b|}{\mathrm{\Lambda }}p)e^{i\frac{\pi a}{b}(\tau -\frac{|b|}{\mathrm{\Lambda }}p{)}^2}.$$

The PSD calculated by Eq. (28) is shown in Fig. 9(b).

Next, by multiplying the PSD shown in Fig. 9(b) by a chirp with opposite phase distribution, $e^{-i\frac{\pi a}{b}{\tau }^2}$, a chirp-periodic IAF, $\widetilde{R_f}(x,\tau )$, is obtained as shown in Fig. 10(a), and the term from each period is exactly the same as the PSD of IAF [Fig. 7(c)]. Its expression is given by

Eq. (29)

$$\widetilde{R_f}(x,\tau )=\sum _p{R}_f(x,\tau -\frac{|b|}{\mathrm{\Lambda }}p)e^{i\frac{\pi a}{b}(\tau -\frac{|b|}{\mathrm{\Lambda }}p{)}^2}{e}^{-i\frac{\pi a}{b}{\tau }^2}.$$

Fig. 10

(a) PSD of the chirp-periodic IAF. Its PSD (b) after the sampling process, (c) after the first chirp multiplication, (d) after the DFT, (e) after the scaling operation, and (f) after the second chirp multiplication.

By applying the CM sampling method to the chirp-periodic IAF, $x=nT$ and $\tau =2mT$, its PSD becomes periodic in the $k$ direction with a period $\frac{1}{2T}$ as shown in Fig. 10(b). Furthermore, the discrete version of $\widetilde{R_f}(x,\tau )$ is given by

Eq. (30)

$$\widetilde{R_f}(nT,2mT)=\sum _p{R}_f(nT,2mT-\frac{|b|}{\mathrm{\Lambda }}p)e^{i\frac{\pi a}{b}(2mT-\frac{|b|}{\mathrm{\Lambda }}p{)}^2}{e}^{-i\frac{\pi a}{b}(2mT{)}^2}.$$

After that, by taking the discrete LCT of ${\tilde{R}}_f(nT,2mT)$, a discrete BDF can be calculated. This can be done by modifying the LCT calculation method introduced in Sec. 4.1.

• Step 1: Multiply ${\tilde{R}}_f(nT,2mT)$ with a sampled chirp signal,

  Eq. (31)

  $$\begin{gathered} \tilde{g}(nT,2mT)=\sum _p{R}_f(nT,2mT-\frac{|b|}{\mathrm{\Lambda }}p)e^{i\frac{\pi a}{b}(2mT-\frac{|b|}{\mathrm{\Lambda }}p{)}^2}{e}^{-i\frac{\pi a}{b}(2mT{)}^2}{e}^{i\frac{\pi a}{b}(2mT{)}^2} \\ =\sum _p{R}_f(nT,2mT-\frac{|b|}{\mathrm{\Lambda }}p)e^{i\frac{\pi a}{b}(2mT-\frac{|b|}{\mathrm{\Lambda }}p{)}^2}. \end{gathered}$$

  Its PSD will skew in the $k$ direction. However, according to Eq. (31) and its PSD shown in Fig. 10(c), the constructed chirp-periodic structure cancels this skew and the shifts of the terms in the $k$ direction caused by the chirp multiplication. It can also be observed that we have successfully obtained a periodic version of Eq. (20) in Sec. 4.1 so that we can apply the DFT to it in the next step instead of the discrete-space Fourier transform.

• Step 2: Take the DFT of $\tilde{g}(nT,2mT)$ with respect to $m$. The resultant $G(nT,\frac{\mathrm{\Lambda }}{|b|}j)$ will be discrete in both directions, i.e.,

  Eq. (32)

  $$G(nT,\frac{\mathrm{\Lambda }}{|b|}j)=\sum _p{G}_c(nT,\frac{\mathrm{\Lambda }}{|b|}j-\frac{1}{2T}p),$$
  where $G_c(nT,k)$ is the Fourier transform of $R_f(nT,\tau )e^{\frac{i\pi a}{b}{\tau }^2}$ with respect to $\tau$ and $\frac{\mathrm{\Lambda }}{|b|}$ is the sampling frequency of the variable $k$. Furthermore, its PSD will rotate by 90 deg clockwise as shown in Fig. 10(d).

• Step 3: Scale the $\frac{\mathrm{\Lambda }}{|b|}j$ in Eq. (32) by a factor $|b|$,

  Eq. (33)

  $$G(nT,\mathrm{\Lambda }j)=\sum _p{G}_c(nT,\mathrm{\Lambda }j-\frac{|b|}{2T}p).$$

  The width of the PSD in both directions will change as shown in Fig. 10(e), and the sampling frequency in the $j$ direction becomes $\mathrm{\Lambda }$.

• Step 4: Finally, the discrete LCT of $R_f(nT,2mT)$ with respect to $m$ can be obtained by a multiplication between the $G(nT,\mathrm{\Lambda }j)$ in Eq. (33) and a discrete chirp signal $e^{\frac{i\pi d}{b}(\mathrm{\Lambda }j{)}^2}$, which is defined as the discrete BDF ${\tilde{B}}_f(nT,\mathrm{\Lambda }j)$. Its PSD is shown in Fig. 10(f),

  Eq. (34)

  $$\widetilde{B_f}(nT,\mathrm{\Lambda }j)=e^{\frac{i\pi d}{b}(\mathrm{\Lambda }j{)}^2}\sum _p{G}_c(nT,\mathrm{\Lambda }j-\frac{|b|}{2T}p).$$

The calculation process for a discrete BDF is established as above, and the sampling frequency in the $x$ and $u$ directions is $T$ and $\mathrm{\Lambda }$, respectively. It can be observed that all the terms in the PSD of the discrete BDF in Fig. 10(f) are orthogonal with each other in both directions so that the aliasing problem can be avoided by applying the sampling requirement proposed in Sec. 4.2, which is $f_s\ge W_k+\frac{1}{|b|}{W}_u$.

Finally, we will show the derived discrete BDF is the generalization of the discrete WDF proposed in the previous literature. If we set the LCT parameter as $M_F=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)$, the term $\widetilde{R_f}(nT,2mT)$ in Eq. (30) simply reduces to a sampled version of periodic IAF without any chirp multiplication, $\sum _p{R}_f(nT,2mT-\frac{|b|}{\mathrm{\Lambda }}p)$. Similarly, the chirp signals in Steps 1 and 4 also reduce to 1, and the scaling factor in Step 2 becomes $b=1$ as well. Therefore, in such a situation, the discrete BDF, ${\tilde{B}}_{f,M_F}(nT,\mathrm{\Lambda }j)$, becomes the DFT of the periodic and sampled IAF, $\sum _p{R}_f(nT,2mT-\frac{|b|}{\mathrm{\Lambda }}p)$. If we assume the length of one period in the $m$ direction is $N$, then this DFT can be expressed as

Eq. (35)

$${\tilde{B}}_{f,M_F}(nT,j\mathrm{\Lambda })=\sum _{m}f(nT+mT)f^{*}(nT-mT)e^{-i4\pi j\mathrm{\Lambda }\frac{mT}{N}}.$$

It can be observed that the above equation has the same form as Eq. (8), which means the discrete WDF proposed by Claasen and Mecklenbrauker¹⁵ is the special case of the derived discrete BDF.

5.

Simulation Results

In this section, several simulation results of the discrete BDF for a test signal with different LCT parameters will be presented and compared with its analytical results to verify the correctness of the sampling requirement derived in Eq. (26). The test signal $f(x)$ is chosen as an addition of two Gaussian functions for its ease of integration and distinctive shape, which is given by

To calculate its BDF sampling requirement when the LCT parameters are $a=1$, $b=3$, and $d=1$, the signal $f(x)$’s Fourier transform and its LCT when $a=4$, $b=3$, and $d=1$ are shown in Figs. 11(a) and 11(b). According to Figs. 11(a) and 11(b), the signal’s bandwidth $W_k$ is 1, and its LCT width when $M=\left(\begin{array}{cc}4a& b\\ c& d\end{array}\right)$, $W_u$, is 36. By applying the sampling requirement proposed in Sec. 4.2, $f_s\ge W_k+\frac{1}{b}{W}_u$, the sampling frequency $f_s$ should be larger than $1+\frac{1}{3}\times 36=13\mathrm{Hz}$.

Fig. 11

(a) Spectrum of the Gaussian test signal. (b) The LCT of the Gaussian test signal when $a=4$, $b=3$, and $d=1$.

By applying this sampling method, the simulated discrete BDF and the analytical BDF of the test signal $f(x)$, as well as the error between them, are shown in Figs. 12(a)–12(c), where the analytical result is obtained by directly sampling a continuous BDF calculated using Mathematica. This Mathematica Notebook file has been uploaded to our GitHub repository.²⁹ It can be observed that the simulated result basically agrees with its analytical value, and the error is under ${10}^{-3}$, which is acceptable. If we reduce the sampling frequency to 6.5 Hz, interference will appear on the simulated BDF, and the error will become much larger than before as shown in Figs. 13(a)–13(c). The BDF is in general complex; absolute values are plotted here, but similar agreement is reached for amplitude/phase or real/imaginary plots, which are omitted for brevity.

Fig. 12

(a) Simulated discrete BDF when $a=1$, $b=3$, $d=1$, and $f_s=13\mathrm{Hz}$. (b) Analytical discrete BDF. (c) Error between the simulation and analytical solution.

Fig. 13

(a) Simulated discrete BDF when $a=1$, $b=3$, $d=1$, and $f_s=6.5\mathrm{Hz}$. (b) Analytical discrete BDF. (c) Error between the simulation and analytical solution.

Fig. 14

Error curve of the discrete BDF for different sampling frequencies and different $b$ values.

When $a=1$ and $d=1$ remain unchanged, the mean-squared error for different $b$ values at different sampling frequencies is shown in Fig. 14. It can be observed that each error curve for each $b$ value decreases rapidly first until it reaches a certain range, and the accuracy will improve much slower after that transitional range. The error of the BDF with a higher $b$ value decreases faster than the BDF with a lower $b$ value, which means the BDF with a higher $b$ value requires a lower sampling frequency to maintain the accuracy. The required sampling frequencies for $b=0.75$, $b=1$, $b=1.5$, $b=3$, and $b=10$ are 43.5, 33, 23.6, 13, and 5 Hz, respectively, and these values have been marked using dots with corresponding colors in Fig. 14. It can be seen that all the proposed sampling frequencies are higher than and very close to the transitional point, which means that the proposed sampling requirement can maintain the accuracy of the discrete BDF as well as the computational efficiency.

Furthermore, when $b=10$ and $d=1$ remain unchanged, the required sampling frequencies for $a=10$, $a=7.5$, $a=5$, $a=2.5$, and $a=1$ are 36.8, 26.8, 18.8, 9.8, and 4.8 Hz, respectively. The mean-squared error for different $a$ values at different sampling frequencies is shown in Fig. 15, and the proposed sampling frequencies have been marked using dots with corresponding colors as well.

Fig. 15

Error curve of the discrete BDF for different sampling frequencies and different $a$ values.

6.

Conclusion

Recently, considerable efforts have been devoted to developing novel tools in optical signal processing, especially applying the LCT. Based on our previous analysis,¹⁴ the BDF stands out as the most promising among these tools. It has been demonstrated to be useful for extracting linear frequency modulation parameters, and we have shown that it can model partially coherent systems. To apply this tool in simulation and to digital data, it is crucial to establish its discrete version and corresponding calculation method and code. These will enable people to investigate further potential applications and apply it to experimental data.

We have presented in this paper a sampling analysis of the BDF and established a framework for calculating the BDF of a discrete signal. Our analysis is different, but we arrive at a discrete BDF that reduces to the discrete WDF proposed by Claasen and Mecklenbrauker¹⁵ in the special case. This is a particularly important discretization as most of the discrete WDFs in the literature are linear combinations of it and one other.

Our sampling analysis indicates a sampling rate appropriate for a signal for our discrete BDF to approximate the continuous case. Our analysis ignores the effects of cross-terms, but these are known to have a modest impact on similar calculations.⁹ This is backed up by the error curves in Figs. 14 and 15, which show greatly reduced (though non-zero) improvements for sampling rates over and above those identified by our analysis. Finally, the consistency with Claasen and Mecklenbrauker’s discrete WDF suggests we have made reasonable assumptions.

Our discrete BDF uses a chirp-periodic IAF, which is known from the theory of discrete LCTs.³⁰

In the future, a deeper investigation of the BDF’s applications will be conducted, and its robustness for experimental data will also be examined.