Matthieu Gallet

Ammar Mian

Senior Member, IEEE Guillaume Ginolhac

Senior Member, IEEE Esa Ollila

Nickolas Stelzenmuller

New Robust Sparse Convolutional Coding Inversion Algorithm for Ground Penetrating Radar Images

Keywords: Ground Penetrating Radar, Sparse Inversion, Convolutive Model, Robust methods

In this paper, we propose two algorithms to enhance the interpretability of the hyperbola in B-scans obtained with a Ground Penetrating Radar (GPR). These hyperbolas are the responses of buried objects or cavities. To correctly detect and classify them, a denoising is typically necessary for GPR images as the signal-to-noise ratio is low, and the various interfaces naturally present in the earth have a strong response. Both algorithms are based on a sparse convolutional coding model plus a low rank component. It is solved through an Alternating Direction Method of Multipliers (ADMM) framework. In order to take into account the presence of outliers and the artifacts caused by the acquisition, the second algorithm is based on the Huber norm instead of the classic L2-norm. These algorithms are tested on a real dataset labeled by geophysicists. The results show the denoising efficiency of this approach, and in particular the robustness of the second algorithm.

I. INTRODUCTION

Ground Penetrating Radar (GPR) consists in transmitting an electromagnetic wave in the ground [START_REF] David | Ground Penetrating Radar[END_REF]- [START_REF] Benedetto | An overview of ground-penetrating radar signal processing techniques for road inspections[END_REF]. This wave can be reflected by buried objects or cavities to be then received by the radar. The linear displacement of the radar antenna during a measurement allows the system to create a 2D image which is called a radargram or B-scan. Since the radar is moving and several acquisitions are made, the response of an object has a form of hyperbola. Material interfaces, such as air-earth, rocksoil, among other naturally occurring layers in the earth also reflect the radar signal. These refelections are visibile in the B-scan, often appearing as quasi-horizontal lines. An example of B-scan with one main object, the response of the ground and the response of a layer is shown in Figure 1.

A typical objective in the analysis of a GPR image is to classify the underground objects. The main information about the object is contained in the shape of the hyperbola. For this task, several classification methods [START_REF] Besaw | Deep convolutional neural networks for classifying GPR B-scans[END_REF]- [START_REF] Yelf | Classification system for GPR parameters[END_REF] have been proposed that are based on supervised approaches like deep neural networks. These methods require a relatively large amount of labeled data to achieve an acceptable classification accuracy. However, the accurate labeling of B-scans is often difficult, even for trained users of GPR systems, due to the low Signal-to-Noise Ratio (SNR) of the objects of interest. Moreover, the responses of the various interfaces naturally present complicate the interpretation of the shape of the Matthieu Gallet, Ammar Mian and Guillaume Ginolhac are with LISTIC (EA3703), University Savoie Mont-Blanc, FRANCE. Esa Ollila is with Aalto University, Finland. Nickolas Stelzenmuller is with Geolithe, FRANCE. This work has been done thanks to the facilities offered by the Univ. Savoie Mont Blanc -CNRS/IN2P3 MUST computing center.

hyperbola. The improvement of the quality of these images would allow a more accurate interpretation and classification of objects found in GPR images.

One strategy to improve image quality involves increasing signal diversity in order to better distinguish the contribution of the various sources. For example, a stepped frequency GPR is used in [START_REF] Gurbuz | A compressive sensing data acquisition and imaging method for stepped frequency GPRs[END_REF] which improves the SNR and increases the effective depth of the measurement. Higher SNR has also been reported from measurements using novel GPR systems such as those using polarization diversity (like in Polarimetric SAR) or a Multiple Input Multiple Ouput radar [START_REF] Zeng | Improving target detection accuracy based on multipolarization MIMO GPR[END_REF]. The later option greatly increases the resolution of the B-scans. Unfortunately these novel systems are often impractical given cost, weight, and size constraints. For example, drone-based measurements are increasingly common, and impose strict size and weight constraints on GPR systems. To improve the quality of the images, it is then necessary to resort to signal or image processing techniques. For example, Kalman filters are often used signal processing technique for detection/localization of buried landmines as in [START_REF] Zoubir | Signal processing techniques for landmine detection using impulse ground penetrating radar[END_REF]. An improvement of this method using a particle filter is proposed in [START_REF] Ng | Particle filtering based approach for landmine detection using Ground Penetrating Radar[END_REF]. A landmine detection method based on hidden Markov model has been proposed in [START_REF] Manandhar | Multipleinstance hidden markov model for gpr-based landmine detection[END_REF]. A main problem in B-scan images is the responses of the ground and the layers which are removed by a correlation-based approach in [START_REF] Kovalenko | A novel clutter suppression algorithm for landmine detection with GPR[END_REF] or a PCA-based approach [START_REF] Li | Gpr signal denoising and target extraction with the ceemd method[END_REF]. Statistical approaches have been also proposed in the literature. A seminal work [START_REF] Brunzell | Detection of shallowly buried objects using impulse radar[END_REF] developed an energy detector. A more elaborated detector based on a GLRT is given in [START_REF] Ho | A linear prediction land mine detection algorithm for hand held ground penetrating radar[END_REF]. This method takes into account the correlation of the noise but it works column by column and assumes that the noise is Gaussian. A detector working on the full image and considering non-Gaussian noise has been developed in [START_REF] Hoarau | Robust adaptive detection of buried pipes using GPR[END_REF]. This method improves the object detection and allows to remove clutter. These methods have been shown to be effective in improving the SNR and detection performance, but they are often targeted to a specific type of object, such as landmines or pipes and do not allow a better classification since the information regarding the properties of the hyperbolas is lost.

In order to develop approaches that are applicable to a wide range of objects, we investigate inversion algorithms. In the GPR domain, most inversion methods are based on a precise physical model by computing for example the Green's function. These approaches consist then in applying a full waveform inversion algorithm which is based on a simple L 2 norm and a Thikonov regularization [START_REF] Meincke | Linear gpr inversion for lossy soil and a planar air-soil interface[END_REF]- [START_REF] Wang | Multiparameter full-waveform inversion of 3-d on-ground gpr with a modified total variation regularization scheme[END_REF] or an interpolation [START_REF] Wang | Efficient implementation of gpr data inversion in case of spatially varying antenna polarizations[END_REF]. Approaches based on compressive sensing [START_REF] Ambrosanio | A compressive-sensing-based approach for the detection and characterization of buried objects[END_REF] have also been proposed. However, all of these models require a complex computation of the direct model and do not take into account layer suppression. For the latter, a spatial-averaging-based pre-processing step or PCA-based approaches as in [START_REF] Li | Gpr signal denoising and target extraction with the ceemd method[END_REF] are often used before inversion. Methods using deep learning have also been developed in [START_REF] He | Gpr image denoising with nsst-unet and an improved bm3d[END_REF], [START_REF] Dai | Dmrf-unet: A two-stage deep learning scheme for gpr data inversion under heterogeneous soil conditions[END_REF] [29]- [START_REF] Liu | Gprinvnet: Deep learning-based groundpenetrating radar data inversion for tunnel linings[END_REF] but these types of approaches require a learning step which is non-robust to datasets of different modalities (frequencies of the radar, elevation, ...). In order to use the inversion approach but with a simple model for the dictionary, the sparse coding representation is a well-suited model to correct the shape of the hyperbola in B-scan images [START_REF] Shao | Sparse representation of gpr traces with application to signal classification[END_REF], [START_REF] Hing | Gpr signal classification with low-rank and convolutional sparse coding representation[END_REF]. Indeed, the dictionary could be easily constructed from a family of hyperbolas and even estimated from the data, as in [START_REF] Giovanneschi | Dictionary learning for adaptive gpr landmine classification[END_REF]. But as for the previous algorithms, these methods do not take into account the responses of the ground and the layers. A seminal work combining the suppression of the clutter and the enhancement of the hyperbola was proposed in [START_REF] Terrasse | Sparse decomposition of the GPR useful signal from hyperbola dictionary[END_REF]. The method is based on the Robust Principal Component Analysis (R-PCA) initially proposed in [START_REF] Candès | Robust principal component analysis?[END_REF] and extended to the convolutional model in [START_REF] Wohlberg | Convolutional sparse representations with gradient penalties[END_REF]. The proposed algorithm is effective because it removes the contribution of the ground and the layers while the hyperbola is enhanced, allowing better interpretation of the B-scan image. But some drawbacks remain such as the construction of the dictionary which does not explicitly use the physics of the wave propagation. This can lead to physically unrealistic hyperbola which can reduce the performance of the final algorithm.

In this paper, we propose a realistic way to build our dictionary with limited size while keeping a good representation of all possible hyperbolas. One drawback of algorithm proposed in [START_REF] Terrasse | Sparse decomposition of the GPR useful signal from hyperbola dictionary[END_REF] is that it is not very robust to noise and other artifacts caused by the acquisition of the data (such as outliers, disalignment, ...). To take such effects into account, we propose a new inversion approach based on the Huber norm [START_REF] Huber | Robust Statistics[END_REF]. This new operator is well known to be more robust than the classical norm in optimization algorithms with regards to outliers and high noise levels. To the best of our knowledge, this is first time that this operator is used in a sparse convolutional coding estimation with the presence of a low rank contribution. In this paper, we develop an Alternating Direction Method of Multipliers (ADMM) framework to solve the optimization problem using L 2 -norm and robust Huber norm. Both algorithms are presented and differences in their steps are underlined. In particular, we discuss the choice of the different hyper-parameters of the algorithms as well as how to improve its convergence.

Finally, both algorithms are tested using a dataset provided by the Geolithe company. This dataset consists of several Bscan images composed of different objects of interest. The main advantage of this dataset is that it is labeled by geophysicists from the company, who provided the corresponding mask images locating and labeling the different hyperbolas. We use these data to demonstrate the effectiveness of the robust approach in removing noise and acquisition artifacts from the images.

We summarize the main contributions of the paper in the following: i) we propose a realistic way to build our dictionary with limited size while keeping a good representation of all possible hyperbolas. This dictionary is based on physical

. ∥ • ∥ 2 is the L 2 -norm, ∥ • ∥ 1 is the L 1 -norm, ∥ • ∥ * is the nuclear
norm and ∥ • ∥ H is the Huber norm defined at [START_REF] Ng | Particle filtering based approach for landmine detection using Ground Penetrating Radar[END_REF] . * (resp. ⊛) is the convolution operator for functions (resp. matrices) . ⊙ denotes the Hadamard product . for a vector x of dimension d, diag(x) is an operator that returns a diagonal matrix of size d × d. . for a matrix X, vec(X) yields its vectorization while unvec(•) allows to get back the matrix . F (•) is the Fast Fourier Transform of a matrix. F -1 (•) is the inverse transform . given a multivariate function f with scalar output, ∇xf denotes the gradient of f with regards to the elements of vector x . δ(x) is the Dirac function equal to 1 for x = 0 and equal to 0 everywhere else . I is the identity matrix and 1 d is the vector of dimension d with all elements equal to 1.

Inversion model:

. Nt, Nx denote the number of samples in dimensions t and x. K is the number of atoms in the inversion model .

