A fast least-squares migration with ultra-wide-band phase-space beam summation method

Least-squares migration (LSM) problems are often formulated as iterative schemes. At each iteration, traditional LSM methods require solving the wave equation several times to compute the gradient and single scattering (Born) predicted data, which makes LSM computationally expensive. We use a beam approach to overcome this problem, based on the ultra-wide-band phase space beam summation (UWB-PS-BS) method. We use the beams to expand and propagate the measured data. The beams calculation is performed via a small number of Green’s functions (GFs). These GFs are calculated once and then stored in the computer’s disk memory. Thus, wave equation calculations are avoided within each iteration. We use a beam transformation of shot gathers which can be considered as the data obtained by a pair of source and receiver beams. This data has a physical interpretation. The beam amplitudes extract the medium reflectivity generated from a localized region in subspace at a specific angle (local Snell’s law). This data calculation is performed before the iterative scheme. Thus, we may a-priori threshold beams with low amplitude beams, as they contain less essential information. Using the beam propagators spectral locality, we may further reduce the number of beams used to construct the image. In this work, we establish a new local beam-based data domain LSM in the frequency domain (FD). We use the properties described above to reduce the computational complexity of the LSM problem. We demonstrate the advantages of the proposed LSM algorithm via numerical examples.


Introduction
Migration techniques are important tools to recover subsurface structures from observed seismic data [1]. These formulations approximate the adjoint of the forward modeling operation. An image obtained by the migration operation contains errors due to inaccuracy of the estimated velocity model, data limitation, and noise [2]. The resulting image can be verified by comparing the predicted data obtained from the migrated image, to the observed data. The error caused by the migration operation can be reduced by formulating the migration problem as a least-squares inverse problem [3,4]. Given a sufficiently accurate velocity model, this problem is often formulated as least-squares migration (LSM) [5].
LSM approaches are linearized inversion methods which minimize the misfit between the observed data and estimated data. Early implementation of LSM has been demonstrated in [5], using crosswell data. In [2], this approach has been applied to surface data, where Kirchhoff operators have been used to propagate the wavefield. It has been shown that the use of the LSM method improves the resolution, and reduces migration artifacts. The method has then been extended to a one-way wave equation [6,7]. To enhance the image quality, two-way wave equation and reverse-time LSM methods have been explored [6,8,[9][10][11].
Generally, LSM can be implemented in either image domain or data domain. The image domain LSM is based on explicitly constructing the Hessian matrix [12][13][14][15]. The LSM result is obtained by applying the inverse Hessian on a pre-computed migration image (also known as the gradient of the misfit function). Constructing the full Hessian matrix is computationally expensive for realistic problems. Data domain LSM utilizes iterative schemes to reduce the misfit function, where the explicit Hessian calculation is avoided. In this paper, we will employ the data domain LSM scheme.
In data domain LSM formulations, L 2 norm misfit function is often minimized by utilizing gradient-based optimization methods. The data domain LSM approach comprises of two steps: a migration part where the gradient of the misfit function is calculated via an adjoint-state method, and a Born modeling (demigration) part that generates estimated data with Born linear forward mapping [12]. These two steps make the implementation of LSM computationally expensive. In a time domain LSM implementation, wavefield propagations are necessary to compute Born data and the gradient. The number of wavefield propagations can be large when there is a large number of shot records. The computational cost of time domain LSM formulations often greatly increases with an increase in the amount of data. To reduce the computational cost, phase encoding formulations have been suggested in [16][17][18][19][20].
Frequency domain (FD) LSM methods have been studied by [21][22][23][24][25][26]. It has been shown that only a subset of the entire frequency components is needed to successfully recover subsurface images. The use of a frequencydependent grid has been proposed in [27]. Using this approach the computational complexity of the LSM problem is reduced. In [25] it has been shown that the Born data and the gradient can be obtained by using precomputed source and receiver Green's functions (GFs) that are stored in the computer's disk or memory. While this scheme works well for small datasets, implementing the scheme for large datasets becomes infeasible since the number of needed GFs would be tremendously large. Therefore, the adjoint-state method is often implemented for shot domain frequency LSM. In [26], a FD double plane wave LSM method has been proposed to improve the computational efficiency of the FD LSM problem. This method is based on the double plane wave transform of [28,29] (see also [30] and [31]) for the problem of reverse time migration. They have shown that only a small number of plane wave GFs are needed to perform LSM, which substantially reduces the wavefield calculations in the iterative updating process. The method works well for large datasets as well.
Beam summation (BS) methods play an important role in the study of wave propagation. In these formulations, the propagated field is represented as a superposition of beams emanating from a source, propagating along certain directions in space. The beams can be tracked locally in an inhomogeneous medium. In addition, these formulations are localized in space, since only beams that pass near a subsurface point need to be taken into account in the field calculation [32]. Due to these advantages beam-based migration techniques have been successfully applied [33][34][35][36][37][38]. To increase the image resolution and reduce the computational cost of LSM formulations, beam-based LSM approaches have been suggested in [39][40][41][42]. These beam-based formulations are generally based on point source expansions methods [43][44][45], where the source and receiver GFs are expanded using superposition of beams. Due to the spectral locality of beam propagators, only beams that contribute to an imaging point need to be used in the LSM, which makes it efficient for target-oriented LSM approaches.
There is another category of beam summation method based on extended (distributed) sources 1 . In the extended-source BS formulations, the field is expressed as a discrete superposition of shifted and tilted beams that emanate from a discrete set of points with initial ray parameters over the surface (see [46], figure 1(b)). The method has been structured initially upon a Gabor-series expansion of surface sources [47][48][49]. It has then been extended by using Gabor frames in the context of radiation from surface field distribution [50][51][52][53] and from volume sources [54,55]. For the seismic survey, the field is measured over a wide aperture. For this reason, we use the extended sources approach. Specifically, we use the ultra-wide-band phase space beam summation (UWB-PS-BS) method [51]. In this method, one utilizes windowed Fourier transform (WFT) frame decomposition of a surface field distribution to propagate the surface field. The frame elements are Gaussian windows, centered in the spatial domain around surface points, and the spectral domain, around a set of ray parameters. The propagated field is given as a superposition of weighted-beams, which emanate from the surface. In this theory, the beam amplitudes are calculated by projecting the surface field distribution on a set of dual frame functions. Due to the frame functions' structure, the projection operation extracts the local radiation spectrum (windowed Fourier transform) of the surface field near a surface point and specific direction (rayparameter). In [46] the UWB-PS-BS method has been used for seismic migration applications. The image is given as the sum of the source and receiver pair of beams and beam transformation of the shot gathers. Using the physical interpretation of the beam amplitudes (local radiation spectrum), it has been shown that the number of beams used for migration can be a priori reduced, and the migration method can be used to image specific targets over the subsurface.
In this research, we formulate the UWB-PS-BS migration method into LSM scheme. In our approach the beam propagators are solutions for the wave equation. They are calculated using numerical integration with subsurface GFs. The number of GFs needed to construct the beam propagators is relatively small. One can precompute FD GFs and store those in the computer's disk or even memory. Thus, in the iterative updating process, solving the two-way wave equation is no longer necessary, pre-computed GFs can be repeatedly used to construct the beam propagators and hence update images. Comparing with shot domain LSM methods, we significantly reduce the computational cost of the LSM implementation by eliminating the time-consuming wave-equation solving process. Additionally, since original shot gathers are decomposed into spatial and spectral localized beam domain data, one may filter out low amplitude beams that correspond to insignificant scattering data. With reduced input beam data, we should obtain an image close to the one received by migrating all the beam data but with reduced computational effort. Comparing with traditional shot domain LSM methods, the proposed beam-based LSM greatly improves the computational efficiency. is computationally infeasible, we often formulate the inverse problem as an iterative scheme. We can iteratively minimize the difference between observed data and predicted data via misfit or objective function with F is the misfit function, the superscript † represents the conjugate transpose operator, and the superscriptd enotes FD quantities. In equation ( ( ) represents the gradient of the misfit function. The best set of parameters is obtained when the misfit function reaches its minimum.

Frequency domain data under the Born approximation
Considering a 3D coordinate frame r = (x, z) = (x 1 , x 2 , z) with velocity v(r), the FD acoustic wave equation due to a point source illumination at x = x s with frequency signature of fˆis given by Here ω is the angular frequency, û is the propagated field and δ denotes the Kronecker delta function. Given a sufficiently accurate background velocity model, v B , the velocity and field can be written as Here Δ v is the velocity perturbation while u r B ( ) and u r ŝ ( ) are the background and scattered fields, respectively. The background field is given by . For this case, the scattering data collected at x = x r can be written as Thus, the weak scattering approximation of the propagation operator L of equation (1) is given by The gradient of the misfit function is given by [56] ] being the gradient of the misfit function with respect to the medium reflectivity. Note that the gradient g is in the length of number of model points r. Writing it explicitly in terms of equation (7), Typically a preconditioner which is usually the approximation of the Hessian matrix is applied to the gradient to improve the convergence of the iterative scheme. Following [57], the approximation of the main diagonal is chosen as This preconditioner is referred to as the illumination compensation.

