Mohamad Mortada

Abdallah Makhoul
email: abdallah.makhoul@univ-fcomte.fr

Chady Abou Jaoude
email: chady.aboujaoude@ua.edu.lb

Hassan Harb
email: hassanharb@auce.edu.lb

David Laiymani
email: david.laiymani@univ-fcomte.fr

A Distributed Processing Technique for Sensor Data Applied to Underwater Sensor Networks

Keywords: Underwater Sensor Networks, Periodic Applications, Euclidean Distance, Data Aggregation, Real Sensor Data

Wireless sensor networks (WSN) present a low cost solution to enhance our lives. They allow a large variety of applications. One of the major challenges faced by WSN is that of energy saving. A well known efficient way to reduce energy consumption is data reduction. It consists in reducing the amount of data sensed and transmitted to the sink. Consequently, sensor data communication should be minimized in order to increase network lifetime. In this paper, we propose an energy-efficient data reduction technique based on a clustering architecture. Our objective is to identify and study data similarities at both sensor and cluster-head (CH) levels. At the first level, each sensor sends a set of representative points to the CH at each period, instead of sending the raw data. When data points are received by the CH, it uses the Euclidean distance in order to eliminate redundant data generated by neighboring sensor nodes, before sending them to the sink. To validate our approach, we applied our techniques on real underwater sensor data and we compared them with other existing data reduction methods. The results show the effectiveness of our technique in terms of improving the energy consumption and the network lifetime, without loss in data fidelity.

I. INTRODUCTION

Wireless sensor networks (WSN) provide a low cost solution to enhance our lives. They allow a large variety of surveillance applications (medical, environmental, smart city, etc.). Their main advantages are fast, easy deployment and low maintenance cost. One of the major challenges faced by wireless sensor networks is that of energy saving. Indeed, data transmission consumes most of the available sensor's energy [START_REF] Raghunathan | Emerging techniques for long lived wireless sensor networks[END_REF]. Furthermore, the periodic data collection in surveillance applications, produces a huge amount of data, which are usually redundant [START_REF] Harb | An enhanced k-means and anova-based clustering approach for similarity aggregation in underwater wireless sensor networks[END_REF]. Then, the transmission of such amount of data is very expensive in terms of energy. In this way, reduction of sensed data becomes an efficient way to reduce energy consumption in WSNs.

In this paper, our aim is to propose a new distributed and low complex sensor data processing technique. We adopt a cluster based network's topology, where, sensor data are processed and aggregated at intermediate nodes, i.e. called Cluster-Heads (CHs), before sending them to the sink. This architecture has been proved to be efficient in terms of scalability and energy saving. Therefore, we studied a two data reduction levels technique. It aims at optimizing the volume of transmitted data at each cluster by achieving aggregation at both sensor nodes and CH levels in a periodically way. At the first level, each sensor node transforms its set of collected data to a reduced set of representative points. Then, it sends the set of points to its CH at the end of each period. After receiving the sets of points from all its sensors, the CH searches the similarity between each pair of data points coming from two sensors, with a technique based on the Euclidean distance concept. This phase allows to eliminate the redundancy between the received set of points. Finally, each CH sends the set of final data points to the sink. It is important to notice, that from the representative set of point the sink will be able to reconstruct the whole set of data with minimal errors.

To evaluate our approach, we choose to apply our techniques to underwater sensor networks. We believe that monitoring the aquatic environments is becoming a requirement for offering a better understanding of marine life. In such networks, the main objective is to monitor and observe the different kinds of aquatic environments, then periodically send the collected data to the end user for analysing and studying purposes. Furthermore, although the acoustic technology used in underwater ensures a long distance data communication [START_REF] Javaid | Chain-based communication in cylindrical underwater wireless sensor networks[END_REF], it consumes most of the available sensor's energy which is usually a limited battery power. For all these reasons, we chose real underwater sensor data to be the test-bed for our proposal. In addition, we compared our techniques with existing data reduction methods. We show how the effectiveness of our approach in reducing data and saving energy while guaranteeing high level of information integrity.

The rest of the paper is organized as follows: Section II gives an overview on related works reported on data aggregation in UASNs. Section III presents our sensor data processing technique, where each sensor node computes its representative data vector. In section IV we provide a multisensor data aggregation technique at the cluster head level. In Section V we detail the simulations we have conducted on real underwater readings data with a discussion of the obtained results. Finally, we conclude our paper and provide our directions for future work in Section VI.

II. RELATED WORKS

In the literature, one can find various data reduction approaches based on in-network processing [START_REF] Harb | Energy-efficient sensor data collection approach for industrial process monitoring[END_REF]- [START_REF] Harb | A distance-based data aggregation technique for periodic sensor networks[END_REF], data compression or data prediction methods [START_REF] Ali | Approximate replication of data using adaptive filters in wireless sensor networks[END_REF]- [START_REF] Mirco | Processing continuous join queries in sensor networks: A filtering approach[END_REF]. They are based mainly on algorithms like least mean square [START_REF] Ankur | Adaptive stream resource management using kalman filters[END_REF], [START_REF] Mirco | Processing continuous join queries in sensor networks: A filtering approach[END_REF] and Kalman Filter [START_REF] Ali | Approximate replication of data using adaptive filters in wireless sensor networks[END_REF], [START_REF] Abdolee | Diffusion lms algorithms for sensor networks over non-ideal inter-sensor wireless channels[END_REF]. Although these approaches predict sensed values and allow efficient data reduction, however they present several disadvantages. They are computationally complex, sometimes they generate communication overhead, and the sink may need some transmissions to detect failures.

To reduce the amount of data transmitted in the network, other techniques have been proposed and consist in data aggregation. Data aggregation aims at eliminating redundancy in data collected and minimizing the number of transmissions, thus saving the overall network energy. Recently, the majority of the proposed data aggregation techniques have been built with the clustering scheme which has been proven to be an efficient way in terms of scalability and data traffic [START_REF] Harb | Comparison of different data aggregation techniques in distributed sensor networks[END_REF].

In [START_REF] Wang | Energy-efficient distributed compressed sensing data aggregation for cluster-based underwater acoustic sensor networks[END_REF], [START_REF] Lin | Energy-efficient compressed data aggregation in underwater acoustic sensor networks[END_REF], the authors study the data aggregation in UASNs as a compression scheme for data generated in each cluster. The authors in [START_REF] Wang | Energy-efficient distributed compressed sensing data aggregation for cluster-based underwater acoustic sensor networks[END_REF] and [START_REF] Wang | Energy-efficient data aggregation scheme for underwater acoustic sensor networks[END_REF] propose two data aggregation schemes, namely block diagonal matrix and block upper triangular matrix, for cluster-based UASNs inspired by the Distributed Compressed Sensing (DCS) technique . The main objective of such schemes is to generate RIP-preserving (Restricted Isometric Property) measurements of sensor readings by taking multi-hop underwater acoustic communication cost into account. Finally, a distributed compressed sensing reconstruction algorithm, called DCS-SOMP, is adopted to recover raw sensor readings at the fusion centre.

Some works in data aggregation in UASNs, such as [START_REF] Tran | Uwsns: A round-based clustering scheme for data redundancy resolve[END_REF], are based on the formation of clusters and the selection of cluster-heads. The authors in [START_REF] Tran | Uwsns: A round-based clustering scheme for data redundancy resolve[END_REF] propose a data aggregation round-based clustering scheme in order to reduce the transmission of redundant data in UASNs. The proposed scheme works in rounds where each round consisting in four main phases: initialization, cluster-head selection, clustering, and data aggregation. By applying some mechanisms in each round, the proposed scheme reduced the energy consumption in the network and minimized data redundancy, while still guaranteeing data accuracy. In [START_REF] Kui | Energy balanced clustering data collection based on dominating set in wireless sensor networks[END_REF], the authors propose EBDSC, a distributed Energy-Balanced Dominating Set-based Clustering scheme, to extend the network lifetime by balancing energy consumption among different nodes. In EBDSC, a node becomes a cluster head candidate if it has the longest lifetime among its neighbours.

Other data aggregation and collection studies, such as [START_REF] Lloret | An underwater wireless group-based sensor network for marine fish farms sustainability monitoring[END_REF], [START_REF] Goyal | Fuzzy based clustering and aggregation technique for under water wireless sensor networks[END_REF], have been proposed. In these works the computation of statistical means and moments summarize the data obtained by the UASNs. In [START_REF] Lloret | An underwater wireless group-based sensor network for marine fish farms sustainability monitoring[END_REF], the authors propose an analytical model group-based sensor network in order to monitor the accurate amount of pollution that is deposited on the seabed. The objective is to study the effects produced by feed loss in the marine fish cages and its environment impact. After searching the best location to place the sensor nodes, the proposed model determines the amount of food that is wasted while it measures the amount of generated deposits. The authors in [START_REF] Goyal | Fuzzy based clustering and aggregation technique for under water wireless sensor networks[END_REF] propose to design a fuzzy based clustering and aggregation technique for UWSN. In this technique the parameters such as the residual energy, the distance to sink, the node density, the load and the ilnk quality are considered as input to the fuzzy logic. Based on the output of fuzzy logic module, appropriate cluster heads are elected and act as aggregator nodes.