Y, N, L, H k , C k , S k and M are matrices of size Nt × Nx . the notation {• k } is used to simplify {• k } k=1,...,K

Physics:

. c is the speed of light, ϵr stands for relative permittivity, σ stands for conductivity and ω is the RADAR frequency.

Optimization and ADMM framework:

. λ is a sparsity favoring hyperparameter . L and S k are the dual variables of L and S k . L is the augmented Lagrangian . ρ L and ρ S are the augmented Lagrangian penalty parameters regarding L and S . • (i) is the value of a parameter at iteration i of the global ADMM algorithm. . for matrices Y, H k , C k , S k , S and L, we define y, h k , c k , s k , s k and l by applying the operation vec(F (•)) . Sγ (•) is the proximal soft thresholding operator [START_REF] Dai | Dmrf-unet: A two-stage deep learning scheme for gpr data inversion under heterogeneous soil conditions[END_REF] . Tγ (•) is the proximal singular value thresholding operator [START_REF] Giovanneschi | Dictionary learning for adaptive gpr landmine classification[END_REF] . PH γ,δ (•) is the Huber proximal operator [START_REF] Wang | Efficient implementation of gpr data inversion in case of spatially varying antenna polarizations[END_REF] considerations and could be built for any dataset by using the RADAR parameters and a linescape on the permittivity and the radius of the object. ii) we develop a new robust inversion algorithm by replacing the L 2 -norm of [START_REF] Terrasse | Sparse decomposition of the GPR useful signal from hyperbola dictionary[END_REF] by robust Huber norm. The ADMM framework is used to solve the optimization problem with new updates for the low-rank contribution and the coefficient maps. Note that this approach is nonsupervised and do not need any training steps and labeled data. iii) The proposed algorithms shows a real improvement on a real dataset.

The article is structured as follows. Section II presents the data model and the procedure for building the dictionary dedicated to GPR images. In the Section III, the algorithm inspired by [START_REF] Terrasse | Sparse decomposition of the GPR useful signal from hyperbola dictionary[END_REF] is given. The main contribution of the paper is given in Section IV where the new optimization problem with Huber norm and the corresponding algorithm are described. Finally both algorithms are tested on a real dataset in the section V. Table I summarizes the notations used in throughout the paper.

II. DATA MODEL OF GPR

In this section we introduce a modeling of the measured GPR data by splitting the backscattering of buried objects from clutter and noise components. From this model, we obtain a formulation involving a sparse convolutional coding of probable targets backscattering. Finally, a physical modeling is proposed to construct a dictionary of possible target responses that will be subsequently used in section III.

A. Data model

A B-Scan radargram can be decomposed as follows:

y(x, t) = s(x, t) + cl(x, t) + n(x, t) (1)
where:

• n(x, t) is a noise component relative to the acquisition • s(x, t) is the backscattering of the different targets. Given the acquisition paradigm of GPR, a single target signal is hyperbola-shaped. In the presence of multiple targets, we have a sum of different hyperbolas [START_REF] Terrasse | Sparse decomposition of the GPR useful signal from hyperbola dictionary[END_REF]:

s(x, t) = K k=1 c k (x, t) * h k (x, t) (2)
where h k (x, t) is the image of a centered hyperbola and c k (x, t) the corresponding coefficient map. • cl(x, t) is a clutter component representing all the elements which cannot be represented by eq.( 2) such as the response of the different layers in the ground. Since the signal y(x, t) is sampled in both dimensions, we obtain a matrix which can be written in the form:

Y = K k=1 C k ⊛ H k + L + N, (3)
where the dimensions of the matrices are equal to those of the B-scan image, i.e. N x × N t and ⊛ is the matrix convolution product. The matrix L contains the responses of the clutter such as horizontal layers or acquisition artifacts. This matrix is assumed to be low-rank since those clutter signals are similar between themselves. The matrix N is mainly composed of the contributions of the different noises. Each H k represents the hyperbola signal of each individual target in the scene and C k is its corresponding mapping. In this scenario, C k is expected to be very sparse, i.e., it has a few values where each target is located. When the targets are not precisely known, this model allows to consider a large dictionary of possible target signatures {H k } k=1,...,K for which even fewer associated mappings {C k } k=1,...,K are non-zeros. We propose hereafter to construct such a dictionary of target signatures based upon a physical modeling of the interaction between a buried object and the emitted electromagnetic wave.

B. Construction of the dictionary {H

k } k=1,...,K
In [START_REF] Terrasse | Sparse decomposition of the GPR useful signal from hyperbola dictionary[END_REF], the dictionary is not completely related to the physical behavior of the buried objects, which leads to a dictionary with outlier elements. Figure 2 displays two elements of the dictionary developed in [START_REF] Terrasse | Sparse decomposition of the GPR useful signal from hyperbola dictionary[END_REF]. We notice that the element shown in the left image is a hyperbola with a physically-plausible shape, but the element on the right panel has an unrealistic shape that will never be present in real GPR images. To overcome this issue, we propose a new way to build our dictionary based on the physics of the GPR signal. Fig. 2: Two elements of the dictionary developed in [START_REF] Terrasse | Sparse decomposition of the GPR useful signal from hyperbola dictionary[END_REF]: the left one has a correct shape whereas the right one is an outlier.

First, we recall the expression of the Ricker wavelet which is the transmitted signal by the GPR [START_REF] Ramadhan | The effect ricker wavelet of duty cycle adjustment on gpr detection result[END_REF], [START_REF] Wang | Frequencies of the Ricker wavelet[END_REF]:

r(t|ω max ) = (1 - ω 2 max t 2 2 ) exp(- ω 2 max t 2 4 ), (4)
where ω max is the maximal pulsation (in rad • s -1 ) of the RADAR. By using the Ricker formulation (4), we are now able to construct the atom of our dictionary by using the following expression:

h(x, t) = A(x, t) • (r(t) * δ dirac (t -g(x))), (5)
where A(x, t) is an attenuation value which ensures the finite dimension of the hyperbola and g(x) is the hyperbola defined by:

g(x|x 0 , t 0 , a, p) = a p 2 + (x -x 0 ) 2 + t 0 -pa (6)
where (x 0 , t 0 ) is the top of the hyperbola. It is useful to connect (x, t) to the true dimensions of a GPR image and we propose the following interpretation:

• the number of pixels along the x axis and the distance covered by the radar in meters,

• the number of pixels along the t axis and the propagation time (one way) of the electromagnetic wave.

Thus the dimensions x and t become physical (respectively meter and second). We denote the sampling intervals in both dimensions (∆ x , ∆ t ). The choice of x 0 and t 0 can be made to maximize the size of the hyperbola in the matrix, and we choose x 0 = max(x)/2, t 0 = max(t)/4. A precise position does not matter since the convolution product allows us to correctly place the hyperbola in the final image. An example of Ricker wavelet and the corresponding attenuation matrix A(x, t) are shown in Fig. 3. To obtain all the atoms of the dictionary, we have to compute several hyperbola g(x) of ( 6) by choosing different values of p and a. These parameters are linked to the physical properties of the materials crossed by the transmitted signal of the GPR and could be written as follows:

p = t 0 v 2 + R a = 2 v , (7)
where:

• R is the radius of the target illuminated by the RADAR,

• v is the wave speed and is expressed w.r.t to physical parameters of the traversed material:

v = c µrϵr 2 1 + ( cr ω ) 2 + 1 1/2 , (8)
with c being the speed light, ω the RADAR pulsation, µ r the relative permeability of the material, ϵ r the relative permittivity of the material and c r the conductivity of the material. Since the materials are mostly non-magnetic and nonconductive, µ r and c r will be the same for all experiments (µ r = 1 and c r = 0) 1 which leads to a simplification of (8):

v = c/ √ ϵ r (9)
Therefore, we are going to choose different relative permittivity ϵ r besides R to build our dictionary. We propose two choices to find an interval of relative permittivity: either to have several relative permittivity ϵ r depending of the scene illuminated by the RADAR or to compute an effective permittivity ϵ that takes into account a possible target of relative 1 Because of the small conductivity to the target materials such as wood, silica and plastic permittivity ϵ i embedded in a material of relative permittivity ϵ m . This latter is obtained by the Maxwell-Garnett model:

ϵ = ϵ m 2δ i (ϵ i -ϵ m ) + ϵ i + 2ϵ m 2ϵ m + ϵ i -δ i (ϵ i -ϵ m ) (10)
where δ i is the fraction of the inclusion of the target in the material.

To conclude, we resume the way to easily build a dictionary. First, we simply take some parameters of the RADAR like f max , ∆ x and ∆ t . Finally, by varying ϵ (or simply ϵ r ), and R, we obtain a collection of hyperbolas giving us the dictionary {H k } 1,K . The interval of ϵ and R can be large but it is important to avoid having a dictionary that is too large, especially for computational complexity and memory issues. An example of four atoms are shown in Fig. [START_REF] Besaw | Deep convolutional neural networks for classifying GPR B-scans[END_REF]. We notice that the convolution of g(x) by the Ricker wavelet function r(t) allows to obtain the triple polarity (-/+/-) which is a main characteristic of GPR images as shown in Fig. 1.

III. CLASSICAL CONVOLUTIVE INVERSION

In Section II we proposed a model of the GPR signal involving an useful signal, clutter elements and noise. The useful signal is described using a sparse coding of possible target signatures in dictionary related to physical parameters of the scene. We consider the problem of estimating the mapping {C k } k=1,...,K and clutter matrix L allowing us to study the composition of a GPR image in terms of buried targets and layer composition.

A. Posing the inversion problem

In order to invert the GPR model, the following optimization problem function is proposed in [START_REF] Terrasse | Sparse decomposition of the GPR useful signal from hyperbola dictionary[END_REF]:

({C k }, L) = argmin {C k },L ∥Y - K k=1 C k ⊛ H k -L∥ 2 2 data fidelity + λ K k=1 ∥C k ∥ 1 sparsity promotion + κ∥L∥ * rank minimization , (11)
where λ and κ are penalty parameters. This cost function is designed to take into account three effects:

• data fidelity term measures the discrepancy of the model to the measured image signal,

• a sparsity promoting term by using a L 1 -norm on each coefficient map C k ,

• a rank minimization term by using a nuclear-norm on the clutter matrix. While not strictly minimizing the rank, it is known to promote low-rank solutions while retaining convexity [START_REF] Recht | Necessary and sufficient conditions for success of the nuclear norm heuristic for rank minimization[END_REF].

B. Solving using an ADMM formulation

In order to solve the sparse deconvolution problem (11) with additional variables, several approaches can be leveraged as discussed in [START_REF] Garcia | Convolutional dictionary learning: A comparative review and new algorithms[END_REF]. The ADMM framework in particular is reportedly the most efficient approach compared to others such as FISTA [43, Sec. IV.B]. The idea behind the ADMM is to split the optimization problem into sub-problems involving only one unknown at a time [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF], making it easier to handle problems having as many unknowns as in GPR inversion problem.