Data domain least-squares migration using WFT frame expansions
In this section, we introduce our beam-based LSM method. The formulation is based on the UWB-PS-BS method of [51]. The main details for UWB-PS-BS theory are summarized in appendix.

Beam representation of the scattered data
Given measured or predicted shot gather data d x x , s r 0 ( ), the beam domain representation of the data is given by [46], equation (22) ò ò , which corresponds to the source and receiver arrays, respectively.
The Born approximation of the scattering data is obtained by first multiplying both sides of (6) by j μ (x s ) and j m¢ x r * ( ), and then changing the order of integration, The operation in equation (11) projects the shot gathers data on localized spatial-spectral window functions. Following equation (A.6) this operation identified as WFT of the data near the phase-space region ). The Born approximation of this data is given by equation (12). Following the discussion after [46], equation (22) the beam representation of the data has an analytic model. It is given by the multiplication of the imaging kernels and the reflectivity of the medium. Each beam propagator propagates along a certain direction in space, and decays away from it. The multiplication Y Y m m¢ r r b* ( )ˆ( ) is localized around a scattering cell (see also [46], figure 2). In addition, since each propagator propagates along a certain direction defined by the initial ray-parameter, each imaging kernel is oriented along the direction defined by the incident and scattered beams ray parameters. It follows from this relation that the Born approximation of the beam domain data is given by the local reflectivity of the subsurface, in an angle defined by the incident and backpropagated beams ray parameters.