Finally, works in [START_REF] Bahi | A two tiers data aggregation scheme for periodic sensor networks[END_REF], [START_REF] Tran | Well-suited similarity functions for data aggregation in cluster-based underwater wireless sensor networks[END_REF] are dedicated to periodic applications in sensor networks. The authors use some similarity functions to aggregate data generated in the networks. The main objective is to eliminate redundancy and reduce the size of data transmitted thus, optimizing the energy consumption and reducing overload on the network level. Further to a local processing at the sensor node level, the authors in [START_REF] Bahi | A two tiers data aggregation scheme for periodic sensor networks[END_REF] propose a prefix frequency filtering (PFF) technique at the CH level. PFF uses Jaccard similarity function to identify similarities between near sensor nodes at each period and integrates their sensed data into one record. Then, several versions of PFF, i.e. PSFF [START_REF] Harb | A suffix-based enhanced technique for data aggregation in periodic sensor networks[END_REF] and KPFF [START_REF] Harb | Kmeans based clustering approach for data aggregation in periodic sensor networks[END_REF], have been proposed in order to optimize the data latency. On the other hand, the authors in [START_REF] Tran | Well-suited similarity functions for data aggregation in cluster-based underwater wireless sensor networks[END_REF] use distance functions, such as Euclidean and Cosine, at the CH level to build an efficient underwater network by reducing packet size and by minimizing data redundancy. However, although the proposed techniques eliminate the similar data, some redundancy still remain in the final data sets sent to the sink.

In this paper, we propose a new less complex data reduction and aggregation technique suitable for low resources sensor networks. In this technique, data aggregation are performed at both sensor and CH levels where we transform the raw data to a set of representative data points after eliminating redundancy among them. Compared to the existing techniques, our technique is more efficient to reduce the redundancy among raw data and thus, to preserve the energy in the network.

III. SENSORY DATA PROCESSING

In this section we present our sensor data processing technique, to be executed by sensor nodes in order to find the representative data points. First, we present the network's topology that we consider in our approach.

A. WSN topology

In this paper, we consider a cluster-based architecture for the network. The cluster-based architecture is based on grouping sensor nodes into clusters, and assigning for each cluster a super node (intermediate node), the cluster head (CH). The CH is elected after the network deployment and can be changed dynamically during the network lifetime. It can be a regular node or a specific more powerful one. Data transmission between sensor nodes and their appropriate CH is based on single-hop communication as presented in Fig. 1. In this work, we use the periodic data collection approach, in which each sensor node sends periodically (period p) its data to the appropriate CH, which in his turn sends them to the sink. Then, Algorithm 1 Local Aggregation Recursive Algorithm.

Require: Vector of readings: R i , distance threshold: t d , start index: start = 1, end index: end = τ . Ensure: List of representative points of R i : P i .

1: P i ← ∅; // list of empty points