In order to apply this framework, we introduce auxiliary variables {S k } k=1,...,K2 . This permits the same formulation in the augmented Lagrangian of [START_REF] Zeng | Improving target detection accuracy based on multipolarization MIMO GPR[END_REF]. In this case we discard κ since the interplay between λ and ρ S parameters can be used to favor rank-minimization:

   minimize ∥L∥ * + λ K k=1 ∥S k ∥ 1 s.t Y = K k=1 C k ⊛ H k + L ∀k, S k = C k (12)
The augmented Lagrangian of this problem is as follows:

L {C k }, {S k }, L, { S k }, L = ρ S 2 ∥ K k=1 C k -S k + S k ∥ 2 2 + ∥L∥ * + ρ L 2 ∥Y - K k=1 C k ⊛ H k -L + L∥ 2 2 + λ K k=1 ∥S k ∥ 1 - ρ S 2 ∥ K k=1 S k ∥ 2 2 - ρ L 2 ∥ L∥ 2 2 , (13)
where { S k } k=1,...,K and L are the dual variables of {S k } k=1,...,K and L, and, ρ S and ρ L are Lagrangian penalty parameters 3 . The ADMM framework solves separably each optimization problem for one variable (by assuming the remaining constants) and iterates until convergence. It is applied here by solving the augmented Lagrangian iteratively for variables {C k }, then {S k } and L (by assuming the remaining ones constants). At each iteration, we also update the dual variables thanks to the residual of the equality constraints. The algorithm to solve [START_REF] Zeng | Improving target detection accuracy based on multipolarization MIMO GPR[END_REF] has been proposed in [START_REF] Terrasse | Sparse decomposition of the GPR useful signal from hyperbola dictionary[END_REF]. A summary is given in A with the notations used in this paper.

In next section we propose an improvement of this approach by replacing the L 2 -norm by a robust norm to obtain better performance of denoising.

IV. ROBUST CONVOLUTIVE INVERSION

In previous section, we have solved the optimization problem [START_REF] Yelf | Classification system for GPR parameters[END_REF], where the attach to the data is performed by exploiting a L 2 -norm on the residual between the reconstructed image and the original one. This approach is efficient in the absence of outliers or when the original B-scan image is subjected to low-level of noise. But when some outliers, such as discontinuities in the radargram due to poor pre-processing, are present or when the measurement is very noisy this method performs poorly. Indeed, in those scenarios, the use of the L 2 -norm will favor solutions that replicate those outliers and noise components. For such cases, it is possible to mitigate these effects by employing a robust cost function that finds a compromise between outlier rejection and closeness to the original problem. In this section, we propose to employ the well-known Huber norm [START_REF] Huber | Robust Statistics[END_REF] used in linear regression tasks. In particular, we consider the modified Huber function proposed in [START_REF] Ollila | Multichannel sparse recovery of complex-valued signals using huber's criterion[END_REF] to handle the case of complex-valued data in the Fourier domain when handling the sub-problem to find {C k }:

H δ (x) = |x| 2 , |x| ≤ δ 2δ|x| -δ 2 , |x| > δ , ( 14
)
where δ is a parameter allowing to perform a compromise between a cost close to the L 2 -norm (high value of δ) and a cost reducing the effect of outliers (lower value of δ). In practice, it can be estimated using the empirical distribution of the image pixel values and selecting a quantile. Based on this function, we define the Huber norm of a vector x ∈ C p :

∥x∥ 2 H = p i=1 H δ (x i ). (15)
The derivative function of (14) will also be needed and is equal to:

Ψ δ (x) = x , |x| ≤ δ δ sign(x) , |x| > δ , (16)
where sign(x) is equal to 0 when x = 0 and equal to x/|x| otherwise.

Using Huber's norm we have the following optimization problem:

argmin {C k },L ∥Y - K k=1 C k ⊛ H k -L∥ 2 H + λ K k=1 ∥C k ∥ 1 + κ∥L∥ * , (17)
that we rewrite into the form given below in order to exploit the ADMM framework:

       minimize ∥Y - K k=1 C k ⊛ H k -L∥ 2 H + κ∥M∥ * + λ K k ∥S k ∥ 1 subject to M = L ∀k, S k = C k , ( 18
)
where M is an auxiliary variable introduced to decouple the minimization of the rank of the clutter from the minimization of the Huber-norm. This implies that a new step for the variable M needs to be introduced to the ADMM iterations.

The augmented Lagrangian of this new problem is as follows:

L {C k }, {S k }, L, M, { S k }, L = κ∥M∥ * + λ K k=1 ∥S k ∥ 1 + ∥Y - K k=1 C k ⊛ H k -L∥ 2 H + ρ S 2 ∥ K k=1 C k -S k + S k ∥ 2 2 - ρ S 2 ∥ K k=1 S k ∥ 2 2 + ρ L 2 ∥M -L + L∥ 2 2 - ρ L 2 ∥ L∥ 2 2 . (19)
We describe below the derivation of the updates for each variable. Main contributions for this algorithm are in the updates of {C k } and L. A value of a variable at a given iteration i is denoted • (i) . A summary of the algorithm is presented in Algorithm 1.

A. Update {C k }

As this problem involves a heavy convolution, we can transpose this problem to the frequency domain with a Fast Fourier Transfrom (FFT) as in [START_REF] Wohlberg | Efficient convolutional sparse coding[END_REF]. We vectorize the resulting matrices. For matrices Y, H k , C k , S k and S, of size N t ×N x , we define vectors y, h k , c k , s k and s k of size N = N t N x by applying the operation vec(F(•)). Note that these vectors are now complex-valued due to the FFT operation.

We also define the following variables z k = s k + s k and x = yl, and rewrite minimization of (19) w.r.t {C k } into:

{c k } (i+1) = argmin {c k } ∥ K k=1 h k ⊙ c k -x∥ 2 H + ρ S 2 ∥ K k=1 c k -z k ∥ 2 2 ,
where ⊙ is the Hadamard product and ∥ • ∥ 2 is the L 2norm on complex vectors. Replacing the Hadamard product with a matrix product by concatenating the K different vectors to form the new

KN vectors c = [c T 1 , . . . , c T K ] T and z = [z T 1 , . . . , z T K ] T and defining a N × KN matrix H = [diag(h 1 ), . . . , diag(h K )],
Algorithm 1: Robust Inversion of GPR model with ADMM Data: Image Y ∈ R Nt,Nx , dictionary {H k ∈ R Nt,Nx } k=1,...,K Parameters: Sparsity parameter λ, Low-rank parameter κ, lagrangian penalties ρ S , ρ L , stopping criterion ξ Result:

Coefficients maps {C k ∈ R Nt,Nx } k=1,...,K , clutter matrix L ∈ R Nt,Nx Initialize {C (0) k } k=1,...,K , {S (0) k } k=1,...,K , L (0) , M (0) , { S (0)
k } k=1,...,K and L (0) with zeros Precompute Fourier elements y, {h k } as in IV-A i ← 0 while η (i) > ξ as in [START_REF] Maas | Using pattern recognition to automatically localize reflection hyperbolas in data from ground penetrating radar[END_REF] do Compute c

(i) k , s (i) k , s (i)
k and l (i) as in IV-A Update {C [START_REF] He | Gpr image denoising with nsst-unet and an improved bm3d[END_REF], M (i+1) with ( 29) and L (i+1) with ( 26) ∀k, S

(i+1) k } with (22), {S (i+1) k } with
(i+1) k = S (i) k + S (i+1) k -C (i+1) k L (i+1) = L (i) + K k=1 C (i+1) k ⊛ H k + L (i+1) -Y i ← i + 1
so that the optimization problem becomes:

c (i+1) = argmin c ∥Hc -x∥ 2 H + ρ S 2 ∥c -z∥ 2 2 , (20)
We have to compute the gradient of this equation w.r.t. c. For the second term, we have the same result as in the classical inversion (see [START_REF] Terrasse | Sparse decomposition of the GPR useful signal from hyperbola dictionary[END_REF] or A), but for the first term, the derivative is given by H † Ψ δ (Hcx) where the function Ψ δ is applied element-wise. We introduce the new vector r δ = [Ψ δ (Hcx)]. By using the structure of matrix H and using linear algebra identities, the total gradient is given by:

∇ c L = h † ⊙ (1 K ⊗ r δ ) + ρ S (c -z), (21)
where h = [h T 1 , . . . , h T K ] T . This equation does not admit an analytical solution and we use a simple gradient descent approach method to estimate c (i+1) :

c (i+1) 0 = c (i) , c (i+1) j+1 = c (i+1) j - ζ c j ∇ c (i+1) j L, (22)
where ζ c is the initial step-size of the gradient. While a stopping criterion can be used, experiments show that a satisfactory solution can be obtained with a few iterations of the gradient descent. Finally, {C

} is obtained by splitting c into the set of vectors {c (i+1) k } and by applying the operation F -1 (unvec(•)).

B. Update L

In this case, the optimization problem becomes:

L (i+1) = argmin L ∥Y- K k=1 C k ⊛H k -L∥ 2 H + ρ L 2 ∥M-L+ L∥ 2 2 ( 23
)
This problem is simplified by using the new variables

A = K k=1 C k ⊛ H k -Y and V = M + L: L (i+1) = argmin L ∥A + L∥ 2 H + ρ L 2 ∥V -L∥ 2 2 . (24)
For this optimization problem we can use the Huber proximal [START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF] defined by:

PH γ,δ (x) = x γ+1 if |x| < δ(γ + 1) x -δγ sign(x) otherwise . (25)
Then using proximal identities [START_REF] Parikh | Proximal algorithms[END_REF] we have the update:

L (i+1) = -A + PH ρ -1 L ,δ (V + A), (26)
where the proximal is applied element-wise.

C. Updates {S k } and M

In this case, we need to solve:

{S k } (i+1) =argmin {S k } λ K k=1 ∥S k ∥ 1 + ρ S 2 ∥ K k=1 S k -C k + S k ∥ 2 2 .
This problem is well known and the solution is given by [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]:

S (i+1) k = S λ/ρ S (C k -S k ) (27)
where S λ/ρ S is the soft thresholding operator defined for a real-valued vector u and scalar value γ by:

S γ (u) = sign(u) ⊙ max(0, |u| -γ), (28)
where the sign and max operations are applied element-wise. For M, the derivation is also very similar and leads to:

M (i+1) = T κ/ρ L (L -L). (29)

D. Algorithm tuning

Convergence and stopping criterion: The optimization problem is convex over all the variables and constraints and the convergence is ensured by the ADMM framework. Nevertheless, this convergence can be slow and can be improved by employing an over-relaxation technique with a relaxation term [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF], [START_REF] Eckstein | J. parallel alternating direction multiplier decomposition of convex programs[END_REF] α > 1 in the update C k .

To stop the algorithm, we use a distance between two iterations of the ADMM over their reconstructed image thanks to L and {C k } 1,K :