Phase-space representation of the data misfit and the adjoint operator
In this work, we use a conjugate gradient algorithm to update the model. To use the conjugate gradient method, with With this normalization, the Born data of equation (12) can be written as

Numerical example
In this section we demonstrate the proposed LSM algorithm with the 2D Marmousi model.

The physical configuration
The model contains 460 and 150 grid points in the horizontal and vertical directions, respectively, with 0.02 km grid spacing in both directions. The velocity model is plotted in figure 1(a). We subtract a slightly smoothed velocity model from the true velocity model to obtain the reflectivity model plotted in figure 1(b).

Choosing the beam parameters
The first step in selecting the beam parameters is to set the beam frequency independent parameter, b. As demonstrated in [51], this parameter is chosen such that we have collimated beam propagators through the propagation domain, while also maintaining w b v 1 min 0  . Following [46] we set b = 20 km which leads to collimated beams for the subsurface. Note that the beam width is determined by [34]. We use the frequency band f = [5,25] Hz, with 0.25 Hz frequency interval for imaging. Following equation (A.10) we obtain the phase space lattice dimensions » = x p 0.6 km, 0.02 s km¯. We also obtain 16 beam initiation points with 43 propagation directions, giving the ray-parameter range p n ≈ [ − 0.42, 0.42] s/km. As illustrated below, we do not have to use all the beams for the proposed LSM to achieve good results.

Imaging with full data set and all GFs
We begin by exploring the proposed LSM algorithm using Born data as the observed data, which is calculated using equation (15). For forward modeling, we use a Ricker wavelet with central frequency f = 15 Hz. We use the frequency range f = [5,25] Hz, with Δf = 0.25 Hz. The beam propagators are calculated using equation (A.5). We use here numerical integration with pre-computed GFs. We use the frequency domain finite difference method to generate 92 GFs, with a horizontal spacing of 0.1 km. Note that the GFs spacing does not satisfy the horizontal spacing Nyquist sampling rate over the surface for frequency f > 7.5 Hz. This relatively small number of GFs allows us to calculate the GFs only once and store them in the computer's memory. Hence no wave propagation is needed at each iteration.
The image obtained by the standard beam-based migration algorithm of [46] is plotted in figure 2(a). As can be seen, most of the subsurface is reconstructed. The sides of the image, especially at the bottom parts are not clearly recovered. In figure 2(b) we plot the image obtained by the proposed LSM algorithm after 30 iterations.
As can be seen migration artifacts are reduced, and the image resolution is substantially improved. Figure 2(c) shows the convergence of the misfit function. We also plot the reflectivity profiles at x = 4.6 km in figure 2(d), to get a better comparison. We also plot the migrated image of [46] over this line. The proposed LSM algorithm successfully retrieves the true reflectivity value for most parts of the model.
In figure 3 we compare the input data and the estimated data obtained by the estimated reflectivity model after 30 iterations. We plot the data obtained by a source located at x = 4.6 km. For clarity all plots are scaled with the same colorbar. The data obtained by the proposed method show a good agreement with the input data.

Imaging with a reduced number of GFs
The number of GFs used for constructing the beam propagators can be reduced. In figure 4(a) and 4(c) we plot the migrated image and LSM image obtained using 46 GFs, with a spacing 0.2 km. As can be seen, the iterative LSM scheme helps to suppress migration artifacts and balance the reflectivity value for the entire model even though with a small number of GFs. In figure 4(b) we plot the image obtained by the beam-based migration algorithm, where we use 23 GFs to construct the beam propagators. For this case, the GFs spacing over the surface is Δx = 0.4 km while the Nyquist sampling rate is 0.15 km for = f f min . Therefore, the image obtained by the beam-based migration algorithm produces an image with poor resolution and spatial artifacts. In figure 4(d), we plot the image obtained by our LSM algorithm. Using the iterative scheme we manage to reconstruct most parts of the subsurface. The spatial resolution is greatly improved with a very small number of GFs. On the other hand, the image quality with 23 GFs is much degraded compared with that of 46 GFs mainly due to the fact that the Nyquist theorem is severely violated. In practice, one might want to use a reasonable size of GFs to reduce the computation cost in the meantime have most of the frequency components satisfy the Nyquist theorem.

Imaging with sparse data
The use of the UWB-PS-BS method also allows us to a-priori reduce the computational cost by using only the data that significantly contribute to the subsurface image. Following the discussion after equation (12), the beam domain data is governed by the subsurface local Snell's law reflection. Thus, we may neglect imaging kernels with low amplitudes as they do not contain significant scattering events. In [46] this relation was extended for the problem of scattering under the Born approximation, where it was shown that a pair of beams with high amplitudes corresponds to reflection events over the subsurface (see discussion after equation (12) and [46], equation (20)). Here, we use this property to reduce the number of elements that construct the propagation operator.
In figure 5(a) we plot the absolute value of the beam domain data generated with equation (15) at frequency f = 10 Hz. The data is plotted in the source and receiver beams domain. Each axis is arranged as a concatenated vector with each initiation point x m and all the propagation directions p n . The data is tagged in colors from the brightest to the darkest. The data set is sparse, and we can use only the pair of beams with high amplitudes, as explained above, to carry out the proposed LSM algorithm. Note that the beam amplitudes calculation is performed over the surface field, as a pre-processing phase, before the beam propagators computation. Thus, one may a-priori reduce the number of elements used to construct the propagation operator in the LSM scheme.
We use this property in the imaging phase. For each frequency we threshold the highest 20% beam amplitudes, and use only the imaging kernels corresponding to these m m¢ , values to construct the propagation operator. For comparison, we plot the corresponding data at f = 10 Hz in figure 5(b). One can see that the two data sets agree with each other. The images obtained by the sparse data set are plotted in figure 6. In figure 6(a) we plot the image obtained by the beam-based migration algorithm, while figure 6(b) illustrates the image obtained with the proposed LSM method after 30 iterations. We have 688 source and receiver beams at each frequency for the full data set, which requires 473 344 imaging kernels and cross-correlation operations to form the image at each frequency. With the sparse data set, only 94 668 cross-correlation operations are needed to calculate the gradient and the data residual at each frequency. We also compare the image results at x = 4.6 km in figure 6(c). In figure 6(d) we subtract the final reflectivity models obtained by the full and sparse data sets to compare between the results. For clarity, we normalize the results to the values of the reflectivity model in figure 1(b). As can be seen, even though we use only a sparse data set, the iterative scheme enables us to achieve a good result that is comparable with that obtained with the full data set, while reducing the computational complexity at each iteration.
As above, we compare the estimated and input data in figure 7 for a source located at x = 4.6 km. For clarity all plots are scaled with the same colorbar. Note that although we used a sparse data set the final result matches the input data.

Imaging with a full wavefield synthetic data
After exploring the algorithm via a Born data set, we explore the LSM algorithm via a full wavefield synthetic data set. The data set is generated using the time domain rapid expansion method [58]. The physical configuration is the same as the Born data set.
There are 460 shots, with 460 receivers for each shot. Shot-profile methods require solving and storing 460 GFs at each frequency. Here, as above, we use only 92 GFs to construct the beams.
The image results are plotted in figure 8. As can be seen, the proposed LSM algorithm substantially improves the subsurface image resolution, image amplitudes are balanced from the shallow to the deep part of the model. The dipping events at the bottom of the model, which are not visible in the regular migration image, can be clearly identified in the LSM result. Note that the full wave data doesn't fully obey the Born approximation. A different propagation operator is used to form the full wave data than for imaging. For this reason we see less accurate results than the image in figure 2(c).
To increase the accuracy of the full wavefield synthetic data LSM image, we use 460 GFs to construct the beam propagators. The image is plotted in figure 9(a). For comparison we also added the data obtained from a source located at x = 4.6 km. As can be seen, the two data sets agree.

Discussion
We introduced here a data domain LSM algorithm using beam propagators. The proposed approach is based on the migration operation of [46]. With this migration algorithm, we use WFT frame expansions to transform the LSM problem directly to the beam domain. We establish the data transformation (Born modeling) and migration (image reconstruction) operators and show that they are conjugate operators. This property allows us to use a gradient scheme to solve the LSM problem.
Data domain beam-based LSM algorithms are iterative schemes that may be implemented in the timedomain or FD. In a traditional shot-profile time-domain LSM formulation, the adjoint-state method is used to compute the gradient. These time-domain implementations require solving the wave equation twice for every shot-one for the forward and one for the backward propagation. Denoting by n as the number of iterations, n s the number of shots, the number of wave equations that needs to be solved is 2 × n s × n.
The current implementation is formulated using FD. In FD LSM implementations, we may calculate the gradient in two ways. The first is to use the adjoint-state method. In this approach, we multiply the inverse of the left hand side of Helmoltz equation to the right hand side twice for each shot over the entire frequency range. When the number of shots and iterations is large, the multiplication of the inverse of the left hand side to the  right hand side can becomes time-consuming. An alternative method is to use the scattering integral approach. There, we save GFs for all source and receiver in the computers memory or disk. Constructing the forward and backward wavefields repeatedly at each iteration is avoided. However, this implementation is storageconsuming and thus limited, especially when the number of sources and receivers is large, storing the GFs for every shot and receiver becomes prohibitive.
In the proposed method, we only need to compute a certain small number of GFs even with large numbers of shots and receivers. The GFs are used to numerically construct the beam propagators (see equation (A.5)). The number of GFs is determined to obtain a sufficient sampling interval over the surface. This sampling interval is generally determined by the Nyquist sampling interval (Δ x < λ/2). Note that this number only gives a bound to the number of GFs that have to be used. As shown in the numerical examples, we can use fewer GFs to obtain the image. The number of GFs is relatively small. Note that the same GFs are used to represent the source and receiver beams. We use here the same GFs throughout iterations. Thus there is no need to solve the wave equation or construct the forward and backward wavefield in each iteration. These points greatly reduce the computational cost of the LSM implementation.
The proposed method is different from the conventional GB LSM algorithms discussed in the introduction. Most of these formulations are based on the localized sources approach [33,34]. These methods are based on expanding each point source GF as a superposition of GBs. In the localized sources GB approaches, one constructs the beam propagators using a paraxial wave equation near the beam axis. These are ray-based methods that require solving a Riccati equation along the axis. On the other hand, as solutions to the wave equation, beams are insensitive to caustics and reduce multipathing. Our approach is based on a full wavefield solution, which improves the accuracy of wave propagation, hence improve the accuracy for the reflectivity estimation.
The use of the WFT expansions allows to a-priori reduce the number of propagators that construct the propagation operator L in the LSM problem. With the WFT expansions, we decompose the shot gathers into a beam to beam scattering amplitudes (beam domain data). Only pairs of beams with high amplitudes have to be taken into account when constructing L. This filtering step is performed as a pre-processing phase before the iterative scheme, and a priori reduces the number of beam cross-correlation operations used to form the image. Comparing to localized source GB methods, the steps described above reduce the computational complexity of the iterative scheme.
The UWB-PS-BS utilizes an overcomplete, frequency independent, beam lattice. This use allows us to transform the theory directly to the time-domain [59], while maintaining the physical properties of the beam amplitudes. On the other hand, due to the overcompleteness of the frame, the UWB-PS-BS method suffers from low overcompleteness at low frequencies (see [51] and also comment after (A.3)). This problem arises mainly in applications with a very wide frequency range (see [60]). Note that similar results for the oversampling criteria of the beams over the surface obtained in [33]. Thus the same overcompleteness arises in formulations based on the latter theory. To overcome this problem we use here a sparse data set, and reduce the number of imaging kernels used for imaging (see section 4.3.3). Alternatively, one may reduce the overcompleteness at low frequencies by using the multiband UWB-PS-BS algorithm introduced in [60]. In this approach, the frequency band is divided into sub-bands, and the UWB-PS-BS algorithm is implemented in each sub-band. As a result, we obtain fewer beam elements at the lower frequency bands while keeping the overcompleteness in all the bands below a given level.

Conclusions
We introduced here a new approach of beam-based data domain LSM using the UWB-PS-BS method. This method is based on WFT frame expansion of a surface field. The beams are used to decompose the surface field and then to propagate it. In the current implementation, we use this method to reduce the computational complexity of the data-domain LSM problem.
In the proposed method, the beams are calculated using the medium GFs. We use here a relatively small number of GFs to construct the beams, which can be calculated only once and store at the computer memory. Thus, no wave propagation needs to be performed at each iteration. Compared to conventional shot domain LSM implementations, we reduce the computational complexity of the problem. For a given subsurface point, only beams passing near it need to be considered. Our approach also allows us to further reduce the computational complexity of the LSM problem, in contrary to the localized sources beam-based LSM formulations. The beam amplitudes correspond to scattering events over the subsurface. Thus, we may filter out pairs of beams with relatively low magnitudes, which reduces the data size and the number of wave fields used to reconstruct the subsurface, while obtaining a good agreement with the true model. Since the beam amplitudes are calculated before the migration phase, this approach a-priori reduces the computational complexity of the LSM implementation. Finally, since the beam propagators are calculated via GFs obtained by solving the twoway wave equation, and not solving ray-based equation, we also enhance the image accuracy comparing with ray-based LSM methods. These unique properties are utilized here to improve the computational efficiency for the LSM problem.

Data availability statement
The data generated and/or analysed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request.
Appendix A. The ultra-wide-band phase-space beam summation method In this work we are apply a beam approach for data-domain LSM, utilizing the UWB-PS-BS method. The key features of the method are summarized below. Further details of the method, with application to seismic wave propagation are given in [46]. with y x ( ) is localized window function, typically obtained by multiplying 2 1D window functions. The frame elements are localized around the phase-space grid of points, defined by the 4-index μ via (see also [46], figure 1 . A . 2 m n¯( which simplifies the calculations. This choice also leads to stable and localized coefficients ( [51]). Equation (A.4) represents the radiated field as superposition of beam propagators, emanate from the points x = x m , and propagate along a direction defined by the initial ray parameter p = p n (see [46], figure 1), and decay away from it. The direction of propagation also known as the beam axis. The beam amplitudes are obtained by projecting the aperture field distribution on the dual frame functions. Thus, they may be interpreted as extracting the local radiation spectrum of u 0 around the point x m in the direction p n . The representation in equation (A.4) is a priori sparse both in the propagators and amplitudes. Giving a point over the subsurface only the beams which pass near this point have to take into account in the summation. As discussed after equation (A.6), m â extracts the local spectrum of u x 0 ( ). Thus, we may filter out beams with low amplitudes, which reduces the number of beams that need to be taken into account in the summation.
A.2. Special case: Iso-diffracting Gaussian windows As was shown in [32,51], a favorable choice for the function is the use iso-diffracting Gaussian windows. These are windows where their width is scaled with frequency by the frequency independent parameter b. These windows are given by For this case the dual window function is given by