2: E d = end index=start (index-start)×(R[end]-R[start]) end-start + R[start] -R index 2 3: if E d ≤ t d then 4: P i ← P i ∪ {(start, R[start])} ∪ {(end, R[end])} 5: else 6: P i ← P i ∪Local Aggregation(R, t d √ 2 , start, start+end
P i ← P i ∪Local Aggregation(R, t d √ 2 , start+end 2 , end) 8: end if 9: return P i D i = τ i=1 r 2 i , where r i ∈ R i (3)
Finally, each sensor node S i sends its list of representative points P i and its radius D i to the CH, at the end of each period. In the next section, we describe how the CH will aggregate the data coming from its member nodes before to send them to the sink.

IV. MULTI-SENSOR DATA SIMILARITY SEARCHING

At the end of each period, the CH receives the sets of points with their corresponding radiuses coming from its member nodes. The objective is then to identify all pairs of member nodes that generate redundant sets in order to eliminate duplication before sending them to the sink. In the previous section, we considered that two sets are similar if the Euclidean distance between them is less than the threshold t d . However, applying the Euclidean distance for every pair of sets is very expensive in terms of computation since it generates O(n 2 ) number of comparisons, where n is the number of received sets. In addition, the computation will be more complex for large data sets as in the case of dense sensor networks. Therefore, in order to reduce the number of comparisons, it is mandatory to search the pairs of redundant sets. This search will be performed in two phases. In the first phase we compute a list of pairs which are "candidates" to be similar. A pair is candidate if it satisfies some conditions and it means that the two sets composing this pair may be similar. However, a pair is not candidate means that it is for sure not similar. To ensure the similarity of candidate pairs, we need a verification phase. This verification is necessary since different sets of points coming from different sensor nodes may be of different size (see sub-section IV-B).

A. The candidate pairs generation phase

In this phase, each CH computes the pairs of sets (set of points or vectors) which are "candidate to be similar". Our intuition is that if the distance between the radiuses of two sets of points is less than the threshold t d then, the Euclidean distance between the two sets of points is candidate to be less than t d . Therefore, in our work, we prove that two sets of points P i and P j are candidates if and only if the distance between their corresponding radiuses is less than t d as shown in the following lemma:

Lemma 2: Consider two sets of points P i and P j with their corresponding radiuses D i and D j respectively. Assume that R i and R j are the initial data vectors of P i and P j respectively. Thus, if the Euclidean distance between R i and R j is less than the distance threshold t d then, the distance between their corresponding radiuses should be also less than t d . Therefore:

P i and P j are considered candidates ⇐⇒ |D i -D j | ≤ t d
Proof 2: This lemma can be simply demonstrated based on the proof of lemma 1.

Algorithm 2 describes how each CH searches the set of candidates for each sensor. It takes as input a collection of data points sets with their radiuses coming from different sensor nodes. It scans sequentially each set of points P i (line 2) and selects the candidates based on the lemma 2 (line 5). Afterwards, P i and all its candidates will be verified according to the Euclidean distance threshold (line 7 and sub-section 4.2). Finally the algorithm returns a list of n-uplet where for each points sets P i are associated the sensors S j "candidate" to be similar to S i . F i ← S i ; // list of similar sensors to S i

4:

for each set of points P j ∈ P such that j > i do 5:

if |D i -D j | ≤ t d then 6:
// P i and P j are candidates 7:

if Euclidean Distance(P i , P j ) ≤ t d then 8:

F i ← F i ∪ {S j } 9:
P ← P -{P j } // remove P j from P S ← S ∪ {(P i , F i )} 14: end for 15: return S

B. Candidate pairs' verification

As previously exposed, two sets of points P i and P j in a candidate pair are considered similar if their distance is less than the distance threshold t d (line 7 in Algorithm 2). However, P i and P j can have different sizes, i.e. number of points. This property makes the computation of the Euclidean distance not trivial.

Algorithm 3 Euclidean Distance Computation Algorithm.

Require: P i = {p i1 , p i2 , . . . , p in i }, P j = {p j1 , p j2 , . . . , p jn j }, p k = (x k , y k ). Ensure: Euclidean distance between P i and P j : E d (P i , P j ).

1: distance = 0 2: for each points {p i , p i+1 } ∈ P i do 3:

if {p j , p j+1 } ∈ P j exists such that x i = x j and x i+1 = x j+1 then 4:

// find equations of the two lines based on equation 45:

(L i,i+1 ) :

y i = a i × x + b i 6:
(L j,j+1 ) :

y j = a j × x + b j 7:
// calculate the Euclidean distance between lines based on equation 58: // create a new point then add it to P j 12:

e d = [E d (L i,i+1 , L j,j+1 )]
p jnew = {x i+1 , y j (x i+1 )} 13:
insert p jnew just before p j+1 in P j 14:

// find equations of the two lines based on equation 415:

(L i,i+1 ) :

y i = a i × x + b i 16:
(L j,jnew ) : with 120 sensors. Thus, data collected by the sensors are sent to their appropriate CHs which are located geographically at the centre of the clusters. We compare our results to those obtained with the technique proposed in [START_REF] Tran | Well-suited similarity functions for data aggregation in cluster-based underwater wireless sensor networks[END_REF], which we will refer as EuDi, used for UASNs. We choose to compare our works to EuDi because the two architectures are the same and because the results obtained by EuDi are well positioned in the state of the art.

y j = a j × x + b j 17: e d = [E d (L i,i+1 , L j,jnew )]
In our simulations, we evaluated the performance using the following parameters:

• the period size, τ , takes the following values: 128, 256, 512 and 1024.

• the similarity between two readings, δ, takes the following values: 0.001, 0.005, 0.01 and 0.025.

A. Percentage of data sent periodically from sensors to CH

Fig. 3 shows the percentage of data sent from each sensor to its CH at each period, after transforming its raw data into a set of points. Compared to the EuDi technique, these results show that, by sending its set of points, a sensor can reduce the amount of transmitted data by 60% and up to 97% . Therefore, our technique can successfully minimize the data transmission to the CH by eliminating redundancy among sensor's raw data.

Furthermore, we observe that sensor nodes send less data when δ or τ increases. This can be explained by the fact that when the similarity between collected readings increases, the number of representative points is also be reduced.

B. Number of iterations required to compute a set of points

Fig. 4 presents the average number of iterations required at each period to find the final set of representative points obtained in algorithm 2. It is important to recall that a high number of iterations can increase the complexity of the proposed algorithm as well as the data latency at the sensor. The obtained results show that, the number of iterations in our technique is almost less than 20, except when δ = 0.001 where it exceeds 40 iterations. We think that these results are suitable for the most kinds of sensors where the parameters values should be determined by the decision makers depending on the application requirements.

C. Number of candidates / comparisons

Fig. 5 shows the number of compared sets without applying our technique (i.e. with naïve comparisons between every pair of sets), the number of candidates generated by our technique and the results obtained after applying the Euclidean distance function (the real number of similar sets).We fixed the period size to 1024 and we varied δ as shown in the figure. We notice that, the number of comparisons in our technique is largely minimized compared to the naïve comparison. This is due to the lemma 2 which prune out the infeasible nonsimilar data sets and limits the number of comparisons to the candidates sets. Moreover, it is important to notice that the number of comparisons in our technique is more minimized when δ increases, thus, the number of candidates tend to match the exact number of similar sets. D. Percentage of final data sent to the sink Fig. 6 shows the percentage of final data sent from each cluster to the sink after aggregating data at sensor and CHs levels. We fixed the value of δ to 0.005 while we varied τ from 128 to 1024. The obtained results show that the data collected in each cluster have been largely reduced using our technique and compared to EuDi, for all values of τ . For instance, using our technique, the percentage of data sent CH 1 CH 2 does not exceed, in the worst case, 4% of the whole collected data. Otherwise, EuDi technique can reduce, in the best case, up to 50% the data collected in each cluster. Therefore, our technique can efficiently minimize the overload in the network and send only the useful information to the sink, without loss of the data integrity.

E. Aggregation process time at the CH

In this section, we compare the execution time required for the aggregation process at the CHs level when fixing δ to 0.005 and varying the period size (Fig. 7). The test machine consisted in a standard Core i5 (1.6Ghz) laptop running Windows 7 operating system with 4Gb of RAM. Clearly this configuration appears to be more powerful than a typical CH. Nevertheless, the aim of these tests is to evaluate and compare our technique. In this way, compared to EuDi the obtained results show that our technique can accelerate the aggregation process from 3 (when the period size is big) to 20 times (when the period size is small). This is due to two reasons: first, the Euclidean distance is applied, in our technique, over the candidate pairs only while it is applied, in EuDi, over all pairs of sets; second, the Euclidean distance is calculated, in our technique, between the representative points while it is calculated, in EuDi, between all the readings of two sets. Therefore, we can deduce that our technique can minimize the data latency at the CHs and deliver data to the sink in a faster way.

F. Energy consumption

In this section, we study the energy consumption of our technique at both sensors and CHs levels, for different values of δ and τ . For this, we use the energy model described in [START_REF] Harb | An enhanced k-means and anova-based clustering approach for similarity aggregation in underwater wireless sensor networks[END_REF]. In Fig. 8, we show the energy consumption in each sensor when fixing δ and varying τ every time. The obtained results show that our technique greatly outperforms EuDi in terms of preserving the energy in the sensors in all tested cases. As we can see, the energy consumption in the sensor can be minimized, using our technique, up to 96% compared to the EuDi technique. Such energy optimization is obtained since the amount of data transmitted has been largely reduced in our technique (see Fig. 3). The energy consumption at the CH level, i.e. CH 2 , using our technique and EuDi is presented in Fig. 9. The results are dependent on, first, the amount of data received from the sensor members (Fig. 3) and, second, the amount of final data sent to the sink after eliminating the redundancy among them (Fig. 6). We can show clearly that our technique significantly reduces, e.g. up to 91%, the energy consumption in CH 2 compared to EuDi. Therefore, our technique can be considered as a very efficient technique since the energy consumption is highly reduced in the network while the information integrity is fully preserved at the sink node.

VI. CONCLUSIONS AND FUTURE WORK

In this paper, we have proposed an energy-efficient data reduction technique dedicated to surveillance applications. After reducing the size of the collected data, each sensor node sends a set of representative points to the CH at the end of each period. When these sets arrive to the CH, this last uses the Euclidean distance to eliminate redundant data generated by neighbouring sensor nodes, before sending it to the sink. Comparing to other existing data aggregation techniques, simulation results on real sensor data show the effectiveness of our technique in terms of energy consumption, while still keeping a high quality of the collected data. We have two major directions for our future work. First, we plan to improve the polygonisation technique to achieve the search of minimum line number corresponding to a given data curve. Thus, sensor nodes will conserve more energy and network lifetime will be extended. Second, we seek to adapt our technique to take into consideration reactive periodic sensor networks, where sensor nodes operate with different sampling rate.

2

 2

Algorithm 2

 2 Candidate Pairs Searching Algorithm. Require: Sets of points: P = {P 1 , P 2 , . . . , P n }, Set of radius: D = {D 1 , D 2 , . . . , D n }, distance threshold: t d . Ensure: List of data points sets with duplicated data sensors for each one. 1: S ← ∅ 2: for each set of points P i ∈ P do 3:

2

 2

2

 2

19 :

 19 distance = distance + e d 20: end for 21: return √ distance

Fig. 3 .

 3 Fig. 3. Percentage of data sent periodically from each sensor to CH 1 .

Fig. 4 .

 4 Fig. 4. Average number of iterations required to compute a set of points.

Fig. 5 .

 5 Fig. 5. Data sets comparison at CH 1 , τ = 1024.

2 Fig. 6 .

 26 Fig. 6. Percentage of final data sent to the sink, δ = 0.005.

2 Fig. 7 .

 27 Fig. 7. Aggregation process time at CHs, δ = 0.005.

Fig. 8 .

 8 Fig. 8. Energy consumption in each sensor.

Fig. 9 .

 9 Fig. 9. Energy consumption in CH 2 .

ACKNOWLEDGEMENT

his work has been supported by the EIPHI Graduate School (contract "ANR-17-EURE-0002") and partially funded with support from the Hubert Curien CEDRE programme, Project RESEAUCORP, Number: 40283YK.


Student member, IEEE Weike Feng

Jean-Michel Friedt

Giovanni Nico
email: g.nico@ba.iac.cnr.it

Fellow, IEEE Motoyuki Sato
email: motoyuki.sato.b3@tohoku.ac.jp

Three-dimensional

come



Abstract-We designed a ground based radar system with a C-band two-dimensional cross multiple input multiple output (MIMO) array for three-dimensional imaging and displacement estimation purpose. For this system, we developed a far-field pseudo-polar image format algorithm using pseudo-polar spherical coordinate. The use of a tensor compressive sensing technique allows to focus undersampled raw data and to optimize the data acquisition time and memory usage. A novel algorithm, named as tensor based iterative adaptive approach, is proposed for the effective and efficient reconstruction of sparse targets with a reduced level of sidelobes. Experiment results validate the designed radar system and the proposed algorithms.

Index Terms-Three-dimensional (3D) imaging, cross multiple input multiple output (MIMO), compressive sensing, pseudo-polar spherical coordinate, iterative adaptive approach.

I. INTRODUCTION

ROUND based synthetic aperture radar (GBSAR) is a powerful technique to monitor displacement, deformation, and vibration of terrain and infrastructure with an accuracy of sub-millimeter [START_REF] Di Pasquale | Monitoring strategies of earth dams by ground-based radar interferometry: How to extract useful information for seismic risk assessment[END_REF][START_REF] Nico | Measurement of Pier Deformation Patterns by Ground-Based SAR Interferometry: Application to a Bollard Pull Trial[END_REF]. In GBSAR, a one-dimensional (1D) synthetic aperture is normally formed by sequentially moving the transceiver on a rail to get a high angle resolution. In recent years, linear multiple input multiple output (MIMO) array based methods have been proposed [START_REF] Tarchi | MIMO radar and ground based SAR imaging systems: Equivalent approaches for remote sensing[END_REF][START_REF] Feng | Near range radar imaging based on block sparsity and cross-correlation fusion algorithm[END_REF], which can reduce the data acquisition period of GBSAR, benefiting some specific applications where it is necessary to do many measurements within a short time. Besides, a shorter data acquisition time is also useful to mitigate the phase artifacts caused by microwave propagation in atmosphere on the displacement estimation.

Characterized by a highly suboptimal aperture compared with the airborne/spaceborne SAR, several imaging algorithms were developed to focus GBSAR data, including time-domain back projection algorithm (TDBP), frequency-domain back projection algorithm, range migration algorithm, and far-field This work was supported by JSPS Grant-in-Aid for Scientific Research (A) 26249058.

Weike Feng is with the Graduate School of Environmental Studies, Tohoku University, 980-8579, Sendai, Japan (e-mail: feng.weike.q4@dc.tohoku.ac.jp).

Jean-Michel Friedt is with the FEMTO-ST, Time & Frequency department, Besancon, France (e-mail: jean-michel.friedt@femto-st.fr).

pseudo-polar format algorithm (FPFA) [START_REF] Guccione | Focusing algorithms analysis for ground-based SAR images[END_REF]. For MIMO array based radar, these algorithms can be adopted with simple modifications [START_REF] Zhuge | Frequency-wavenumber domain focusing under linear MIMO array configurations[END_REF][START_REF] Feng | Near range radar imaging by SFCW linear sparse array[END_REF]. Compared with other imaging algorithms, despite being only applicable in the far field, FPFA has the great advantage of computation complexity, especially when the aperture length is comparable to the range resolution [START_REF] Fortuny-Guasch | A fast and accurate far-field pseudopolar format radar imaging algorithm[END_REF][START_REF] Han | Generalized pseudopolar format algorithm for radar imaging with highly suboptimal aperture length[END_REF].

Since only 1D aperture is synthesized, GBSAR can only get a two-dimensional (2D) focused SAR image. In such a case, the 3D scene is projected onto a 2D imaging plane, resulting in some information losses. Normally, for the applications of GBSAR to dam and landslide monitoring, an accurate digital elevation model (DEM) is needed for the accurate positioning of targets that have large displacements [START_REF] Di Pasquale | Monitoring strategies of earth dams by ground-based radar interferometry: How to extract useful information for seismic risk assessment[END_REF]. However, the unambiguity requirement always makes the projection of the 2D GBSAR image onto the DEM a complicated task. Moreover, an accurate DEM is not always available. Therefore, a GBSAR system that can get the 3D focused image is desired, which requires a 2D synthetic aperture. Two simple ways are possible to get a 2D aperture: 1) mechanically move the transceiver on a 2D scanner, and 2) similarly with the MIMO-SAR system [START_REF] Zhuge | A sparse aperture MIMO-SAR-based UWB imaging system for concealed weapon detection[END_REF][START_REF] Zhu | Range migration algorithm for near-field MIMO-SAR imaging[END_REF], using a linear MIMO array moved on a mechanical rail. However, the data sampling time of these two approaches is normally too long for many applications.

In this paper, we present a 3D volumetric imaging GBSAR system with a 2D cross-MIMO array. Compared with other 2D MIMO configurations, cross-MIMO is easier to implement in practice. According to the signal model in far-range condition, we established a pseudo-polar spherical coordinate, based on which the FPFA algorithm is extended to its 3D version for effective 3D imaging. When the imaging scene is sparse, i.e. only several targets have significant reflection coefficients, we proposed a tensor compressive sensing (T-CS) [START_REF] Qiu | Three-dimensional sparse turntable microwave imaging based on compressive sensing[END_REF][START_REF] Caiafa | Multidimensional compressed sensing and their applications[END_REF] based imaging algorithm to shorten data acquisition time, reduce memory usage, and improve imaging quality. Furthermore, a tensor based iterative adaptive approach [START_REF] Yardibi | Source localization and sensing: a nonparametric iterative adaptive approach based on weighted least squares[END_REF][START_REF] Jahromi | Two-dimensional iterative adaptive approach for sparse matrix solution[END_REF][START_REF] Feng | Jointly iterative adaptive approach based space time adaptive processing using MIMO Radar[END_REF] (T-IAA) is proposed to solve the T-CS problem efficiently and accurately. Experiment results are presented to show the effectiveness of the designed radar system and the proposed algorithms.

II. RADAR SYSTEM AND SIGNAL MODEL

The configuration of the designed radar system is shown in Fig. 1, which mainly includes a vector network analyzer (VNA) to generate the frequency stepped continuous wave (FSCW) signal, four 1×8 switches to realize the time-divided signal transmitting and receiving, two micro-controllers used to control the switches, a uniform receiving array with 16 Vivaldi antennas, a uniform transmitting array with 16 spiral antennas, and a processing PC used to save and process the VNA data and to remotely control the micro-controllers and the VNA. The imaging geometry of the designed system is shown in Fig. 2, where the i-th receiver is at (xi, 0, 0), the j-th transmitter is at (0, 0, zj), i=1,2,…,16, and j=1,2,…,16. With a target located at (x0, y0, z0), the distance between the i-j-th sampling point and the target is given by

0 2 2 2 2 2 2 , 0 0 0 0 0 0 ( ) + ( ) i j i j R x x y z x y z z        (1)
Fig. 2. Imaging geometry of the designed radar system.

Then, with the FSCW transmitted signal, the received signal at i-j-th sampling point and q-th frequency can be expressed as

0 0 0 , 0 ( , , ) exp( j2 / ) ( , , ) q i j s i j q f R c n i j q     (2)
where fq=f0+(q-1)Δf, q=1,2,...,Q, f0 is the start frequency, ∆f is the frequency step, α0 is the reflection coefficient of the target, c is the velocity of light, and 0 ( , , ) n i j q is the additive noise. Therefore, the received signal from the whole illuminated scene is given by , j2

( , , ) / ( , , ) ( , , ) d d + ( , , )

q i j f R x y z c s i j q x y z e x ydz n i j q       (3)
where , ( , , ) ij R x y z is the distance between the target at (x, y, z)

and i-jth sampling point, ( , , )

x y z  is the reflection coefficient, and  denotes the illuminated scene.

According to [START_REF] Feng | Near range radar imaging based on block sparsity and cross-correlation fusion algorithm[END_REF], the reflection coefficient of target at (x, y, z) can be estimated by the BP algorithm as follows:

, +j2 ( , , ) / 1 1 1
( , , ) ( , , )

q i j Q IJ f R x y z c i j q x y z s i j q e        (4)
However, the estimation of the entire illuminated scene using (4) has a high computation complexity, which prevents it from practical applications. Therefore, a faster algorithm is needed.

III. FAST 3D IMAGING ALGORITHM

A. 3D Pseudopolar coordinate

In Fig. 2, let

i j i j R R x z   (5)
Therefore, the reflection coefficient of the target at (R, θ, φ) can be estimated by

,, sin sin 2 j2 j2 j2 j 1 1 1 ( ) ( , , ) ,, ij q i j q cc R Q IJ xz f c i j q R s i j q e e e e               (6)
where ,, 2

(sin sin ) / i j q i j q f x z c      , qc q f f f , c
f is the center frequency, and

/ cc cf  
is the wavelength. The absolute value of ,,

i j q  in (6) is bounded by ,, (sin sin ) 2
(sin sin ) / 4

xz i j q i j q LL f x z c R             (7)
where

/2 R c B 
is the range resolution, B is the bandwidth, Lx and Lz are the synthetic aperture length in x direction and z direction, respectively. As analyzed in [START_REF] Han | Generalized pseudopolar format algorithm for radar imaging with highly suboptimal aperture length[END_REF], if ,, i j q



satisfies the condition (8), the last exponential term of ( 6) can be ignored. max max ,, ( sin sin ) 1 22

xz i j q LL R          (8)
For the designed radar system, max max

( sin sin ) xz LL   is close to 2 R  .
Besides, since the mean value of xi (zj) is zero and due to the factors sinθ and sinφ, ,, i j q  is smaller than those given in (8) [START_REF] Fortuny-Guasch | A fast and accurate far-field pseudopolar format radar imaging algorithm[END_REF]. Therefore, the reflection coefficient of the target at (R, θ, φ) can be approximately estimated by

sin sin 2 j2 j2 j2 1 1 1 ( ) ( , , , ) , ij q cc R Q IJ xz f c i j q s i j q e e R e            (9)
It can be observed that the exponential terms in (9) form the kernel of a 3D Fourier transform. Therefore, a 3D pseudopolar spherical coordinate ( , , )

   is defined as

2 / , sin / , sin / c c Rc          
(10) Based on this pseudopolar coordinate, the received signal cube and the 3D FPFA algorithm can be formulated as the following multilinear expression:

1 1 2 2 3 3      S F F F (11) and   H H H 1 1 2 2 3 3 3D FFT       S F F F S ( 12
)
where (• ) H denotes conjugate transpose,   denotes the mode- tensor by matrix product [START_REF] Kolda | Tensor decomposition and applications[END_REF],

S is the received signal cube,  is the reflection coefficient tensor of the imaging scene, F1, Compared to ( 4) and ( 6), the proposed 3D FPFA focusing algorithm [START_REF] Qiu | Three-dimensional sparse turntable microwave imaging based on compressive sensing[END_REF] can significantly reduce the computational cost by using the fast Fourier transform (FFT). However, due to the intrinsic limitation of Fourier transform based methods, highlevel sidelobes will be generated by [START_REF] Qiu | Three-dimensional sparse turntable microwave imaging based on compressive sensing[END_REF], which will blur the weak targets and thus introduce negative influence on the following displacement estimation.

B. Tensor based IAA

When the imaging scene is sparse, i.e. only few strong targets exist, T-CS based method can be used for 3D imaging as shown in [START_REF] Caiafa | Multidimensional compressed sensing and their applications[END_REF], which can reduce the sidelobe level and achieve higher resolution than Fourier transform based imaging methods [START_REF] Qiu | Three-dimensional sparse turntable microwave imaging based on compressive sensing[END_REF].

0 1 1 2 2 3 3 =min || || s.t. || || F          S F F F ( 13
)
where  denotes the noise level, ||• ||0 denotes the number of nonzero elements in a vector or a tensor, and the Frobenius norm of a tensor Y is defined as

1 2 3 1 2 3 2 ,, || || | | F i i i i i i y    Y .
One advantage of T-CS based imaging method is that it can still work well with under-sampled data, which can help to reduce the data acquisition time, the computation complexity, and the memory usage of [START_REF] Caiafa | Multidimensional compressed sensing and their applications[END_REF]. For the designed radar system, to reduce the spatial sampling points in x and z directions, a data acquisition program is used to control the switches to randomly turn on several transmitter-receiver pairs of the cross MIMO array. Since the frequencies of the FSCW signal cannot be randomly generated by VNA, the frequency under-sampling is realized in the data processing step. The random selection of frequencies by using a discrete frequency synthesizer driven by a set of random indication numbers is under study.

The under-sampled signal model, i.e. the compressive measurement of target reflection coefficient tensor, is shown in Fig. 3, where Then, the focused image of targets can be obtained by Eq. ( 14) is NP-hard but can be solved by some efficient algorithms, such as the Kronecker-orthogonal matching pursuit (Kron-OMP) algorithm and the N-way Block OMP algorithm (NBOMP) in [START_REF] Caiafa | Multidimensional compressed sensing and their applications[END_REF], and the 3D smoothed l0 norm algorithm (3D-SL0) in [START_REF] Qiu | Three-dimensional sparse turntable microwave imaging based on compressive sensing[END_REF]. However, Kron-OMP and NBOMP need to know the number of targets or the noise power, which is commonly unknown in practice. For 3D-SL0 algorithm, several parameters need to be tuned, which requires a lot of efforts in some cases. Therefore, based on the principle of IAA algorithm [START_REF] Yardibi | Source localization and sensing: a nonparametric iterative adaptive approach based on weighted least squares[END_REF] and its efficient 2D version [START_REF] Jahromi | Two-dimensional iterative adaptive approach for sparse matrix solution[END_REF], we propose a nonparametric and user parameter-free algorithm, called T-IAA algorithm, to solve [START_REF] Yardibi | Source localization and sensing: a nonparametric iterative adaptive approach based on weighted least squares[END_REF].

Conventionally, to fit the classical CS theory,  is vectorized and estimated by the following 1D minimization problem. . The minimization of ( 16) yields [START_REF] Kolda | Tensor decomposition and applications[END_REF] Due to the unknown R, IAA algorithm operates iteratively to get [START_REF] Kolda | Tensor decomposition and applications[END_REF]. For most practical applications, convergence occurs after no more than 10-15 iterations [START_REF] Jahromi | Two-dimensional iterative adaptive approach for sparse matrix solution[END_REF]. The initialization of IAA algorithm is

     Q s ( 16
)
where b  is the b-th element of  , b  is the b-th column of  , 1 2 H 1 || || b b   Q x x Q x ,
H 1 H 1 = un b b b b       R s R
0 HH = un b b b b     s , and the o-th (o=1,2,…,O) iteration is carried out by H 1 H 1 H 1 H 1 = ( ) ( ) o o un o b b b b             s ( 18
)
where

1 1 2 {| | } o o diag     
and (• ) -1 denotes matrix inverse. Based on the Capon filter property, the vector form of ( 18) is given by [START_REF] Yardibi | Source localization and sensing: a nonparametric iterative adaptive approach based on weighted least squares[END_REF]

Input: un S , 1 un F , 2 un F , 3 un F , 1 o
 , and the iterations number T.

1) Initialization:

0 0 00 0, 0, , 0 un U R P      S , 2 0 0 || F R  
, and 1 t  ; 2) Update:

1 H H H 1 1 1 1 1 1 2 2 3 3 2 1 2 1 1 2 2 3 3 1 1 2 12 1 1 1 1 , [ ( ) ( ) ( ) ], , / , , , | | ,
U P R R T R                                      F F F F F F
3) Iteration:

1 tt , if tT  go back to 2), otherwise, stop.

Output:

1 oT   UU .

Since ()

oo vec   , the tensor based IAA algorithm can be derived as shown in Table . 2, in which the initialization step is

f x z 2R/c sin / c   s in / c   1 1 un Q I J C    S L M N C     1 1 Q L un C   F 1 2 I M un C   F 1 3 J N un C   F
obtained by making use of the minimum l2-norm solution and † ()  denotes pseudo-inversion. 1) Initialization: 3) Update:

0 † † † 1 1 2 2 3 3 ( ) ( ) ( ) un un un un      S F F F , 0 0 2 ||  ,
1 1 H H H 1 1 2 2 3 3 [ ( ) ( ) ( ) ]

IV. EXPERIMENT RESULTS

In this section, we present some experiment results to illustrate the effectiveness of the designed system and the proposed imaging algorithms. It should be pointed out that, to improve the imaging quality, the antenna direct coupling signal of the designed system is measured in an anechoic chamber and then subtracted from the real measurement data. The working frequency of the designed system is from 4.75 to 5.25 GHz with 201 frequency steps, giving a range resolution of 0.3 m. Due to the antenna size limitation, the inter-element spacing of the transmitting array and the receiving array are set to be 5 cm and 6 cm, respectively, giving an unambiguous angle about 60 degrees. The resolution in x direction is about 0.067 rad, and the resolution in z direction is about 0.080 rad. The point spread function of the designed cross-MIMO radar system for a target located at (10, 0, 0) can be obtained as shown in Fig. 4. Firstly, two trihedral corner reflectors (CRs) with size 40 cm and 30 cm are used as the point targets. They are located at (0, 6.5, 0) m and (1, 9.5, 0.5) m approximately. The experiment setup is shown in Fig. 5 and the focused images obtained by BP algorithm, 3D FPFA algorithm, 3D FPFA algorithm together with spatially variant apodization (SVA) [START_REF] Stankwitz | Nonlinear apodization for sidelobe control in SAR imagery[END_REF] approach, and T-IAA algorithm (for which we set O=T=10) with 1/8 data are shown in Fig. 6, where the coordinate is same as Fig. 2 and the dynamic range is 20 dB.

It can be seen from Fig. 6 (a) that two CRs can both be effectively imaged. Compared with BP algorithm, the proposed 3D FPFA algorithm can achieve similar result, as shown in Fig. 6 (b), but with much reduced computing time. Besides, by using the SVA method, the high-level sidelobes can be effectively suppressed without reducing the resolution, as shown in Fig. 6 (c). With 64 randomly selected pairs of transceivers and 101 frequencies, the imaging result of the proposed T-IAA algorithm is shown in Fig. 6 (d), from which it can be observed that the sidelobes are much reduced. Furthermore, measured by TIC and TOC instructions in MATLAB and averaged by 50 trials, the running time of BP algorithm, 3D FPFA algorithm, 3D FPFA algorithm with SVA approach, and T-IAA algorithm are 7.76 s, 0.12 s, 0.48 s, and 1.50 s, respectively. This indicates that BP algorithm is most time-consuming and the proposed algorithms can reduce the computing time. Besides, since only 1/8 data is used by the proposed T-IAA algorithm, the memory usage is also reduced. Apart from the proposed T-IAA algorithm, other T-CS algorithms can also be used to solve [START_REF] Yardibi | Source localization and sensing: a nonparametric iterative adaptive approach based on weighted least squares[END_REF]. For example, the imaging results obtained by Kron-OMP algorithm and NBOMP algorithm [START_REF] Caiafa | Multidimensional compressed sensing and their applications[END_REF] with 1/8 data are shown in Fig. 7. By comparing with Fig. 6 (d) and Fig. 7, it can be observed that the results of these three T-CS algorithms are similar, strong sidelobes can be effectively reduced. However, as mentioned in Section-III, OMP-type algorithms need to get the target number or the noise level to terminate the iterations, which is difficult to estimate. On the contrary, the T-IAA algorithm is nonparametric and can always converge within 10-15 iterations. Therefore, T-IAA algorithm maybe more suitable in practical applications.

Then, one of the CRs, which is mounted on a linear scanner, as indicated by the solid red circle in Fig. 5, is moved away from the radar system from 2 mm to 20 mm with a 2 mm step along the line of sight (LOS). The displacement estimation results obtained by BP algorithm, FPFA algorithm, and T-IAA algorithm are shown in Fig. 8. It can be learned that all these methods can achieve accurate displacement estimates. We note that, due to the influence of random selection, the joint T-IAA algorithm, which can be easily derived from the joint 2D IAA algorithm proposed in [START_REF] Feng | Jointly iterative adaptive approach based space time adaptive processing using MIMO Radar[END_REF], is used in this case.

At last, a building with four floors is used as the target, as shown in Fig. 9. The purpose of this experiment is to validate the imaging performance of the designed system for distributed target 3D imaging. The imaging result obtained by BP algorithm and 3D FPFA algorithm with SVA approach are shown in Fig. 10 with a dynamic range of 25 dB. For BP algorithm, strong sidelobes make the imaging result difficult to understand. However, by using SVA approach, four layers and main parts of the building can be clearly observed by 3D FPFA algorithm. This result further demonstrates the effectiveness of the designed C-band cross-MIMO radar system and the proposed imaging algorithm.

CONCLUSIONS

A new methodology for real time 3D imaging and displacement measurement based on GBSAR interferometry has been presented. It is based on a 2D cross-MIMO radar system which allows an effective 2D aperture synthesis and a fast data acquisition. The SAR focusing algorithm is based on a pseudo-polar spherical coordinate. Experiments have been carried out to prove the concept of the proposed methodology using a C-band prototype. The advantages of the proposed methodology in terms of focusing accuracy has been studied using trihedral corner reflectors (CRs) deployed in the scene. The proposed cross-MIMO radar has also been tested in an interferometric configuration to measure the displacements of CRs, using the proposed focusing algorithms, and comparing the results with those provided by other traditional focusing algorithms. Future work will be devoted to the application of the cross MIMO prototype and focusing algorithms to the monitoring of structures, glaciers and landslides.

Fig. 1 .

 1 Fig. 1. Configuration of the 3D imaging radar system.

.

  from the center of radar system to the target, θ0 denote the angle between op and the y-o-z plane (i.e. azimuth angle), and φ0 denote the angle between op and the x-o-y plane (i.e. elevation angle), we have Based on the cosine theorem, the distance between i-j-th sampling point and the target located at far field can be approximated as 0

F2,

  and F3 are the Fourier transform matrices.

unSF , 2 un F and 3 unF

 23 is the under-sampled received signal cube, 1 un are the under-sampled Fourier transform matrices.

Fig. 3 .

 3 Fig. 3. Compressive measurement of the target reflection coefficient tensor.

  denotes Kronecker product, and vec(• ) denotes the vectorization operation. IAA attempts to solve[START_REF] Jahromi | Two-dimensional iterative adaptive approach for sparse matrix solution[END_REF] by minimizing the following weighted least-squares cost function.

Input: un S , 1 unF , 2 unF , 3 unF

 123 , and the iterations number O and T.

1 oo

 1 , if oO  go back to 2), otherwise, stop.

Fig. 4 .

 4 Fig. 4. Point spread function of the designed system for a target at (10, 0, 0).

Fig. 5 .

 5 Fig. 5. Experiment with two trihedral corner reflectors as targets.

Fig. 6 .

 6 Focused images obtained by (a) BP algorithm, (b) FPFA algorithm, (c) FPFA algorithm with SVA, and (d) T-IAA algorithm with 1/8 data.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Focused images obtained by (left) Kron-OMP algorithm and (right) NBOMP algorithm with 1/8 data.

Fig. 9 .

 9 Fig. 9. Experiment with a four-floor building as target.

Fig. 10 .

 10 Fig. 10. Focused image of the building as distributed target obtained by (left) BP algorithm and (right) 3D FPFA algorithm with SVA approach.

  that can be obtained by the tensor version of CG algorithm, as shown in Table.1, (• ) * denotes conjugate, and denotes Hadamard product.

			o	1 H      1 H 1 () o o  	1 H    o su un o 	1 	(19)
	where	1  us is obtained by conjugate gradient 1 H 1 () o o un    
	(CG) algorithm to reduce the computation cost.
	Assuming	1 || 2 o    1  o	, it can be derived that
					1    { ( o diag vec	o	1 	)}	(20)
	where diag(⋅) denotes diagonalization.
	Therefore, based on	3 un  F 		2 F	un		1 un F and the properties
	of Kronecker product and Hadamard product, we can obtain
	(  o vec	11 [ ooun 1 1 ( )   U F	H	2 	(	2 F	un	)	H	3 	(	3 un F	H ) ])	(21)
	where	1 o U is the tensor form of	1 o

u

TABLE I .

 I TENSOR BASED CG ALGORITHM.

  / .

	t P R	t	t	t P T	o	P	t	un		un	un
	T T		un		un	un	tt			P	t	T
	U	tt	t	t	t	t	t	t	t	F	t	t	t

TABLE II .

 II TENSOR BASED IAA ALGORITHM.


Jean-Michel Friedt

Weike Feng

Stéphane Chrétien

Gwenhaël Merou

Motoyuki Sato

Gwenhael Goavec-Merou

Passive radar for measuring passive sensors: direct signal interference suppression on FPGA using orthogonal matching pursuit and stochastic gradient descent

Keywords: Passive sensor, surface acoustic wave, short range passive RADAR, orthogonal matching pursuit, stochastic gradient descent, interference suppression

INTRODUCTION

Passive RADAR provides an attractive solution to multiple remote sensing issues including not relying on a dedicated emitter to illuminate targets bathed in electromagnetic smog, [START_REF] Griffiths | An introduction to passive radar[END_REF] and being able to detect signals returned by a cooperative target whose signature is representative on a physical quantity (identification, temperature, strain ...). [START_REF] Feng | Passive radar delay and angle of arrival measurements of multiple acoustic delay lines used as passive sensors[END_REF] The basic signal processing principle of passive RADAR is to collect on the one hand a reference signal of the pseudo-random signal -most passive radar systems use existing commercial communication or navigation signals with appropriate random properties -illuminating the target, and on Further author information: send correspondence to J.-MF, E-mail: jmfriedt@femto-st.fr 1 the other hand the surveillance signal reflected by the target, possibly frequency shifted by a Doppler offset frequency if the target is moving 3 (Fig. 1). Unlike active RADAR, the emitted continuous signal is not known so that must be observed by a reference channel (or reconstructed when a single receiver is used [START_REF] Searle | Signal processing considerations for passive radar with a single receiver[END_REF] ) rather than locally memorized. The emitted signal might leak on the surveillance channel, inducing a strong background correlation signal that might be broad enough to overwhelm the signals reflected by the targets. The context of this study focuses on probing one particular case of targets designed as passive wireless sensors with echoes delay by 1 to 3 µs. A cooperative target illuminated by a non-cooperative emitter returns echoes whose delays are dependent not only on the bistatic range but also on some physical quantities of the sensor environment. Such a result is achieved by converting the incoming radiofrequency signal into an acoustic wave by patterning electrodes on a piezoelectric substrate in the classical Surface Acoustic Wave (SAW) transducer operation. [START_REF] Morgan | Surface acoustic wave filters: With applications to electronic communications and signal processing[END_REF] By a differential measurement in which multiple echoes induced by the acoustic transducer and returned by the sensor are sampled, subtracting delays between multiple echoes gets rid of the bistatic range delay dependence and provides means for remote sensing physical quantities using passive, cooperative targets probed by a passive RADAR receiver. The design of the SAW transducer aims at delaying echoes beyond clutter and environmental reflections: considering the receiver sensitivity, delaying echoes by 1 to 2 µs (150 to 300 m monostatic range) guarantees that the sensor response is well beyond any non-cooperative target reflection (Fig. 2).

Since the classical estimator for detecting echoes in a RADAR signal is the correlation, passive RADAR target identification heavily relies on correlating an unknown but measured reference signal from a noncooperative illuminator with the surveillance signal collected by an antenna (or antenna array in our case) directed towards the sensor. However, any structure in the illuminating signal as often found in communication signals will depart from the random structure whose correlation would be a Dirac function and will introduce some long term memory effects with strong correlation side lobes. Thus, even if clutter has faded out by the time the sensor echoes are detected, the illuminating signal source autocorrelation might still overwhelm the sensor signal as will be demonstrated in the case of WiFi. This Direct Signal Interference (DSI) removal is the core requirement to extract the target (here cooperative with the sensor designed on purpose) signature from the correlation between the reference signal and the cleaned surveillance signal.

The classical approach to suppress this DSI signal is to use a least square weighing approach where time shifted copies of the reference signal are sought in the surveillance signal. [START_REF] Garry | Evaluation of direct signal suppression for passive radar[END_REF] Assuming the echoes generated by targets -whether non-cooperative or cooperative sensor targets for which the echoes are designed to lie beyond clutter -are delayed by more than a minimum time offset T , all copies of the reference signal delayed between 0 and T in the surveillance signal with estimated weights are subtracted before cross-correlating the reference and the surveillance signal to identify target. Practically, this operation is achieved by creating a matrix X with time delayed copies of the reference signal ref and identifying the weight w of the time-delayed copies of ref in the surveillance signal sur as

X =       ref 1 0 0 ... 0 ref 2 ref 1 0 ... 0 ref 3 ref 2 ref 1 ... 0 ... ... ... ... ... ref N ref N -1 ref N -2 ... ref N -T       : w = pinv(X) • sur (1)
where the pseudo-inverse pinv(X) is (X H • X) -1 • X H with H the conjugate transpose of the matrix. In this definition of X, signals before the first sample (ref 1 ) are assumed to be zero. A typical case is N large (numerous samples to improve signal to noise ratio) and T small (e.g. in a cooperative target short range RADAR system sampling the baseband at 100 MHz, removing clutter delayed by up to T = 200 ns (30 m monostatic range) is achieved with discrete delays from 0 to 20). The challenge of practically achieving this computation (Eq. 1) is that X H • X is an (T + 1) × (T + 1) matrix which must be inverted before being multiplied by the (T + 1) × N X H matrix, resulting in a (T + 1) × N matrix multiplied by the N element long surveillance signal to yield the T + 1 weights w of the contribution of the delayed reference signal in the surveillance (Fig. 3). Subtracting the delayed reference from the surveillance is finally achieved with sur -X • w. For example, 7 uses a Cholesky decomposition to compute the inverse, and hence does not benefit from the history of the matrix assembly.

One initial approach to reduce the computational burden is to reduce the number N of samples by selecting a random subset of N lines: the resulting under-sampled copy of X is named X . This number of sub-samples is dependent on the targeted accuracy and rises as the logarithm of the number of delays involved in the DSI. [START_REF] Kane | Sparser Johnson-Lindenstrauss transforms[END_REF][START_REF] Donoho | Compressed sensing[END_REF] In the experimental dataset shown in Fig. 3, only 3 time-delay copies (or 93 ns considering our sampling rate of 125/4=31.25 MS/s) of the reference signal need to be removed from the surveillance signal for the sensor to become dominant in the correlation between the reference and cleaned surveillance signal, while considering additional delays (Fig. 3, bottom chart) does not improve the sensor signal. Hence, only very few weights need to be computed to cancel DSI and allow for the sensor response to be analyzed.

Fig. 4 illustrates the strategy followed in this paper based on the DSI matrix representing the time-delayed copies of the reference signal sought in the surveillance signal. Here only positive delays are involved, but the practical algorithms consider positive as well as negative delays to account for correlation in the non-random reference signal. On the one hand a measurement set of the surveillance and reference channels has been collected and time-delayed copies of the reference are sought in the surveillance: the number of samples (lines) is fixed and the Orthogonal Matching Pursuit (OMP) algorithm iterates along the columns to sequentially identify and subtract the weighted copy of each time-delayed reference. On the other hand, the number of delayed copies of the refer- ence signal is set and as more data are collected, the weights associated with each time-delayed copy of the reference are updated with the new observations: the Stochastic Gradient Descent (SGD) operates along the lines for a fixed number of delays (columns). Both algorithms aim at removing the computation burden of the matrix pseudo-inverse classically used to identify the weight of each time-delayed copy of the reference in the surveillance signal.

observations: SGD         delays: OMP ---------------------→      x 1 0 • • • 0 x 2 x 1 • • • 0 . . . . . . . . . . . . x m x m-1 • • • x m-T     
This approach was used in 10 to demonstrate passive RADAR measurement of passive cooperative targets illuminated by a WiFi signal. However, the measurement update rate is limited by DSI suppression: the data acquisition, frequency transposition and filtering are processed in a pipeline running in the Field Programmable Gate Array (FPGA) of the Zynq system on chip (Fig. 1). Due to the computational burden of DSI suppression, the current implementation transfers the resulting reference and surveillance vectors to the host computer in charge of running the DSI suppression and cross-correlation. Since cross-correlation is readily achieved as a pipelined algorithm on an FPGA, running the full processing algorithm on FPGA for increased data acquisition rate requires solving the DSI suppression following an algorithm compatible with a pipelined approach applicable to a gate array.

ORTHOGONAL MATCHING PURSUIT FOR DIRECT SIGNAL INTERFERENCE REMOVAL

The computational burden of inverting the matrix remains high since it requires collecting the whole dataset before performing the DSI suppression. Here we consider the iterative subtraction process of subtracting sequentially the heaviest contributions of the reference signal in the surveillance channel, using the OMP algorithm. [START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF] Based on compressive sensing (CS) theory, [START_REF] Donoho | Compressed sensing[END_REF] the DSI can be estimated by OMP algorithm with under-sampled measurement, which can help to reduce the computational complexity. Indeed, OMP implementation on FPGA has been demonstrated already [START_REF] Ge | An efficient FPGA implementation of orthogonal matching pursuit with square-root-free qr decomposition[END_REF][START_REF] Morales-Pérez | FPGA implementation of orthogonal matching pursuit algorithm[END_REF][START_REF] Polat | High-speed FPGA implementation of orthogonal matching pursuit for compressive sensing signal reconstruction[END_REF] with reasonable computational power requirement with respect to those provided by the Zynq 7010 system on chip fitted on the Redpitaya platform we are using on our experiment.

In this algorithm, the contribution of each time delayed copy of X in the under-sampled surveillance signal sur is computed as (X ) H • sur (Fig. 5). The time delay with highest contribution is selected at index p 1 and its weight is computed as pinv(X (p 1 )) • sur , where X (p 1 ) is the p 1 th column of X (appendix, first listing). Notice that in this first iteration, p 1 is a scalar so that X (p 1 ) is a vector and

pinv(X (p 1 )) • sur = < sur , X (p 1 ) > ||X (p 1 )|| 2 2 ∈ C (2)
with <, > the dot product n sur n • X (p 1 ) * n since pinv(X (p 1 )) = (X (p 1 ) H • X (p 1 )) -1 • X (p 1 ) H and X (p 1 ) H • X (p 1 ) = ||X (p 1 )|| [START_REF] Feng | Passive radar delay and angle of arrival measurements of multiple acoustic delay lines used as passive sensors[END_REF] 2 . In the next step of the OMP algorithm, the surveillance signal cleaned from this initial contribution of the reference signal sur 1 = sur -X (p 1 ) • (pinv(X (p 1 )) • sur ) is then again analyzed for finding the next heaviest contribution of X in sur 1 . This requires searching for the next heaviest contribution p 2 of X in sur 1 as the column for which (X ) H • sur 1 is maximized. It is known that subtracting only column p 2 of X from sur 1 will only converge to cleaning sur from time delayed copies of the reference if p 1 and p 2 are different enough for the contributions to be independent. Practically, we expect the strongest contributions of the reference signal in the surveillance to lie at time delays close to 0 and not to be independent from each other. Hence, a global weight computation must be addressed by aggregating p 1 and p 2 in a vector p and identifying the optimal weight of both columns as pinv(X (p)) • sur . This computation will grow in complexity towards the initial problem since p will grow as more contributions of the reference signal are sought in the surveillance. However in this iterative approach, pinv(X (p)) can be deduced from the previous iteration by considering the block-wise inverse following the Sherman-Morrison formula:

(X (p) H • X (p)) -1 = X (p 1 ) H X (p 2 ) H • X (p 1 ) X (p 2 ) -1 = X (p 1 ) H • X (p 1 ) X (p 1 ) H • X (p 2 ) X (p 2 ) H • X (p 1 ) X (p 2 ) H • X (p 2 ) -1 (3)

Block matrix inversion

Since the new matrix inverse results from a previous iteration with the addition of one more vector from the reference signal to be subtracted from the surveillance, the problem of finding this new inverse can be expressed with the following shaped matrices:

A B B H C -1 = A -1 + 1 k A -1 BB H A -1 -1 k A -1 B -1 k B H A -1 1 k (4)
where

k = C -B H A -1 B.
Hence, in this equation A -1 is already known from the previous iteration, and only matrix products remain, no inversion (appendix, second listing).

The next section will compare both the least square estimate of the weights of the time-delayed copies of the reference signal in the surveillance signal using the pseudo-inverse (appendix listing 1) and the result of iterative inverse matrix computation (appendix listing 2) in which no explicit matrix inversion is needed but only matrix products. The data have been collected from an experimental WiFi passive RADAR setup with an 8-path antenna array receiver so that eight independent measurements of the the surveillance signals are analyzed. [START_REF] Feng | Passive radar delay and angle of arrival measurements of multiple acoustic delay lines used as passive sensors[END_REF] The reference signal will either be used to synthesize known surveillance signals for algorithm assessment, or will be compared with experimental measurements of surveillance signals for sensor characteristics extraction.

Experimental demonstration

Using OMP and the iterative matrix inversion implementation, the DSI is indeed suppressed as shown in Fig. 6. In this figure, the result of DSI suppression using the pseudo-inverse (top) is favorably compared with the OMP result (middle) as emphasized by the difference between the two curves (bottom). Thus, the iterative matrix inversion implementation of OMP is functional as an algorithm for extracting sensor echoes from the correlation between the DSI-cleaned surveillance signal and the reference signal. The implementation remains computationally intensive and requires keeping a memory of the whole reference and surveillance signals, using most resources on the low-end FPGAs (e.g. Zynq-7010 fitted on the Redpitaya board) considered for assembling an embedded reader with real time sensor echo detection capability.

STOCHASTIC ITERATIVE GRADIENT DESCENT WITH HARD THRESHOLDING

The OMP approach aims at identifying which time-delayed copy of the reference signal lies in the surveillance signal through an iterative process on the time delays, removing from the strongest to the weakest copy of for time delays ranging from -16 to +16 sampling intervals (-0.5 to +0.5 µs) for all 8 antenna measurements. The difference between the two results is shown as the bottom chart. In all cases the echoes of the sensors are visible on all 8-antenna signals between 1 and 3 µs, with some leftover DSI in the 0.5 to 1 µs delay range that does not impact our ability to recover the sensor measurement. The correlation peak at 3.2 µs is intrinsic to the WiFi signal structure and not part of the sensor response.

the reference signal in the surveillance signal. Thus, OMP assumes that the whole surveillance and reference datasets have been collected and stored in memory. Such an approach is poorly suited to FPGA implementation which is best designed for pipelined processing of new samples as they are collected. Furthermore, the iterative matrix inversion algorithm depicted earlier, despite removing the need for matrix inversion, still requires computing the scalar inverse 1/k, a challenging operating on an FPGA.

Thus, such an approach does not address the practical measurement issues that 1/ the time delays and weights of the reference signal copies in the surveillance signal hardly vary from one measurement to another considering the fast (radiofrequency) sampling rate and 2/ data from the surveillance and reference channels are continuously streamed from the two analog-to-digital converter channels. Hence, as shown in Fig. 4, rather than addressing the DSI suppression issue along the time-delay (matrix column) direction, we consider the challenge of identifying time-delayed copies of the reference signal in the surveillance signal along the acquisition (matrix lines) dimension.

The SGD method [START_REF] Frotzscher | A stochastic gradient LMS algorithm for digital compensation of Tx leakage in Zero-IF-Receivers[END_REF] aims at solving the least-squares problem using an iterative optimization scheme which avoids matrix inversions and can proceed incrementally, yielding a highly scalable algorithm. In our setting, identifying the weights of the time delayed copies of the reference signal in the surveillance signal can be rephrased as a least squares problem min

w∈C T +1 1 2 sur -Xw 2 2 (5)
The SGD algorithm [START_REF] Bertsekas | Incremental gradient, subgradient, and proximal methods for convex optimization: A survey[END_REF] performs gradient steps by updating the weight vector w of the form

w (l+1) = w (l) + α (l) X H i (sur i -X t i w (l) ) (6)
where i is randomly selected among the observed time indices. Standard choices for the learning rate α (l) include the constant step strategy and the α (l) = 1/ √ l schedule. Throughout these investigations we have selected α ranging from 5 • 10 -3 to 3 • 10 -2 on synthetic data (Fig. 7) and α = 10 -2 on experimental datasets (Figs. 8 and9).

Our new approach to the problem combines the stochastic gradient scheme (6) with the Iterative Hard Thresholding Scheme of 17 by setting to 0 the weights below a given threshold defined as a fraction of the maximum of the weight values to accelerate the convergence of the algorithm.

w (l+1) = T λ (l) w (l) + α (l) X H i (sur i -X t i w (l) ) (7)
where T λ is the thresholding operator that sets all values below the threshold λ to zero and λ (l) is chosen so as to select the values below 1% of the largest in absolute value.

Not only does the thresholded stochastic gradient descent method provide the weights of the DSI coefficients for future subtraction, it also allows for identifying the sensor echoes if a long enough horizon is allowed for the algorithm search. The assumption for stochastic gradient descent though is the sparsity of the weights, assuming few echoes in the returned signal (a well met hypothesis) spaced apart (an hypothesis to be discussed). . Left: stochastic gradient descend (SGD) applied to the synthetic dataset described in the text, with a hard threshold of λ = 0.05 times the maximum weight value. The red and blue circles indicate the complex weights applied to the delayed copies of the reference signal to synthesize the surveillance signal, emphasizing the excellent match between the targeted value and the weights identified by SGD. Right: evolution of the relative error (logarithmic scale) as a function of iteration number (abscissa -logarithmic scale) for various values of α. The iterations were stopped when the relative error variation reached 10 -8 , indicating no evolution of the weights during the gradient descent iterations.

Fig. 7 demonstrates the application of the stochastic gradient method (appendix, listing 3: note the real implementation of the proposed method for complex values) to synthetic surveillance data generated by adding time delayed copies of experimental WiFi reference signals. In these examples, the time delays introduced were -60, -30, 0, 50 and 70 with arbitrary complex weights of representative of typical SAW echo delays since at a rate of 125/4=31.25 MS/s, an echo delayed by 1 to 2 µs is located at sample numbers 32 to 64. Fig. 7 (right) emphasizes the dependency of the convergence rate on the free parameter α which is here selected arbitrarily to yield convergence as quickly as possible: a rationale approach [START_REF] Chrétien | Hedging parameter selection for basis pursuit[END_REF] to identifying α will be part of further investigations.

Hence, in the case of the iteratively thresholded stochastic gradient method, the objective might no longer DSI suppression to allow for the cross correlation, but sparse distributions of delayed copies of the reference signal in the surveillance signal position and weight identification. However, we have observed that sensor echoes are too weak to be identified by the SGD algorithm and are cancelled by the hard thresholding due to their low power with respect to the DSI. Indeed with a echoes 25 to 50 times weaker than DSI (Fig. 3), either the threshold λ is too high and cancels the weight of the sensor echo copies of the reference signal in the surveillance, or λ is too small at the level of the noise and prevents SGD from converging. Hence, the demonstration on experimental data remains on DSI weight identification (Fig. 8) followed by DSI subtraction and cross-correlation of the cleaned surveillance signal with the reference signal (Fig. 9) for sensor echo extraction. Fig. 8 demonstrates convergence of the SGD algorithm in about 2000 iterations. At a sampling rate of 125/4=31.25 MS/s or 32 ns/cycle, convergence will require 64 µs of data acquisition, well below typical time constants for temperature or chemical sensing using passive wireless SAW sensors. Furthermore, for faster physical quantities such as strain sensing, 19 the sensor environment is expected not to vary significantly over sub-100 µs durations and the DSI identification using SGD considered to be still valid for extracting sensor echoes whose phase is representative of stress levels in the sensor. Fig. 9 (left) compares the cross-correlation of the cleaned surveillance signal with the reference signal either using the least-square solution (LS) obtained by applying the pseudo-inverse of the time-delayed copies of the reference signal matrix to the surveillance signal, or using the SGD algorithm. While both results exhibit excellent match, the SGD algorithm is well suited for implementation in an FPGA. In all these demonstrations, the classical random selection of the sample in the dataset when applying SGD is replaced with a sequential selection of each element of the dataset (cyclic buffer with the sample index selected modulo the dataset length -32768 samples in these examples) representative of a continuous stream of data collected by the analog to digital converters connected to the surveillance and reference channels of the passive RADAR receiver. Running a C++ implementation of 2000 iterations of the SGD algorithm on the Cortex-A9 ARM processing system of the Zynq-7010 fitted on the Redpitaya board lasts 33.5 ms (168 ms for 10000 iterations), emphasizing the need for an FPGA implementation to achieve the 500-fold speed gain targeted for real time processing (Fig. 9, right).

CONCLUSION

We have addressed the issue of real time direct signal interference removal in a passive RADAR setup designed to measure passive sensor responses. Since radiofrequency datastreams are collected using a field programmable gate array, algorithms compatible with such architectures are considered, namely with minimal memory consumption and lack of matrix inversion. We have demonstrated that an iterative matrix inversion applied to Orthogonal Matching Pursuit will remove the need for matrix inversion in the FPGA and restrict computations to vector multiplications. However, a scalar value inverse remains to be computed, a potentially challenging task in the FPGA hardly supporting division. As an alternative, we have explored the Stochastic Gradient Descent algorithm with Hard Thresholding which demonstrates comparable performance with an algorithm suitable for processing streams of data and adapting continuously the weights of the Direct Signal Interference for subtraction and correlating the cleaned surveillance signal with the reference signal.

Figure 1 .

 1 Figure 1. Left: data acquisition and processing principle: a Redpitaya board based Zynq System on Chip collects data from 125 MS/s analog to digital converters (ADC) following a radiofrequency to baseband frequency transposition (hardware mixer and filters). The sampled data are frequency transposed and filtered for adjacent WiFi channel data collection by the Programmable Logic (PL) part of the Zynq. At the moment, the resulting data are transferred to the Processing System (PS) to be sent to the host computer for further processing, limiting the measurement refresh rate. The present work aims (red) at including the DSI suppression and correlation in the PL. Right: experimental setup using an 8-patch antenna array as receiver for direction of arrival computation, making the DSI suppression load even heavier when probing multiple sensors. The antenna array (background) to sensor (foreground) distance is 1 m. The inset on the left is a screenshot of the computer display on the right exhibiting the azimuth-range map with the sensor echo delays (representative of the physical quantity) along the verical axis and the azimuth for spatially separating the sensor response along the horizontal axis.

Figure 2 .

 2 Figure 2. Surface Acoustic Wave (SAW) reflective delay line scattering coefficient response, bottom in the frequency domain and top in the time domain as computed through the inverse Fourier transform of the frequency domain characterization. Echoes are delayed by 1 to 2.5 µs beyond clutter of short range RADAR readers, while the actual delays encode an 8-bit identifier. The precise time delays between echoes informs on the sensor temperature through the acoustic wave velocity. The narrow bandwidth (35 MHz) was selected on purpose as representative of the passive RADAR bandwidth despite preventing the separation of two of the echoes.

Figure 3 .

 3 Figure 3. Passive RADAR interrogation of the SAW sensor characterized in Fig.2using a WiFi non-cooperative source, with the removal, from top to bottom, of none or 1 (sensor signal hidden under the sidelobes of DSI), 2 (sensor signal becoming visible between 1 and 3 µs), 3 (rejection of DSI and only sensor response left) and 4 (no improvement with respect to 3-sampling time delay subtraction) sampling delay copies of the normalized reference signal from the normalized surveillance signal. Such few delayed copies is consistent with the systematic analysis shown in Fig.5. Note that, in the top three figures, the peaks at about 3.2 µs are generated by the properties of WiFi signal.

Figure 4 :

 4 Figure 4: Direct Signal Interference matrix shaped by delaying the reference signal and search for the T weights (one for each column) of each copy in the surveillance signal.

Figure 5 .

 5 Figure 5. Experimental dataset analysis exhibiting the normalized weights of the delayed copies (delay in sample number as abscissa) of the reference in the surveillance channel. Notice that, as expected, weights are negligible in the negative delay and sharply drop to negligible values in the positive delays above 5 timesteps, consistent with Fig.3.

Figure 6 .

 6 Figure 6. Comparison of the weight estimated through the pseudo inverse function (top) and through OMP (middle)for time delays ranging from -16 to +16 sampling intervals (-0.5 to +0.5 µs) for all 8 antenna measurements. The difference between the two results is shown as the bottom chart. In all cases the echoes of the sensors are visible on all 8-antenna signals between 1 and 3 µs, with some leftover DSI in the 0.5 to 1 µs delay range that does not impact our ability to recover the sensor measurement. The correlation peak at 3.2 µs is intrinsic to the WiFi signal structure and not part of the sensor response.

Figure 7

 7 Figure 7. Left: stochastic gradient descend (SGD) applied to the synthetic dataset described in the text, with a hard threshold of λ = 0.05 times the maximum weight value. The red and blue circles indicate the complex weights applied to the delayed copies of the reference signal to synthesize the surveillance signal, emphasizing the excellent match between the targeted value and the weights identified by SGD. Right: evolution of the relative error (logarithmic scale) as a function of iteration number (abscissa -logarithmic scale) for various values of α. The iterations were stopped when the relative error variation reached 10 -8 , indicating no evolution of the weights during the gradient descent iterations.

Figure 8 .

 8 Figure 8. Stochastic Gradient Descent (SGD) applied to experimental WiFi passive RADAR signals collected as a SAW sensor was illuminated by a non-cooperative emitter. Left: the weights of DSI identified by SGD match the least-square (LS) analysis found by applying the pseudo-inverse of the time delayed copies matrix to the surveillance signal. Right: relative error (logarithmic scale) as a function of iteration number (logarithmic scale) demonstrating convergence in about 2000 iterations.

13 Figure 9 .

 139 Figure 9. Left: comparison of the normalized cross-correlation of the raw surveillance signal and reference signal in which the sensor response is hidden below interference and WiFi autocorrelation structure, surveillance signal cleaned from DSI using the least square (LS, pseudo-inverse matrix calculation) classical technique and the proposed Thresholded SGD algorithm. Excellent match is observed on the sensor echoes using the latter techniques. Right: FPGA SGD with Hard Thresholding implementation principle. Colours refer to the line numbers in Listing 3 in the appendix. Iteration control logic is not shown on this chart focusing on arithmetic processing.

ACKNOWLEDGMENTS

This work is supported by the OscillatorIMP ANR (French Agence Nationale de la Recherche) Equipex grant as well as the FAST-LAB Labcom ANR grant. J.-MF is grateful to Tohoku University (Sendai, Japan) for a visiting scientist grant during winter 2017 which prompted this investigation.

Appendix

OMP algorithm implementation using GNU/Octave's pinv() pseudo-inverse function:

Listing 1. Explicit matrix pseudo-inverse computation least square method function [x,Ind]=omp(PSI_und,y,iter_max) % PSI is the pre-defined under-sampled dictionary (delayed copies of the reference signal) % y is the under-sampled received signal % x is the coefficients of the selected column of PSI % Ind is the index of the selected column of PSI SGD method with Iterative Hard Thresholding (IHT):

Listing 3. Stochastic Gradient Descent implementation for complex datasets % assumes vector x to be reference channel and vector y to be surveillance --returns w1 weights for k = 1:number_of_iterations % typically 1000-10000 i = mod(k,N-2*p-1)+1; % should be random, but systematic iteration simulating datastream i = i+p; % allow for negative delays xx = [real(x(i-p:i+p)).',-imag(x(i-p:i+p)).';... imag(x(i-p:i+p)).',real(x(i-p:i+p)).'];% real and imaginay e = [real(y(i));imag(y(i))] -xx*w1; % real and imaginary, since weights w1 are complex w1 = w1 + alpha * xx.' * e; % gradient descent w1_real = w1(1:end/2); w1_imag = w1(end/2+1:end); w1_com = w1_real +1j * w1_imag; I = find(abs(w1_com)<threshold*max(abs(w1_com))); % hard thresholding @ threshold*max(weights) w1_com(I) = 0; w1 = [real(w1_com);imag(w1_com)]; end