η (i+1) = ∥L (i+1) -L (i) + K k=1 C (i+1) k ⊛H (i+1) k -C (i) k ⊛H (i) k ∥ 2 .
(30) Choice of ρ S , ρ L , κ, ζ c and λ: The values of ρ S and ρ L are obtained by design of experiments techniques (see section V). These parameters can remain fixed throughout the iterations of the algorithm. It also possible to update them at each iteration by using the primal and dual residuals to improve convergence rate as discussed in [START_REF] Wohlberg | Convolutional sparse representations with gradient penalties[END_REF], [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. The low rank parameter κ is fixed to 1 because the associated Lagrangian ρ L allows to control the low rank character. The gradient step ζ c is found by a cross validation approach. Actually, all these parameters are mainly chosen to ensure the convergence of the ADMM and to do not need to be next changed even if the dataset is different.

The more important parameter is obviously the sparsity parameter λ as it controls the sparsity of the solution. It can be chosen as a compromise between high-response hyperbolas with associated low-noise and closeness to the original image. In practice, a cross-validation approach on a well-known scene can be used to determine it, or it can be adjusted with regards to geology expertise. In section V, we propose a crossvalidation approach but we also show that the performance of the proposed algorithm remains quasi constant over an interval of λ.

Solution without L: If the layers responses are removed by a pre-processing then one may omit the low-rank matrix L from the model. In this case, we can discard the corresponding update and the dual variable L from the equations. Moreover in the update of C k in Section IV-A, we only have x = y.

Computational complexity of the algorithm: The updates of the variables {S k }, { Sk } and L are relatively light since they are mainly based on thresholding methods and simple matrix operations (Hadamard product, addition,...). The update of {C k } is based on the computation of the gradient [START_REF] Hoarau | Robust adaptive detection of buried pipes using GPR[END_REF] which is based on a Hadamard product which is also light. This step is a little bit more complex since it needs a certain number of iterations. But this number can be chosen small (3 or 5 are enough to ensure the convergence of the ADMM) which limits the complexity cost of this step. The update of L in ( 26) is also light since the proximal is just a thresholding applied element-wise. Finally, the main complexity cost comes from the update of M in (29) which needs a SVD of a matrix of the size of the image. To have a better idea, we give an overview of the computation time. For the computation time, we have for 100 iterations with 3 iterations of gradient descent and a image of (N x , N y ) = (233, 128): 28.1s ± 1.39s. Compared to 24.8s ± 1.35s of the classical inversion algorithm of [START_REF] Terrasse | Sparse decomposition of the GPR useful signal from hyperbola dictionary[END_REF], the increase in computation time is reasonable. This study has been done with a laptop equipped with an 8-core i5-1135G7 @ 2.40GHz processor and 16Gb of RAM.

Another issue to investigate is the RAM memory needed for the updates. For this topic, the main issue comes from the gradient computation in [START_REF] Hoarau | Robust adaptive detection of buried pipes using GPR[END_REF] since it involves the matrix product with the matrix H which is the resultant of the diagonalization of the atoms of the dictionary (i.e. for k atoms of N × N pixels this gives us a matrix of size (k × N 2 , N 2 )). For example, if k = 40 and N = 128 in float32 it takes almost 37.3Go for a (655360 × 16384) matrix. However as only the values of the diagonals of the sub-blocks of the matrix are non null, the matrix calculation H † H becomes simply {h k } k=1,...,K ⊙ {h k } k=1,...,K and gives us a matrix of size (k × N 2 ). In a similar way, with k = 40 and N = 128 in float32 it takes only 2.36Mo.

V. NUMERICAL EXPERIMENTS

In this section we first describe experiments done on Bscans acquired in the field to showcase the benefits of the proposed convolutive inversion, as well as the use of the Huber norm. We also analyze the effect of hyperparameters tuning on the inversion performance. In the following, the classical algorithm developed in [START_REF] Terrasse | Sparse decomposition of the GPR useful signal from hyperbola dictionary[END_REF] and presented in the section III will be referred to as L 2 -GPR. Main details of the updates and the algorithm are given in appendix A. The proposed algorithm described in the section IV with Huber norm will be referred to as Hub-GPR. As a comparison between different methods, we compute Receiver Operator Curve (ROC) plots. In addition to the algorithms proposed in this paper, we have tested a simple SVD (to remove the spatially coherent noise) and a simple L 2 recovery (without the matrix L) combined with the SVD. Both methods will be referred by SVD-GPR and L 2 -SVD-GPR.

The code is available online: https://github.com/ Matthieu-Gallet/GPR-robust-inversion. The implementation of the algorithms was done under Python 3.6. We made use of the packages pyDOE 4 for the design of experiments in Section V-B and GPRPy5 for preprocessing the GPR images.

A. Dataset description and Dictionary parameters

The full dataset provided by Geolithe is composed of 1000 radargrams of (N x , N y ) = (4000, 800) pixels. These radargrams were acquired from a prepared site, in which various objects were buried at known locations and The objects vary in size ( 10 cm to 2 m) and material (wood, plastic, metal). The test site includes two 46m x 7m trenches, one composed of sand and the other of a mixture of sand, gravel and earth, including stones up to 10cm in diameter. GPR acquisitions were performed with two GSSI hyperstacking antennas. In addition to the GPR acquisitions performed by sliding the antennas along the ground( a typical utilization, with the highest SNR), acquisitions were also performed with the antennas held a fixed distances offset from the surface, up to 150 cm, in order to simulate a GPR acquisition by drone. The radargrams were annotated by geophysicists at Geolithe by tracing the signature of the object in the radargram. This pixel-wise annotation produces a labeled mask, in which all image pixels all labeled by class (e. g. object A, object B, no object, etc.) An example of a radargram with an annotation mask is shown in Fig. 6a. Although the extent of an object signature in the radargram is somewhat subjective (e. g. the bottom boundaries of the hyperboles signatures typically fade into the background noise of the radargram), the known position, type, and depth of the object allow us to use these annotations as a reliable ground truth. Many other characteristics are given for each image as illustrated in Table II. We build our dictionary using the following parameters provided by the GSII RADAR:

• ω max = 2πf max where f max is either 200 or 300 M Hz • ∆ x = 0.0101 m and ∆ t = 0.105•10 -9 s

We finally built a representative dictionary of K = 30 atoms by using for values of R ∈ {0.01, 0.1, 1} and ϵ ∈ {5, 6.46, 8.34, 10.77, 13.91, 17.97, 23.21, 29.97, 38.71, 50}.

B. Study of hyperparameters ρ S and ρ L

In order to check the effects of hyperparameters, we perform an experimental cross-validation study. We first used the inversion of the L 2 -GPR algorithm (the behavior is the same for Hub-GPR algorithm). In particular, we want to check the efficiency of the convolutional approach and the fit provided by our dictionary and verify if the choice of ρ S and ρ L is critical. We selected a representative radargram for this study. The original radargram is shown in Fig. 6a, and the mask image showing the true hyperbola is given in Fig. 6b. The acquisition of this radargram was made with an elevation of 75cm. We propose to study the robustness of our approach w.r.t. these hyperparameters (ρ S , ρ L ). To do this, we use the Latin hypercubes approach [START_REF] Franco | Planification d'expériences numériques en phase exploratoire pour la simulation des phénomènes complexes[END_REF], which allows us to carry out a design of experiment by judiciously sampling the parameter space. We transformed the uniform distribution of the Latin hypercubes samples to a log-uniform distribution in order to cover the large scale of both ρ S , ρ L (between 1 to 10 5 ). This approach allows minimizing the number of simulations by reducing the interactions between the parameters. The resulting image given in Fig. 5 allows us to conclude that to avoid a large error, it is sufficient to choose a value of (ρ S , ρ L ) smaller than 10 3 . To conclude, a narrow choice for (ρ S , ρ L ) is not needed to reach correct performance.

C. Performance of Hub-GPR algorithm

In this section, we investigate the potential benefits of using a Huber norm instead of the L 2 -norm in the optimization 0.00 problem. In order to have a better conclusion, we also tests different methods which are summarized in the following:

• Original image,

• SVD-GPR algorithm in which the first rank contribution is removed,

• L 2 -SVD-GPR algorithm where the sparse recovery is applied on the image preprocessed by the SVD (as before the first rank mode is removed),

• L 2 -GPR algorithm of section III,

• Hub-GPR algorithm of section IV. Figures 7a,7b, 7c and 7d respectively show the original radargram along with the results of the SVD-GPR, L 2 -SVD-GPR, L 2 -GPR and Hub-GPR algorithms. The parameter values are ρ S = 500 and ρ L = 250. The threshold for the Huber norm δ in ( 14) is set to 750. We choose the number of iterations of the gradient descent in [START_REF] Meincke | Linear gpr inversion for lossy soil and a planar air-soil interface[END_REF] equal to 5, and with a step size ζ c equal to 5 • 10 -4 . The number of iterations of the ADMM is 100. The sparsity parameter λ could vary between 0.2 and 0.8. The masks in red allow to clearly identify the objects of interest. Compared to the radargram in 6a, the result of all algorithms allows better detection of hyperbolas. The noise is also strongly reduced. This reduction is more visible with the Hub-GPR algorithm. Moreover, it seems that several hyperbolas are better detected in the image obtained by the latter.

For a quantitative study, we calculate ROC plots (true positive rate versus false positive rate parameterized by the detection threshold) 6 on the outputs of the considered methods.

The energy is next computed on each pixel for all resulting images before performing the detection step. For fairness, we have tested several values for λ for the three last methods and chosen the best.

The results for ROC plots are shown in Fig. 8. In all radargrams, we verify that the Hub-GPR algorithm gives the best result. The method seems to work well whatever the elevation (and thus the SNR) and the frequency of the RADAR. We also notice that the Hub-GPR algorithm improves the detection whatever the material of the object to be detected (wood, metal, plastic, ...).

D. Robustness to the sparsity parameter

In this section, we propose to study the robustness of the Hub-GPR algorithm to the choice of the sparsity parameter λ. For each image IRADAR, we compute the ROC plots for λ = 0.2, 0.4, 0.6, 0.8. The results are shown in Fig. 9. We first notice that the results keep close for all values of λ. Actually, the best value is around 0.4 and the performance decrease for 0.8.

E. Robustness to the dictionary

To do that, we made a simulation with different sizes of dictionaries:

• small: 18 atoms • medium: 36 atoms • big: 65 atoms We experiment on a crop 7 of the same image of the subsection V.C. The cropped image corresponds to the left part of the image split in two. The results shown in Fig. 10 showcase the effectiveness of the approach even with lower or bigger (when we are not sure about the objects present in the scene) number of atoms.

F. Robustness to noise of Hub-GPR algorithm

In this section, we measure the robustness of both approaches to noise. We will consider the original radargram acquired with GSSI software. The original as well as the labels mask, are shown in Fig. 11. For the rest of the experiments, we set ρ S to 500 and ρ L to 250. The threshold for Huber norm δ in ( 14) is fixed to 750. We choose the number of iterations of the gradient descent in [START_REF] Meincke | Linear gpr inversion for lossy soil and a planar air-soil interface[END_REF] equal to 5, and step size ζ c equal to 5 • 10 -4 . The total number of iterations of the algorithm is set to 20. This saves computation time, and experiments show that the solution is not so different with this reduced number of iterations. The sparsity parameter is equal to 0.4 as it is a value that gives good performance for the L 2 -GPR and Hub-GPR algorithms.

We consider two types of noise: one additive and one multiplicative. For both noises, we take a Gaussian noise N g of zero mean and variance between 0.01 and 10. For additive noise the noisy image is given by: I noisy = I + N g , and for multiplicative noise: I noisy = I + I * N g . The quantitative study is done on both the original un-processed and the preprocessed images. We compute several metrics: the classical MSE (Mean Square Error), the PSNR (Peak Signal to Noise Ratio) and the SSIM (Structural Similarity Index) [START_REF] Zhou Wang | Image quality assessment: from error visibility to structural similarity[END_REF]. These metrics with respect to the variance of the noise are given in Figs 12 and13. For all cases, the robust method gives better performance for all metrics, and shows its resilience to noise whether it is additive or multiplicative.

VI. CONCLUSION

We developed an inversion algorithm adapted to GPR images by using a convolutional model with a dictionary that is directly linked with to the physics of the GPR. We then generalized the approach by replacing the L2-norm in the problem of optimization by a Huber norm in order to cope with heavy-tailed noise, outliers and/or artifacts (misalignment, ...) often encountered in GPR images. Results of the new robust GPR inversion method showed the advantage of this new algorithm. Furthermore, our studies showed that the method is not overly sensitive to the choice of hyperparameters. All of the experiments were conducted on a real dataset provided by Geolithe, and labeled by geophysicists, which showcased the real world applicability.

Future works will focus on the construction of a method that can avoid gradient descent step and thus save computation times. Majoration-Minization framework [START_REF] Sun | Robust estimation of structured covariance matrix for heavy-tailed elliptical distributions[END_REF] will be investigated to solve this issue. Moreover, this convolutional approach will be extended to dictionary-learning approaches [START_REF] Garcia | Convolutional dictionary learning: A comparative review and new algorithms[END_REF] in order to be able to learn the dictionary directly from the data.

APPENDIX

The ADMM for the classical inversion problem ( 13) is summarized in Algorithm 2. We will next quickly describe the method to obtain the updates of {C k } k=1,...,K and L.

A. Update of {C k } k=1,...,K By using the notations in the frequency domain and the same transformations as IV-A, we obtain from ( 13) the new optimization problem:

c (i+1) = argmin c ρ L 2 ∥Hc -x∥ 2 2 + ρ S 2 ∥c -z∥ 2 2 , (31)
where x = yl + l and H, z, c are as defined in Section IV-A. Optimization problem [START_REF] Liu | Gprinvnet: Deep learning-based groundpenetrating radar data inversion for tunnel linings[END_REF] has been shown in [START_REF] Wohlberg | Efficient convolutional sparse coding[END_REF] to have a closed-form solution:

c (i+1) = (ρ L H † H + ρ S I) -1 (ρ L H † y + ρ S z), (32)
where (•) † is the conjugate transpose operator. For numerical stability, we can use the Sherman-Morrison formula instead of the inversion of (ρ L H † H + ρ S I). And finally, {C

} is obtained by splitting c into the set of vectors {c (i+1) k } and by applying the operation F -1 (unvec(•)).

As for the robust inversion, we obtain the coefficients values by reconstructing each c k and performing the operation F -1 (unvec(•)) on each coefficient map.

B. Update of L

We need to solve the following optimization problem:

L (i+1) = argmin L ∥L∥ * + ρ L 2 ∥ K k=1 C k ⊛ H k -Y + L + L∥ 2 2 ,
k + S (i+1) k -C (i+1) k L (i+1) = L (i) + K k=1 C (i+1) k ⊛ H k + L (i+1) -Y i ← i + 1
the solution of which is given by [START_REF] Candès | Robust principal component analysis?[END_REF]:

L (i+1) = T ρ -1 L (Y -L - K k=1 C k ⊛ H k ), (33)
where T ρ -1

L

is the proximal operator of thresholding the singular values [START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF]. For a matrix M with Singular Value Decomposition (SVD), M SVD = U diag(d)V T , we have:

T γ (M) = U diag(S γ (d))V T . (34)

Fig. 1 :

 1 Fig. 1: Example of a B-scan y(x, t) with t the delay of the wave and x the displacement of the radar.

Fig. 3 :

 3 Fig. 3: Main elements of the dictionary construction: Ricker wavelet (left) r(t|f max ) with f max = 350MHz (blue) and f max = 150MHz (orange) and attenuation matrix right A(x, t) (right).

Fig. 4 :

 4 Fig. 4: Main elements of the physical dictionary proposed in the paper (ω max = 2πf max = 2π • 350MHz ∼ 2.2e9 rad.s -1 , ∆ x = 0.0101 m and ∆ t = 0.105•10 -9 s). A: ϵ r = 1, R = 0.2, B: ϵ r = 1, R = 2, C: ϵ r = 100, R = 0.2 and D: ϵ r = 100, R = 2

Fig. 5 :

 5 Fig. 5: MSE w.r.t. the hyperparameters ρ S and ρ L . The dots display the sampling points used.

Fig. 6 :

 6 Fig. 6: Original radargram IRADAR 00H (6a) with the hyperbolas obtained by true objects labeled by Geolithe (6b). The red boxes show the locations of the labeled hyperbolas.

0.00 5 .Fig. 7 :

 57 Fig. 7: Results of the inversion algorithms for the IRADAR 00H radargram. ρ S = 500 and ρ L = 250. The threshold for the Huber norm δ is 750. The number of iterations of the gradient descent is 5, and with a step size ζ c equal to 5 • 10 -4 . The number of iterations of the ADMM is 500.

Fig. 8 :Fig. 9 :

 89 Fig.8: ROC plots for SVD-GPR, L 2 -SVD-GPR, L 2 -GPR and Hub-GPR algorithms on IRADAR 005, IRADAR 032, IRADAR 819 and IRADAR 00H radargrams. For each method, λ is chosen to give the best area. ρ S = 500 and ρ L = 250. The threshold for the Huber norm δ is 750. The number of iterations of the gradient descent is 5, and with a step size ζ c equal to 5 • 10 -4 . The number of iterations of the ADMM is 500.

Fig. 10 :Fig. 11 :

 1011 Fig. 10: Results of Inversion with Huber norm on a crop of image 00H

  Fig. MSE (top), PSNR (middle) and SSIM (bottom) w.r.t. noise variance with the original image and for the multiplicative noise. λ = 0.4 and the number of iterations is 20.

Fig. 13 :

 13 Fig. 13: MSE (top), PSNR (middle) and SSIM w.r.t. noise variance with the original image and for the additive noise. λ = 0.4 and the number of iterations is 20.

Algorithm 2 :

 2 Classical inversion of GPR model with ADMMData: Image Y ∈ R Nt,Nx , dictionary {H k ∈ R Nt,Nx } k=1,...,K Parameters: Sparsity λ, lagrangian penalties ρ S , ρ L , relaxation α, stopping criterion ξ Result: Coefficients maps {C k ∈ R Nt,Nx } k=1,...,K , clutter matrix L ∈ R Nt,Nx Initialize {C (0) k } k=1,...,K , {S (0) k } k=1,...,K , L (0) , { S (0)k } k=1,...,K and L (0) with zeros Precompute Fourier elements y, {h k } as in IV-A i ← 0 while η (i) > ξ as in[START_REF] Maas | Using pattern recognition to automatically localize reflection hyperbolas in data from ground penetrating radar[END_REF] do Compute c

TABLE I :

 I Notations and definitions used in the paperGeneral:. R is the set of real numbers and C is the set of complex numbers . lowercase (resp. uppercase) bold letter denote vectors (resp. matrices) . (•) T (resp. (•) † ) denotes the transpose (transpose conjugate) operator . | • | is the absolute value (resp. modulus) for a real-valued (resp. complex-valued) scalar

TABLE II :

 II Details of the images used in the study

Note that the term attached to the data has been replaced by the additional constraint

See[44, chap. 3] for details on the algorithm and definition of the terms used.

https://pythonhosted.org/pyDOE/

https://github.com/NSGeophysics/GPRPy

ROC are computed by using the scikit-learn method roc auc score with the masks and the image obtained by the studied algorithm.

To save time and test many sparsity options.


Laurent Nony

Jean-Marc Themlin

An introduction to the Discrete Fourier Transform and its applications in signal processing

This Jupyter notebook is meant to introduce the concepts of Discrete Fourier Transform (DFT) as a fundamental tool of signal processing. The theoretical foundations of the Fourier transform are introduced, however with a minimal mathematical formalism. The reader must have been introduced to the Fourier transform concept and must have some mathematical background before going through this document.

Here, the focus is rather brought on the technical implementation of the DFT, as well as on illustrations of the pratical use of the DFT by means of illustrative codes implemented in Octave 1 , a free programming langage that is compatible with Matlab.

The outline of this notebook is as follows: I-Foundations of the Fourier transform for continuous time signals II-The Discrete time Fourier Transform (DtDFT) III-The Discrete Fourier Transform (DFT) IV-Computation of the DFT V-DtFT & DFT: synopsis VI-Applications VII-DFT implementation using Octave: summary VIII-References

I-Foundations of the Fourier transform for continuous time signals I-1-Definition and conventions

A Fourier transform (FT) is a mathematical transform that decomposes functions depending on time, or space, into functions depending on temporal, or spatial, frequency. The Fourier transform of a function (here a signal) is a complex-valued function representing the complex sinusoids out of which the original function might be reconstructed. For each frequency, the magnitude of the complex value (modulus) represents the amplitude of a constituent complex sinusoid with that frequency, and the argument of the complex value represents that complex sinusoid's phase offset.

Here, the focus is brought to one-dimensional, time-dependent signals only, but the framework can be extended to 2D and space-dependent signals.

Let's assume a continuous time signal (i.e. analog signal) z(t), integrable in time (i.e. presenting no divergence), and its Fourier transform counterpart Z( f ) = FT{z(t)}, forming a Fourier transform pair:

z(t) ⇌ Z( f ) (1)
By definition, Z( f ) and z(t) are respectively derivable from each other according to :

Z( f ) := FT{z(t)} = +∞ -∞ z(t)e -j2π f t dt (2)
and, upon inverse FT operation:

z(t) := FT -1 {Z( f )} = +∞ -∞ Z( f )e +j2π f t d f (3)
Note: Another definition of the Fourier transform using the pulsation ω = 2π f instead of the frequency can be found in the text books:

Z(ω) := 1 √ 2π +∞ -∞ z(t)e -jωt dt (4)
But unless specified, we won't use this definition and we will stick to Z( f ).

I-2-Fourier series & Fourier transform

The FT formalism applies to any kind of signal, including periodic AND non-periodic signals.

For periodic signals however, Fourier series are a valuable tool, whose formalism is easier than the FT one. Besides, since one can always consider a single period T of any T-periodic signal (frequency f := 1/T), the Fourier series formalism obviously also applies to finite-duration signals of duration T.

Fourier series can represent T-periodic signals as the sum of sinusoids whose frequencies are integer multiples of the fundamental frequency f = 1/T of the periodic signal. The modulus and phase of the Fourier series coefficients represent the amplitude and phase of each harmonicallyrelated sinusoid. Hence the Fourier series coefficients spectrum of a T-periodic signal is discrete and its constituting frequencies are harmonics (i.e. integer multiples) of the fundamental frequency f = 1/T, whereas the Fourier transform spectrum of the same signal is continuous.

The table below summarizes the mathematical formalism for the Fourier series coefficients and for the Fourier transform.

II-The Discrete time Fourier Transform (DtFT)

In most applications, the FT is to be computed out of discrete time sampled signals. Therefore, it's natural to introduce the discrete time Fourier Transform (DtFT).

II-1-Time sampling

Let T s be the sampling period of the continuous time signal z(t), hence its corresponding sampling frequency:

f s := 1 T s (5)
The discrete time sampled signal is denoted z s (t). At that stage, we make no assumption on the duration of its time support, that may be considered as infinite, hence an infinite number of samples. z s (t) is derived out of the Dirac comb operator of temporal period T s , which is the sampling operator:

δ T s (t) = +∞ ∑ n=-∞ δ(t -nT s ), (6)
Therefore the continuous time t gets sampled according to:

t → t s [n] = tδ T s (t)| n = nT s , ∀n ∈ Z (7)
and then the discrete time sampled version of the continuous time signal is:

z s (t s ) = z(t)δ T s (t) (8)
Hence, the value of the n th sample, with n ∈ Z, of z s (t s ): z s (t s [n]) → z s [n].

The figure below summarizes the sampling process of a continuous time signal z(t) into its sampled version z s (t s ).

II-2-Formal expression of the DtFT

Let Z s ( f ) be the DtFT of the sampled signal z s (t s ), forming the FT pair:

z s (t s ) ⇌ Z s ( f ) (9)
The formal expression of Z s ( f ) can be derived as follows. As we are about to demonstrate it, a straightforward consequence of the temporal sampling process is that the DtFT is f s -periodic in the frequency space. Besides, it can be shown that the FT of a T s -periodic Dirac comb operator in time (T s = 1 f s ), is a f s -periodic Dirac comb operator in the frequency space, weighted by f s . The corresponding FT pair is:

δ T s (t) = δ 1/ f s (t) ⇌ f s δ f s ( f ) = 1 T s δ 1/T s ( f ) (10)
Thus, using the convolution theorem, the FT of the sampled signal z s (t s ) is:

Z s ( f ) := FT{z s (t s )} = FT{z(t)δ T s (t)} = FT{z(t)} * FT{δ T s (t)} = FT{z(t)} * f s δ f s ( f ) = Z( f ) * f s δ f s ( f ) (11)
The sampling process turns the FT into a DtFT, now defined by the formal expression:

Z( f ) → Z s ( f ) := Z( f ) * f s δ f s ( f ) (12) Because the δ f s Dirac comb replicates Z( f ) every f s , Z s ( f ) is f s -periodic.

II-3-Explicit expression of the DtFT

To derive a more explicit expression of Z s ( f ), we consider the discrete time expression of the FT of a continuous time signal (equ.2):

Z( f ) := +∞ -∞ z(t)e -j2π f t dt → Z( f ) = +∞ ∑ n=-∞ z s [n]e -j2π f nT s T s (13)
That is, according to equ.12:

Z s ( f ) := Z( f ) * f s δ f s ( f ) = ∑ +∞ n=-∞ z s [n]e -j2π f nT s T s * f s δ f s ( f ) = ∑ +∞ n=-∞ z s [n]e -j2π f nT s * δ f s ( f ) = ∑ +∞ n=-∞ z s [n]e -j2πn f fs * δ f s ( f ) (14)

III-The Discrete Fourier Transform (DFT)

The discrete Fourier Transform (DFT) of a generic discrete time signal z s (t s ) is an approximation of its DtFT (equ.14) obtained by :

• truncating the summation over n • discretizing the f parameter • truncating the f s -periodization

III-1-Truncating the summation: temporal windowing

When dealing with realistic discrete time sampled signals, it is impossible to handle an infinite number of samples. The discrete time sampled signal to be analyzed, z s (t s ), always consists of a finite number of samples. Let us assume that the sampling operation yields N samples of the discrete time signal. Hence the truncation of the summation over n in equ.14.

The N samples are obtained by windowing operation of z s (t s ), i.e. the multiplication of the signal by a windowing function w(t) to form the windowed sampled signal:

z sw (t s ) = z s (t s )w(t) (15)
w(t) may have a typical rectangular shape (but note that this definition is not exclusive), whose width sets the windowing duration T w , i.e. the duration of the time support of z sw (t s ). w(t) and its FT W( f ) form a FT pair defined as:

w(t) = 1, ∀t ∈ [0; T w ] 0, anywhere else ⇌ W( f ) = T w sin(πT w f ) πT w f := T w sinc(T w f ) (16)
By definition:

T w := NT s (17)
Now, the n th sample, with n ∈ N, n ∈ [0; N -1], of the windowed discrete time support t s is t s [n] := nT s , and the corresponding sample of

z sw (t s ) is z sw (t s [n]) → z sw [n].
The figure below completes the former figure about the sampling process of a continuous time signal z(t) into its sampled and windowed version z sw (t s ).

Important note: Although the last sample of the windowed discrete time support is written as (N -1)T s , the total duration of the signal is T w = NT s . In other words, each sample covers a duration of T s in time (cf. figure below, item 5).

III-2-Resuming the expression of the DtFT

The effect of the windowing operation is not trivial to the DtFT, which forces us to reconsider its expression (equ.14).

To that purpose, we resume the analysis of sections II-2 and II-3 out of the FT of the windowed continuous time signal Z w (t) := FT{z(t)w(t)}, thus: z(t)w(t) ⇌ Z w ( f ). It is also reminded that: w(t) ⇌ W( f ).

We then have:

Z w ( f ) := FT{z(t)w(t)} = Z( f ) * W( f ) → Z w ( f ) = N-1 ∑ n=0 z sw [n]e -j2π f nT s T s * W( f ) (18)
That is, by analogy with equ.12:

Z sw ( f ) := Z w ( f ) * f s δ f s ( f ) = ∑ N-1 n=0 z sw [n]e -j2π f nT s T s * W( f ) * f s δ f s ( f ) = ∑ N-1 n=0 z sw [n]e -j2π f nT s * W( f ) * δ f s ( f ) = ∑ N-1 n=0 z sw [n]e -j2πn f fs * W( f ) * δ f s ( f ) (19)
The previous equation constitutes the most complete expression of the DtFT of a continuous time signal z(t), whose windowed and discrete time version is z sw (t s ).

The figure below illustrates the influence of the temporal windowing, both in the time space and in the frequency space, by FT pairing. It is important to notice that we have taken a symetric interval of definition for the continuous time signal z(t).

III-3-Frequency sampling

III-3-a Sampled frequency Both the T s -sampling and the windowing in time force us to consider the sampling of the frequency, too. As seen before, the T s -sampling in time imposes a f s -periodization of the DtFT. In order to keep the same number of samples in the frequency domain, N samples will be used to sample the irreducible part (one period) of the DtFT. Therefore, the N samples of z s (t s ) resulting from the windowing operation feature a single f s period of its DtFT.

Hence a frequency sampling rate:

δ f = f s /N.
Therefore the continuous frequency gets sampled according to:

f → f δ f [k] := f δ δ f ( f )| k = kδ f = k f s N , ∀k ∈ Z (20)
To simplify the notations, we will write in the following:

f δ f [k] → f k .
Then the discrete frequency sampled version of the continuous frequency DtFT signal is:

Z sw ( f δ f [k]) = Z sw ( f k ) = Z sw ( f )δ δ f ( f )| k (21)
Hence, out of equ.19:

Z sw ( f k ) = N-1 ∑ n=0 z sw [n]e -j2πn f k fs * W( f k ) * δ f s ( f k ) × δ δ f ( f )| k ( 22
)
where the symbol × represents a product.

III-3-b Spectral resolution

The frequency sampling rate δ f states the formal expression of the spectral resolution of the DFT (or frequency resolution), i.e. the shortest separation in frequency between two subsequent samples in the spectrum:

δ f := f k+1 -f k = f s N = 1 NT s = 1 T w (23)

III-3-c Important consequence on the DtFT

Because δ f = 1 T w , the FT pair δ δ f ( f ) ⇌ 1 δ f δ 1/δ f (t) may be written as:

δ δ f ( f ) ⇌ T w δ T w (t) (24)
Thus, we can establish the FT pair between the discrete frequency DtFT and the windowed discrete time signal as:

[Z sw ( f k ) = Z sw ( f )δ δ f ( f )] ⇌ z sw (t s ) * T w δ T w (t) (25)
The windowed discrete time signal is replicated every T w duration, which means that it is actually T w -periodic! Nevertheless, this operation affects the amplitude of z sw (t s ) as the signal is now scaling with T w . To avoid this, Z sw ( f k ) is to be normalized by T w .

Z sw ( f k ) → Z sw ( f k ) T w (26)
Hence, out of equ.22:

Z sw ( f k ) = 1 T w N-1 ∑ n=0 z sw [n]e -j2πn f k fs * W( f k ) * δ f s ( f k ) × δ δ f ( f )| k (27)
Introducing the explicit expression of W( f k ) for a rectangular windowing (equ.16), we have now:

Z sw ( f k ) = N-1 ∑ n=0 z sw [n]e -j2πn f k fs * sinc(T w f k ) * δ f s ( f k ) × δ δ f ( f )| k (28)
The consequence of the frequency sampling has properly been taken into accout, we get rid from the sampling operator δ δ f ( f )| k to lighten the previous equation, while keeping the frequency sampling in mind:

Z sw ( f k ) = N-1 ∑ n=0 z sw [n]e -j2πn f k fs * sinc(T w f k ) * δ f s ( f k ) (29)
The figure below completes the previous one and illustrates the above elements.

III-4-Truncating the f s -periodization

It was mentioned above that the N samples of z s (t s ) resulting from the windowing operation feature a single f s period of its DtFT. Thus, the f s -periodization is truncated to the frequency range f k ∈ -f s 2 ;

f s 2 .
The discrete frequency support of the DFT is then defined as:

f k := - f s 2 + k f s N , ∀k ∈ N, k ∈ [0; N -1] (30)
or equivalently, using the spectral resolution δ f :

f k := - f s 2 + kδ f , ∀k ∈ N, k ∈ [0; N -1] (31)
Because the frequency range is restricted to f k ∈ -f s 2 ;

f s 2 , the convolution with the Dirac comb operator δ f s ( f ) can be removed from the expression of the DtFT (cf. equ.29):

Z sw ( f k ) = N-1 ∑ n=0 z sw [n]e -j2πn f k fs * sinc(T w f k ), with f k := - f s 2 + k f s N , and k ∈ N, k ∈ [0; N -1] (32)
An equivalent expression is obtained when introducing the dimensionless frequency:

k := f k f s (33) Hence, k ∈ -1 2 ; + 1 2 -1 N . The expression of Z sw ( f k ) now becomes: Z sw ( f k ) = N-1 ∑ n=0 z sw [n]e -j2πn k * sinc(T w f k ), with k := - 1 2 + k N , and k ∈ N, k ∈ [0; N -1] (34)
The previous equation constitues a possible definition of the DFT, however the regular algorithm implemented in most of softwares to compute the DFT does not perform the straight implementation of equ.34. To earn computation time, a further step is considered.

III-5-Pre-implementation step

In equ.34, the summation is performed over negative and positive values of the dimensionless frequency k ∈ -1 2 ; + 1 2 -1 N . To avoid the handling of signed and unsigned numbers, the implementation actually performs a shift of k by 1/2, such that k ∈ 0; 1 -1 N , and therefore:

Z sw ( f k ) = N-1 ∑ n=0 z sw [n]e -j2πn k * sinc(T w f k ), with k := k N , and k ∈ N, k ∈ [0; N -1] (35)
The previous equation now constitues the DFT, as implemented in the regular algorithm "fft" implemented in most of softwares.

Important notes:

• In equ.35, the DFT does not primarily depend on f s . The numerical implementation does not require the value of that parameter to compute the DFT of a sampled signal. However, the quantitative interpretation of the DFT does mandatorily require to know f s ! The f sdependence of the DFT will actually be concealed in the abscissa axis of the spectrum, as well as in its normalization procedure (cf. below). • In equ.35, due to the shift of k that now varies in the range 0; 1 -1 N , the first N 2 samples of the DFT and of the frequency support stand for k ∈

0; 1 2 -1 N , hence f k = k f s ∈ 0; f s 2 -δ f
, that is the positive frequencies range. Because this range is restricted to f s 2 ; + f s 2 , the next N 2 samples of the DFT and of the frequency support k ∈ 1 2 ; 1 -1 N , will stand for the negative frequencies range: -f s 2 ; -δ f . The concept applies with the samples of z sw (t s ) and those of the time support too.

• It is most important to notice that, strictly speaking, the physical unit of Z sw ( f k ), as computed by equ.35, is similar to the physical unit of z s (t s ). Yet, the regular physical unit of the Fourier transform is the one of signal multiplied by a time. This issue has to do with the normalization of the DFT, which will be discussed hereafter. • The convolution with $ sinc(T_w f_k)$ is not effectively realized in a typical computation of the DFT, which generally only uses the sum ∑ _{n = 0} {N-1}z_{sw}[n]e {-j2πn k}. It merely reflects the fact that the (time) windowing of the original signal has two consequences on the computed spectra : it increases the width of the spectral features in the DtFT and the DFT, and it introduces some oscillations in the spectra, corresponding to the sidelobes of the window function. Apart from the natural rectangular window, other window shapes can be used, that can reduce either the widening of the spectral features, either the residual oscillations.

The figure below illustrates the last sequence to calculate the DFT.

IV-DtFT & DFT: synopsis

The figure below summarizes the discretization process of a continuous time signal z(t) yielding its DtFT (step 5) and lastly, its DFT (step 9).

V-Computation of the DFT

In most of softwares, the computation of the DFT, as given by equ.35, is performed by the instruction "fft". E.g. if z sw [n] is a windowed sampled signal (a "vector") consisting of N samples in Octave, its DFT will be given by the instruction fft(z sw ).

We hereafter give couple of examples of such an implementation, compare it to the result returned by the "fft" instruction, and highlight some traps to avoid.

V-1-Raw computation

We examplify the calculation of the DFT with a simple cosine waveform signal of period T 1 (frequency f 1 = 1/T 1 ) over a windowing duration: T w = 1 s:

z(t) = cos(2π f 1 t), ∀t ∈ [0; 1] ( 36
)
whose FT is well-known:

Z( f ) = 1 2 [δ( f -f 1 ) + δ( f + f 1 )] (37)
On the other hand, the signal being T 1 -periodic, it can be expanded in Fourier series, whose set of coefficients are real:

Z 1 = Z -1 = 1 2 (38)
Below is given an Octave code that implements equ.35 and compares the result to the one returned by the "fft" instruction. We discuss the DFT spectra in the monolateral range 0;

f s 2 (i.e. 0;

f s 2 -δ f ) for now.
Illustrative code: The code below computes half the bilateral spectrum of a sinusoidal signal, windowed over an integer number of periods, without normalization.

[1]: [1]:

######################################################## # Illustration of the fft instruction from hard coding # ######################################################## #####
The modulus and phase of the hard-coded DFT from equ.35 and those returned by the "fft" instruction are identical, indeed.

The modulii feature a unique peak at the frequency of the cosine waveform, f 1 (here 12 Hz) in the range 0; f s 2 (here [0; 256[ Hz). However, unexpectedly, the magnitude of the peak does not match the magnitude of the Fourier series coefficient |Z 1 | = Z 1 = 1 2 , but rather N 2 (here 256, since we have set N = 512), which is discussed in the section below.

V-2-Normalization: Fourier series & Fourier transform

The discrepancy between the magnitude of the modulus of the computed peaks and the theoretical, expected, one is due to the lack of normalization of the DFT, as computed by equ.35, or by the "fft" instruction.

The normalization of the DFT is usually performed upon division of the result of the "fft" instruction by the number of samples z s (t s ) consists of (e.g. N).

Doing so, the normalized DFT stands for the Fourier series coefficients of the T w -periodic, or finite, duration sampled signal z s (t s ), whose physical units are those of z s (t s ).

However, to use the FT to evaluate the (continuous) Fourier Transform of an infinite-duration signal, whose physical units are those of the signal times time (e.g. V.s), it is necessary to use the normalizing factor T s = N T w .

To normalize properly the DFT, we will remind:

• If "fft" is meant to compute the Fourier series coefficients of a T-periodic discrete time signal z s (t s ), then the normalization is performed upon division of the result of the "fft" instruction by the number of samples z s (t s ) consists of. The physical units of Fourier series coefficients are consistently those of z s (t s ). • If "fft" is meant to compute the DFT of a discrete time signal z s (t s ) of finite duration T, then the normalization is again performed upon multiplication of the result of the "fft" instruction by the number of samples z s (t s ) consists of. The physical units of Fourier series coefficients are consistently those of z s (t s ). • If "fft" is meant to compute the Fourier Transform of an infinite duration time signal, then the normalization is performed upon division of the result of the "fft" instruction by the sampling period T s = N T w used to sample z s (t s ). The physical units of the DFT are now consistent with those of the FT of the continuous time signal z(t).

• Note that if T w = 1 s, the two normalization processes are equivalent, since 1 T s = N T w = N. The corresponding modified Octave code showing both the properly normalized Fourier series coefficients AND DFT of the cosine waveform signal, calculated out the "fft" instruction only, is given below.

Illustrative code:

The code below computes half the bilateral spectrum of a sinusoidal signal, windowed over an integer number of periods, and compares the two ways of normalizing the DFT.

| = Z 1 = 1 2 .
As to the DFT, its spectrum is now consistent with that of the Fourier series coefficients, too. However, this coincidence between the two normalizations is due to the fact that T w = T w 0 = 1. If you change that value to kT w 0 , the amplitude of the fundamental frequency will increase by k. In fact, since we are clearly dealing with a periodic signal, there is no reason to use the normalisation by T s in this case, and the normalization by N should therefore be used.

V-3-Bilateral spectra

Whether we deal with the Fourier series coefficients of a T-periodic or time-limited signal, or the DFT of an infinite duration signal, their spectrum is naturally bilateral and spreads, as stated above, over -f s 2 ;

f s 2 , or equivalently over -f s 2 ;

f s 2δ f . Several approaches lead to a consistent bilateral representation. Ours is a three-steps procedure:

1. The bilateral frequency range is first defined. 2. It has been stated in section III-4 that the first N/2 samples of the DFT stand for the positive range of frequencies whereas the second half of the samples stand for the negative range of frequencies. In Octave, the instruction "fftshift" allows for swaping the second half of the elements of a vector to the first half part of it and vice versa. Therefore, applying the instruction "fftshift" to the DFT-calculated vector makes its representation consistent over the bilateral frequency range. 3. The bilateral spectrum can be plotted consistlently (modulus, phase, real part, imaginary part).

The figure and the Octave code below illustrate the bilateral representation of the Fourier series coefficients and DFT spectrum of the cosine waveform signal. [6]:

The spectra have been calculated consistently and are bilateral in the range -f s 2 ; + f s

2

(here [-256; +256[ Hz).

V-4-Phase issues

In the spectra calculated so far, the phase trace shows apparent random fluctuations. These stem from the numerical fluctuations of the imaginary and real parts of the DFT. Indeed, the phase is calculated out of the regular definition:

arg(Z sw ) = φ(Z sw ) := atan Im{Z sw } Re{Z sw } (39)
Therefore tiny fluctuations of the real part of the DFT, will introduce large fluctuations of the phase, but fluctuations of the imaginary part of the DFT can cause unwanted fluctuations of the phase too. To overcome that, we state a "zero-phase" criterion: if the modulus of a given sample of the DFT is larger than an arbitrary small value, then the sample is expected to be valuable and its phase must be computed consistently (equ.39), otherwise we force it to be zero. This can easily be implemented in Octave by using a boolean test condition when computing the phase of the DFT vector.

The Octave code below illustrates our "zero-phase" criterion on the Fourier series coefficients and DFT spectrum of the cosine waveform signal. [4]:

The bilateral phase spectrum now features only zeros, as expected for the cosine waveform signal.

We exemplify again the calculation of the phase, but now with a sine waveform signal of frequency f 1 over a similar windowing duration: T w = 1 s:

z(t) = sin(2π f 1 t), ∀t ∈ [0; 1[ ( 40
)
whose FT is well-known, too:

Z( f ) = 1 2 δ( f -f 1 )e -j π 2 + δ( f + f 1 )e +j π 2 (41)
and whose set of Fourier series coefficients are imaginary: [7]:

Z 1 = Z * -1 = 1 2 e -j π
The bilateral phase spectrum features the expected ±π/2 values (±1.57 rad) for the phase at the ∓ f 1 frequencies, as expected from the theoretical calculation of the FT of that signal, which states the consistency of our phase criterion.

V-5 "Seeing" the time & frequency shift

To figure out the time shift as well as the frequency shift effect that is implictly introduced by the DFT computation, we consider a basic signal consisting in a time-shifted Dirac peak by a factor τ, hence:

z τ (t) = z(t -τ) = δ(t -τ) (43)
Considering the property of the FT regarding time-shifted signals, we have the FT pair:

[z τ (t) = z(t -τ)] ⇌ [Z τ ( f ) = Z( f )e -j2π f τ ] (44)
Therefore, here, with our Dirac peak:

[z τ (t) = δ(t -τ)] ⇌ [Z τ ( f ) = e -j2π f τ ] (45)
The FT has a constant modulus of 1 and a phase scaling linearly with the frequency f , however with a slope scaling with -τ. The phase slope is therefore expected to be negative with τ > 0 and vice versa.

The figure and the Octave code below illustrate the time and frequency shift that the DFT computation intoduces on the above-defined time-shifted Dirac signal that we windowed and sampled over T w = 1 s.

##### Displays i=5; # Change that index from 1 to length(tau) to select the desired value of tau figure subplot(3,1,1), plot(f_n, abs(Z_s_DFT(i,:)),'o'), xlabel('f_n [Hz]'), →ylabel('|Z_s|'), title('Normalized DFT') subplot [START_REF] Dft\{signal\} | Theo. DFT') axis([fmin fmax -A*width A*width]) subplot(2,2,4), plot(f_n, imag(Z_s_DFT)), xlabel('f_n [Hz]'), →ylabel('Im\{Z_s\}') hold on, plot(f_n, imag(Z_s_DFT_ns),'r'), plot(f_n,imag(Z_s_DFT_theo),'g') axis([fmin fmax -1e-18 1e-18]) figure subplot(2,2,1), plot(f_n, abs(Z_s_DFT)), xlabel('f_n [Hz]'), ylabel('|Z_s|'), →title('Normalized DFT') hold on, plot(f_n, abs(Z_s_DFT_ns),'r'), plot(f_n,abs(Z_s_DFT_theo),'g') legend('DFT\{fftshifted signal\}','DFT\{signal\}[END_REF]1,2), plot(f_n, Phi_Z_s_DFT(i,:),'o--'), xlabel('f_n [Hz]'), →ylabel ('arg(Z_s) [rad]') legend(strcat('\tau=', num2str(tau(i)))) subplot [START_REF] Dft\{signal\} | Theo. DFT') axis([fmin fmax -A*width A*width]) subplot(2,2,4), plot(f_n, imag(Z_s_DFT)), xlabel('f_n [Hz]'), →ylabel('Im\{Z_s\}') hold on, plot(f_n, imag(Z_s_DFT_ns),'r'), plot(f_n,imag(Z_s_DFT_theo),'g') axis([fmin fmax -1e-18 1e-18]) figure subplot(2,2,1), plot(f_n, abs(Z_s_DFT)), xlabel('f_n [Hz]'), ylabel('|Z_s|'), →title('Normalized DFT') hold on, plot(f_n, abs(Z_s_DFT_ns),'r'), plot(f_n,abs(Z_s_DFT_theo),'g') legend('DFT\{fftshifted signal\}','DFT\{signal\}[END_REF]1,[START_REF] Dft\{signal\} | Theo. DFT') axis([fmin fmax -A*width A*width]) subplot(2,2,4), plot(f_n, imag(Z_s_DFT)), xlabel('f_n [Hz]'), →ylabel('Im\{Z_s\}') hold on, plot(f_n, imag(Z_s_DFT_ns),'r'), plot(f_n,imag(Z_s_DFT_theo),'g') axis([fmin fmax -1e-18 1e-18]) figure subplot(2,2,1), plot(f_n, abs(Z_s_DFT)), xlabel('f_n [Hz]'), ylabel('|Z_s|'), →title('Normalized DFT') hold on, plot(f_n, abs(Z_s_DFT_ns),'r'), plot(f_n,abs(Z_s_DFT_theo),'g') legend('DFT\{fftshifted signal\}','DFT\{signal\}[END_REF], plot(f_n, Phi_Z_s_DFT(i,:),'o--'), xlabel('f_n [Hz]'),

VI-Applications

Although the theoretical framework of the DFT has been established for the general case of noncausal signals, that is signals that exist for positive as well as negative time t ∈ -T w 2 ; + T w 2 , so far we have treated examples of signals built as causal signals, that is signals defined for positive time only.

In the examples given hereafter, we focus on periodic, or non-periodic, non-causal signals.

VI-1 Fourier series coefficients of a periodic, non-causal signal

We start with the case of periodic, non-causal signals.

VI-1-a Problem positioning

Signal 1:

To that end, we consider the former T 1 -periodic sine waveform signal defined over a windowing duration T w = 1 s, but we arbitrary define it as:

z(t) = sin(2π f 1 t), ∀t ∈ - 1 2 ; 1 2 (46)
Hence, the non-causal nature of that signal. The code below reproduces the previous code for that signal, except that, our signal being periodic, we only focus of the Fourier series coefficients spectrum. [10]:

Illustrative

[10]:

We have added a zoom in the [-15; +15] Hz range, where three samples only are visible due to the frequency resolution of the problem δ f = 1/T w = 1/T 1 = f 1 = 12 Hz. The magnitudes of the modulii of the peaks at ± f 1 are 0.5 as expected, but the phase peaks are surprinsigly reversed compared to the former calculation of the same signal. Owing to its periodicity, the spectrum should be identical, though.

VI-1-b Questions & answers So, what is the problem?

Answer:

The DFT does not make any assumption on the causal, or non-causal, character of the sampled signal z s (t s ). It is the user who decides whether the signal is to be plotted against a causal, or a non-causal, time interval. The DFT does not consider the actual time support {nT s } but only the (time) index n, which always begins at n = 0 (Python) or n = 1 (Octave). Therefore, the structural nature of the DFT is to assume that the signal z s (t sw ) is causal, i.e. spreading over a fully positive time interval (cf. section III-5).

In the example above, the DFT algorithm therefore interprets the signal as starting from a "false reference time" t f rt = 0 (the signal "goes down") and lasting a duration T w = T 1 . In this case, doing so, the operation falsifies the phase of the signal at the "true reference time" t trt = 0, as defined by the user for the temporal representation of the signal (the signal "goes up"). If the signal is interpreted as starting from t f rt = 0, it is immediately seen that the regular expression for it is z(t) =sin(2π f 1 t), and NOT z(t) = + sin(2π f 1 t), as we initially believed. Hence the inversion of the phase peaks in the spectrum by a factor of π.

But, why did this effect not occur when considering the intial time interval t ∈ [-1/2; +1/2]?

Answer:

Having a look to that signal, it is readily seen that, over that windowing, the phases of the signal at t f rt = 0 and t trt = 0 are similar. Hence, the consistent calculation of the phase of the DFT.

Is there a way to treat non-causal signals without having to care about the way the windowing is performed?

Answer:

Yes, if we know the signal is truly non-causal and that the way its DFT phase is calculated matters, we can force the non-causal character of the signal to be taken into account, as demonstrated below.

VI-1-c Forcing the non-causality to be considered in the DFT spectrum

One forces the first half of the signal z sw (t s ) (first N/2 samples) standing for the negative part of the time support to be placed at the end of signal vector. Doing so, the "false reference time" of the DFT is forced to match the "true reference time" of the signal, and the DFT can be computed consistently. This is achieved by using again the instruction "fftshift", but that is now applied to the signal z sw (t s ) itself, before the computation of the DFT. The figure and the code below exemplify that.

  figure plot(t_s, z_s, 'o-'), xlabel('time [s]'), ylabel('Signal z_s(t)'), →title('Temporal representation of the signal') [1]:

  figure plot(t_s, z_s, 'o-'), xlabel('time [s]'), ylabel('Signal z_s(t)'), →title('Temporal representation of the signal') [3]:

  figure plot(t_s, z_s, 'o-'), xlabel('time [s]'), ylabel('Signal z_s(t)'), →title('Temporal representation of the signal') [6]:

  figure plot(t_s, z_s, 'o-'), xlabel('time [s]'), ylabel('Signal z_s(t)'), →title('Temporal representation of the signal') [4]:

  figure plot(t_s, z_s, 'o-'), xlabel('time [s]'), ylabel('Signal z_s(t)'), →title('Temporal representation of the signal') [7]:

  figure for i=1:length(tau) stem(t_s, z_s(i,:)), hold on endfor xlabel('time [s]'), ylabel('Signal z_s(t)'), title('Temporal representation of →the signal') ans = -63.248 ans = -48.925 [3]:

  figure plot(t_s, z_s, 'o-'), xlabel('time [s]'), ylabel('Signal z_s(t)'), →title('Temporal representation of the signal') [13]:

  figure subplot(2,1,1), stem(f_n, abs(Z_s_FSC)), xlabel('f_n [Hz]'), ylabel('|Z_s|'), →title('Normalized Fourier series coefs.') subplot(2,1,2), stem(f_n, Phi_Z_s_FSC), xlabel('f_n [Hz]'), ylabel('arg(Z_s) →[rad]')

For values of τ smaller than half of the windowing duration T w 2 = 0.5, the slope of the DFT phase is negative (leaving aside the wrapping of the phase within the region [-π; +π]). And the slope gets more and more negative as τ increases. As τ crosses T w 2 , the slope of the DFT phase gets positive, whereas there should be no reason for that.

The only interpretation is to assume that τ → -τ. The DFT algorithm interprets the second half of the signal (N/2 samples) as standing for negative times. Consequently the second half of the frequency support stands for negative frequencies. Hence the necessity to built the frequency support consistently! Illustrative code: [START_REF] Dft\{signal\} | Theo. DFT') axis([fmin fmax -A*width A*width]) subplot(2,2,4), plot(f_n, imag(Z_s_DFT)), xlabel('f_n [Hz]'), →ylabel('Im\{Z_s\}') hold on, plot(f_n, imag(Z_s_DFT_ns),'r'), plot(f_n,imag(Z_s_DFT_theo),'g') axis([fmin fmax -1e-18 1e-18]) figure subplot(2,2,1), plot(f_n, abs(Z_s_DFT)), xlabel('f_n [Hz]'), ylabel('|Z_s|'), →title('Normalized DFT') hold on, plot(f_n, abs(Z_s_DFT_ns),'r'), plot(f_n,abs(Z_s_DFT_theo),'g') legend('DFT\{fftshifted signal\}','DFT\{signal\}[END_REF] endfor Signal 2: Let us now look at the situation where, instead of a 1 s-windowing, the windowing is now defined as T w = T 1 = 1/ f 1 , symetrically around 0. In other words:

The code below exemplifies that.

Illustrative code:

[10]: # Parameters N = 512; # number of samples f1= 12; # frequency of the cosine waveform signal in Hz Tw = 1/f1; # windowing duration now defined over one period only #In-built variables Ts = Tw/N; # sampling period fs = 1/Ts; # sampling frequency The phase is now represented consistently.

Obviously, that use of "fftshift" on a well-windowed non-causal signal, as initially defined (equ.43), does not modify its DFT, as shown in the code below.

VI-2 DFT of a non-periodic signal

We now consider the more generic case of non-periodic signals, hence the use of the FT formalism.

VI-2-a Causal signal Signal 3:

For that application, we consider a causal, time-limited but non-periodic (because of the random noise) signal that is built as follows. We consider the realistic situation of a relevant electrical signal, perturbed by an ambient electrical noise. The relevant signal is supposed to be a sinusoidal signal of frequency 17 Hz. The noise results from two components: a sinusoidal oscillation, whose frequency and amplitude are arbitrarily set to 50 Hz and 0.1 V, respectively, and a uniform random background (white noise), whose magnitude equals that of the 50 Hz oscillation. The amplitude of the relevant signal is 5 times smaller than that of the 50 Hz oscillation. Hence, the relevant signal is "drowned in noise".

The signal is windowed over T w = 1 s. The code below shows how powerfull the DFT is to detect the spectral components of the noisy signal.

VI-2-b Non-causal signal

Signal 4:

In this application, we consider the important case of the rectangle function, a typical example of non-causal signal.

VII-DFT implementation using Octave: summary

The table below gives the summary of the Octave syntax to compute the mono-, or bilateral, the DFT spectrum of a windowed discrete time signal z sw (t s ) consisting of N samples, whether the signal is periodic of not.

VIII-References

1 -The Fast Fourier Transform and Its Applications, Book by E. Oran Brigham (original ed. 1988) [0]: