A data driven compressive sensing approach for time-frequency signal enhancement

Signal Processing 141 (2017) 229–239

Contents lists available at ScienceDirect

Signal Processing journal homepage: www.elsevier.com/locate/sigpro

A data driven compressive sensing approach for time-frequency signal enhancement Ivan Volaric a, Victor Sucic a,∗, Srdjan Stankovic b a b

Faculty of Engineering, University of Rijeka, Vukovarska 58, HR-51000 Rijeka, Croatia Faculty of Electrical Engineering, University of Montenegro, Dzordza Vasingtona bb, CG-81000 Podgorica, Montenegro

a r t i c l e

i n f o

Article history: Received 16 December 2016 Revised 15 May 2017 Accepted 5 June 2017 Available online 13 June 2017 Keywords: Time-frequency representation Ambiguity function Signal sparsity Compressive sensing Basis pursuit Linear unconstrained optimization

a b s t r a c t Signals with the time-varying frequency content are generally well represented in the joint timefrequency domain; however, the most commonly used methods for time-frequency distributions (TFDs) calculation generate unwanted artifacts, making the TFDs interpretation more difficult. This downside can be circumvented by compressive sensing (CS) of the signal ambiguity function (AF), followed by the TFD reconstruction based on the sparsity constraint. The most critical step in this approach is a proper CS-AF area selection, with the CS-AF size and shape being generally chosen experimentally, hence decreasing the overall reliability of the method. In this paper, we propose a method for an automatic data driven CS-AF area selection, which removes the need for the user input. The AF samples picked by the hereproposed algorithm ensure the optimal amount of data for the sparse TFD reconstruction, resulting in higher TFD concentration and faster sparse reconstruction algorithm convergence, as shown on examples of both synthetical and real-life signals. © 2017 Elsevier B.V. All rights reserved.

1. Introduction In many real-life applications signal frequency content is a key information which needs to be extracted. In order to do this, one can simply use the Fourier transform, but in doing so, signal time attributes are lost. However, by using time-frequency distributions (TFD) one can analyze the evolution of signal energy as a function of both time and frequency, providing additional information about the signal nature. However, if signal has more than one linear frequency modulated (LFM) component, or a non-LFM component, its TFD gets corrupted by highly oscillatory artifacts. Quadratic TFDs (QTFD), most commonly used TFDs in practice, utilise filtering with 2D low-pass filters, which inherently affects time-frequency (TF) localization properties [1]. The need for a trade-off between interference suppression and TF localization has led to a number of TF localization improvement methods, one of which is described in the sequel. Over the last few years, compressive sensing (CS) has been an important research topic [2–5], with applications in medicine [6,7], geophysics [8,9], communication [10,11], etc. Traditionally, compressive sensing (CS) implies signal sampling with sub-Nyquist frequencies, with samples randomly picked [2–4,6]. However, the

∗

Corresponding author. E-mail address: [email protected] (V. Sucic).

http://dx.doi.org/10.1016/j.sigpro.2017.06.013 0165-1684/© 2017 Elsevier B.V. All rights reserved.

samples can be picked to favour specific signal features, while discarding the others [12–14]. This is followed by a signal reconstruction algorithm for solving unconstrained optimization problems (i.e. basis pursuit (BP) [15–22], matching pursuit (MP) [23], orthogonal matching pursuit (OMP) [24,25], etc.). However, the reconstruction algorithm leads to a meaningful result only if the signal is sparse, which means that the signal can be represented in a certain domain with K non-zero coefficients, where K  Nt (Nt being the number of signal samples in the time domain). Most signals are non-sparse in the domain of interest, but can become sparse (or approximately sparse) by applying a domain transformation. For example, a sinusoidal signal can be represented with only one sample in the frequency domain. Ideal TFDs are sparse since they are composed of components instantaneous frequency (IF) trajectories, hence CS can be utilized in such a way to include only the signal components samples, while discarding the interference samples; by applying a reconstruction algorithm, resolution loss is minimised [7,8,13,14,25–30]. Current TF signal processing approaches have focused on various optimization algorithms for signal reconstruction, leaving CS area size and shape selection underutilized. The CS area is usually a rectangle, containing approximately Nt samples inside of it [13,26,27]. In this paper, we propose a method for data driven automatic CS area selection. Our goal is to select CS area as large as possible without artifact inclusion, which then increases the signal

230

I. Volaric et al. / Signal Processing 141 (2017) 229–239

reconstruction algorithm input data amount, hence decreasing its computational requirements. The paper is organized as follows. Section 2 gives a short introduction to TFDs, while Section 3 describes TFDs as a sparsity inducing signal representation, and introduces the CS area selection algorithm. In Section 4 we compare different optimization algorithms performances based on the measure of reconstructed TFD concentration [31] for the case when the CS area is selected both manually and when it is selected automatically by the proposed method. 2. Quadratic time-frequency distributions Let us consider an LFM signal z(t) with a time-varying phase

ϕ (t), and a slowly varying amplitude A(t) of the form: z(t ) = A(t )e jϕ (t ) .

(1)

Its ideal TFD, ρz (t, ω ), is a set of Dirac functions, with a perfect energy localization around the signals instantaneous frequency, ω0 (t):

(2)

with ω0 (t ) = dϕ (t )/dt. The ideal TFD in most cases is impossible to accomplish, since practical TFDs are not perfectly localized, and furthermore, they are corrupted by the cross-terms. The WignerVille distribution (WVD), Wz (t, ω), provides perfect localization for a single LFM component [1]:



∞

−∞

Rz (t, τ )e− jωτ dτ ,

Rz (t, τ ) = z t +

τ  ∗

τ

2

2

t−

z

Rz (t, τ ) =



zi t +

i=1

⎡

.

(4)

τ  ∗

τ

2

2

zi t −

⎤ N N     c c ⎢ τ  ∗ τ ⎥ + zj t − ⎣zi t + 2 ⎦, 2 i=1



(5)

j=1 j=i

Rz (t, τ )e− jν t dt −∞  ∞ ∞ = Wz (t , ω )e− jν t e jωτ dt dω,

(6)

−∞

the cross-terms can be filtered out with a 2D low-pass filter, since they are highly oscillatory, and thus located away from the AF plane origin [1]. However, in doing so, the auto-terms get partially filtered out as well, thus reducing TF concentration of the components. This has led to the formulation of QTFD as [1]:

A z ( ν , τ ) = A z ( ν , τ ) g( ν , τ ) ,

(7a)

ρz (t, ω ) = Wz (t, ω ) ∗ ω∗ γ (t, ω ),

(7b)

t

|ν| < D, |τ | < E, otherwise,

(8)

where the parameter c defines the shape of the kernel, while the parameters D and E specify the spread of the kernel along the respective AF axis. Since there is no single best performing kernel for all signals, the need to adaptively construct kernel has arisen. One of the best known methods is the radially Gaussian kernel (RGK) [33], in which the kernel is obtained by solving the optimization problem with the following objective function:

gopt (ν, τ ) = arg max



∞



∞ −∞

| A z ( ν , τ ) g( ν , τ ) | 2 d τ d ν ,

(9)

with the constraints which ensure that the kernel has properties of the low-pass filter. Since cross-terms do not contain signal energy, maximizing (9) ensures that the auto-term energy is transferred from the AF to the resulting TFD, without attenuation, while avoiding the cross-terms.

3.1. Compressed sensed ambiguity function As mention earlier, ideal TFDs are inherently sparse, since they are composed of the components IFs, thus requiring only Nc Nt <
Az (ν, τ ) = φ (ν, τ )  Az (ν, τ ),

φ (ν , τ ) =

1, 0,

( ν , τ ) ∈ , otherwise.

(10)

(11)

The sensing matrix defined in this way discards highly oscillatory cross-terms located away from the domain origin, while preserving the auto-terms located closer to the origin. In the standard CS notation, matrix multiplication is commonly used to define connection between the observation and solution matrices, thus from (6) and (10):

Az (ν, τ ) = ψ · ϑz (t, ω ),

∞

−∞

cE 2

where the operator  denotes element-by-element matrix multiplication, and φ (ν , τ ) is the sensing matrix which defines Nν × Nτ area around the AF plane origin:

where the factors under the first sum are the auto-terms, and the remaining factors are the cross-terms. As it can be seen from (5), the cross-terms are mathematical byproduct introduced by LAF’s quadratic nature, and will appear midway between each pair of components [1]. In the ambiguity function (AF), Az (ν , τ ), defined as:

Az ( ν , τ ) =

cD2

e2c e ν 2 −D2 e τ 2 −E 2 0,

3. Time-frequency distribution as a sparsity inducing domain

However, if the signal has Nc > 1 components, its LAF becomes: Nc 

g( ν , τ ) =

(3)

where Rz (t, τ ) is the signal localized autocorrelation function (LAF), defined as:



ω

t

frequency, respectively. One of the state-of-the-art TFD is the compact kernel distribution (CKD) defined as [32]:

g(ν,τ ) −∞

ρz (t, ω ) = A2 (t )δ (ω − ω0 (t ) ),

Wz (t, ω ) =

where g(ν , τ ) and γ (t, ω) are filter functions (also known as kernels), while Az (ν, τ ) and ρ z (t, ω) are the filtered AF and TFD, respectively. The symbols ∗ and ∗ denote convolution in time and

(12)

where ψ is a domain transformation matrix, representing a twodimensional Fourier transform, and ϑz (t, ω) is a sparse TFD reconstructed from the CS-AF. However, in most practical realizations, ψ and its inverse are implemented as functions, since matrix free algorithms are significantly faster and require less memory space [17,20]. Since Az (ν, τ ) has cardinality card(Az (ν, τ )) = Nν · Nτ samples, and ϑz (t, ω) has cardinality of card(ϑz (t, ω )) = Nt · Nω samples, where card(Az (ν, τ ))  card(ϑz (t, ω )), the system is underdetermined, thus ϑz (t, ω) can have an infinite number of possible solutions, and the goal of the reconstruction algorithm is to find the optimal solution to:

ϑz (t, ω ) = ψ H · Az (ν, τ ) = ψ H · φ (ν, τ )  Az (ν, τ ),

(13)

I. Volaric et al. / Signal Processing 141 (2017) 229–239

where the operator · H denotes Hermitian transpose, while matrices ψ H and φ (ν , τ ) have the required low mutual coherence [4]. The mutual coherence is measured as the largest correlation between any two elements of sensing and transformation matrices:

μ (φ , ψ ) =



Nt Nω · max |φk , ψ j | ∈ [1,



Nt Nω ].

(14)

If we take samples from a sparse domain, there is high probability that most of them would be equal to zero, which would give very little information about the signal energy. Instead, we want to take samples from a such sensing domain, which when transformed back to the sparse domain, will spread out through entire measurement domain, giving us as much information about the signal energy as possible. It is well known that Fourier basis is indeed incoherent with spike basis [4,5]. 3.2. Compressive sensing area selection Proper selection of Nτ and Nν is crucial for TFD localization and cross-term suppression. Too few samples will reduce TF localization, which results in a blurred TFD. With too many samples, on the other hand, the CS area would contain the cross-terms along with the auto-terms, resulting in a cross-term corrupted TFD. The CS area size should be [27]: card(Az (ν, τ )) ≈ Nc Nt log(Nt Nω ), with the lower limit of card(Az (ν, τ )) ≈ Nc Nt log(Nω /Nc ). However, the CS area selection for TF applications is poorly addressed in the literature, and most authors use a square area with Nτ = Nν ≈ Nt , which gives: card(Az (ν, τ )) ≈ Nt [13,26,27]. However, this size is of the order lower than the theoretical size mentioned above, and it is even smaller than the lower limit. This is because in the TF applications the goal is not to exactly reconstruct starting TFD, since this would result in the cross-term corrupted WVD. Instead, the goal is to obtain a new WVD-like TFD, but with highly suppressed cross-terms. Cross-terms, because of their oscillatory nature, are located away from the AF domain origin with the distance equal to the time-frequency distance between components [1], which leads us to a conclusion that choosing the CS area independent of AF geometry is suboptimal; when dealing with closely located components, the CS area is limited by the cross-terms location; however, when components are more spread-out, the CS area can be significantly larger. This has served as a motivation for development of a method for adaptive CS area selection, with the general goal to capture as large as possible area around the AF origin, without including any cross-terms (if included, they would reappear in the sparse TFD). The reason behind this specific approach is to lower the requirements of the reconstruction algorithm; it is easier to reconstruct a signal having more samples to begin with. The parameters Nτ and Nν in the proposed CS area selection method are picked to be equal to the time and frequency distance between the two closest components in the TF domain, respectively, or equivalently in the AF domain, to the lag and doppler distance between the origin and the first pair of cross-terms. The proposed method is based on the AF domain approach, by searching zero doppler (Az (0, τ )) and zero lag (Az (ν , 0)) slices for the first energy spike located away from the origin. The detected points on the lag and doppler axis are in fact the points where the closest pair of the cross-terms is located. The algorithm steps are as follows: 1. Since the AF domain is symmetrical, just one half of the lag and doppler AF zero slices are extracted: SL (τ ) = Az (0, Nτ /2 → Nτ ), and SD (ν ) = Az (Nν /2 → Nν , 0 ). The notation N1 → N2 is used to define a subset of samples of the respective AF zero slice that are preserved for further processing. 2. In the next step, the derivatives of the zero-doppler and the zero-lag slices are calculated. The AF zero slices are highly os-

231

cillatory, and thus their derivatives would also oscillate around the respectful axis, meaning that the cross-terms location information would be hard to extract from them. This problem is solved by filtering the AF zero slices with a moving average ∗ (ν ), refilter. The filtered slices are denoted as SL∗ (τ ) and SD spectfully. From extensive simulations, we have found out that the moving average filter window needs to be long enough to smoothen up the respectful slice, revealing its global trend-line. The conducted simulations involving various test signals have shown that Nτ /12 and Nν /12 are a good choice for the moving average filters window sizes. 3. Finally, Nτ and Nν are calculated by searching dSL∗ (τ )/dτ and ∗ (ν )/d ν for a location of local maximum after the first crossd SD ing from negative to positive values. By doing so, auto-terms located at the beginning of the slices, can settle down, ensuring that detected peak is caused by the cross-terms. To make the method even more robust, two additional conditions have been included, which lead to the final area size: 1) the next sample ∗ (ν )/d ν must be positive, and 2) it has to of dSL∗ (τ )/dτ and dSD be above a preset threshold value (in the examples to follow, the threshold value is set to 20% of the slices maximum value). The proposed method can detect cross-terms location, even in the case when cross-terms do not intersect with the AF axes, because from the AF volume distribution relationships [34]:



∞ −∞



∞ −∞

|Az ( ν , τ )|2 d τ = |Az ( ν , τ )|2 d ν =



∞

−∞



∞ −∞

|Az (0, τ )|2 e− j2π ντ dτ ,

(15a)

|Az (ν, 0 )|2 e j2π ντ dν,

(15b)

it can be concluded that cross-terms are, in fact, detectable on AF zero slices, even when they are not intersecting with it. Another way to understand this is to look at the AF as a twodimensional auto-correlation function. This approach can be confirmed by combing (4) and (6), and substituting ν = 0, in order to obtain the AF zero doppler slice, Az (0, τ ):

Az ( 0, τ ) =



∞

−∞



z t+

τ  ∗

τ

2

2

z

t−

dτ ,

(16)

which gives the well known expression for the one-dimensional auto-correlation function. Other AF slices can be interpreted as cross-correlation between the signal z(t) and its frequency shifted copies z(t )e− j2π ν t . Furthermore, it is well known that auto-correlation, along with the cross-correlation functions exhibit cross-terms between components. This fact lead us to the same conclusion as before, that the cross-terms do not need to intersect with the AF axis, in order to make their location detectable by the proposed algorithm. 3.3. Sparse signal reconstruction algorithm The problem described in (13) is a well-known unconstrained optimization problem [18], and thus it can be rewritten as:

ϑz (t, ω ) = arg min

1

ϑz (t,ω ) 2

||ϑz (t, ω ) − ψ H Az (ν, τ )||22 + τ c(ϑz ), (17)

where regularization function c (ϑz ) : R2 → R is nonsmooth and nonconvex, and τ is a regularization parameter. Since the problem  (t, ω ) is ill-posed, to calculate it we need to apply of estimating ϑ z a regularization function as some sort of prior information. Regularization functions are usually sparsity inducing, and in TF applications 1 norm is commonly used [8,13,27,29]. The 1 norm based regularization function sums all TFD samples:

c 1 ( ϑz ) =

 i

|ϑzi (t, ω )|,

(18)

232

I. Volaric et al. / Signal Processing 141 (2017) 229–239

Fig. 1. Considered three LFM component signal: (a) WVD (MzS = 3.1258), (b) RGK with kernel volume parameter α = 3 (MzS = 0.3372), (c)AF with manually selected CS-AF area with Nτ = Nν = 15 (gray rectangle), and CS-AF area obtained by the proposed method with Nτ = 73 and Nν = 17 (black rectangle). Table 1 Sparse TFD concentration comparison of the three LFM component signal reconstructed with the various 1 norm based optimization algorithms. IST[15]

TwIST[16]

GPSR[17]

SpaRSA[18]

FISTA[19]

SALSA[20]

YALL1[21]

NESTA[22]

 = 10%

MAUTO MFIX

0.2270 0.3055

0.1553 0.1878

0.2725 0.3647

0.1550 0.2181

0.3393 0.5101

0.2977 0.3815

0.0484 0.0653

0.4798 0.5908

 = 1%

MAUTO MFIX

0.2108 0.2945

0.1553 0.1725

0.1569 0.2656

0.1550 0.2016

0.3393 0.4481

0.2416 0.2980

0.0243 0.0452

0.1668 0.2602

 = 0.1%

MAUTO MFIX

0.0990 0.1556

0.1104 0.0776

0.0856 0.1386

0.0971 0.1128

0.0728 0.4481

0.1913 0.2059

0.0106 0.0227

0.1037 0.1564

The bold value indicates the better performing of the two approaches for the CS-AF area selection (manual or automatic).

Fig. 2. Reconstructed sparse TFD of the three LFM component signal with GPSR algorithm [17]: (a) Fixed CS-AF area with  = 0.1%, MzS = 0.1368; (b) Fixed CS-AF area with  = 1%, MzS = 0.2656; (c) Fixed CS-AF area with  = 10%, MzS = 0.3647; (d) Automatically selected CS-AF area with  = 0.1%, MzS = 0.0856; (e) Automatically selected CS-AF area with  = 1%, MzS = 0.1569; (f) Automatically selected CS-AF area with  = 10%, MzS = 0.2725.

I. Volaric et al. / Signal Processing 141 (2017) 229–239

233

Fig. 3. Bat echolocation signal: (a) WVD (MzS = 2.0656 ); (b) RGK-TFD, with α = 3 (MzS = 0.6051 ); (c) AF with the manually selected CS-AF area with Nτ = Nν = 19 (gray rectangle), and the CS-AF area obtained by the proposed method with Nτ = 71 and Nν = 22 (black rectangle).

Fig. 4. Newborn EEG seizure signal: (a) WVD (MzS = 1.6256 ); (b) CKD, with c = 0.01, D = 0.075, and E = 0.075 (MzS = 0.5265 ); (c) AF with the manually selected CS-AF area with Nτ = Nν = 15 (gray rectangle), and the CS-AF area obtained by the proposed method with Nτ = 23 and Nν = 15 (black rectangle).

and the optimization problem (17) becomes:

ϑz1 (t, ω ) = arg minϑz (t,ω ) subject to:

||ϑz (t, ω )||1 , ||ϑz (t, ω ) − ϒ Az (ν, τ )||22 ≤  ,

(19)

where  is a user defined energy threshold value. Eq. (19) has a unique well-known closed form:

ϑzi1 (t, ω ) = soft(ϑzi (t, ω ), τ ),

(20)

where soft(x, τ ) is a soft-threshold function defined as:

soft(x, τ ) = sign(x ) max{|x| − τ , 0}.

(21)

Over the years, a number of algorithms have been proposed for solving (19). In this paper we will use commonly used state-ofthe-art algorithms based on the basis pursuit: iterative shrinkage thresholding (IST) [15], two-step iterative shrinkage thresholding (TwIST) [16], gradient projection for sparse reconstruction (GPSR)[17], sparse reconstruction by separable approximation (SpaRSA) [18], fast IST (FISTA) [19], split augmented Lagrangian shrinkage algorithm (SALSA) [20], your augmented Lagrangian algorithm for 1 (YALL1) [21], and NESTA [22].

2 z2 (t ) = e j2π (0.00117t ) ,

(22b)

2 z3 (t ) = 1.5e j2π (0.3t+0.00078t ).

(22c)

The first component, z1 (t), has been modeled in such a way that it is almost completely masked by the cross-terms between the other two components, making it almost unidentifiable in the WVD, as shown in Fig. 1(a). Fig. 1(b) shows RGK of the considered signal with the kernel volume parameter α = 3. The current CS area selection methodology gives Nτ = Nν = 15, while the here-proposed method results in a larger CS area with Nτ = 73 and Nν = 17, as shown in Fig. 1(c). Both CS-AFs have been processed by a number of state-of-the-art 1 based optimization algorithms, for three values of the optimization stoping criteria,  , with the resulting TFD concentration being measured by the TFD concentration measure, MzS , introduced in [31]:



MzS

  1/2 1 = ϑzN1 i (t, ω ) Nt Nω

2

,

(23)

i

4. Experimental results





where ϑzN1 (t, ω ) is ϑzi1 (t, ω ) normalized over the total energy. i

Example 1. A multi-component signal with Nt = 256 samples is composed of Nc = 3 LFMs: 2 z1 (t ) = 0.5e j2π (0.15t+0.00098t ),

(22a)

Smaller MzS values correspond to more concentrated TFDs. The resulting concentration measures are given in Table 1. The proposed method resulted in smaller concentration measures for most of the tested algorithms, across all stopping criteria. The

234

I. Volaric et al. / Signal Processing 141 (2017) 229–239 Table 2 Sparse TFD concentration measure comparison of the bat echolocation signal reconstructed with the various 1 norm based optimization algorithms. IST[15]

TwIST[16]

GPSR[17]

SpaRSA[18]

FISTA[19]

SALSA[20]

YALL1[21]

NESTA[22]

 = 10%

MAUTO MFIX

0.0162 0.0229

0.0158 0.0222

0.0164 0.0226

0.0147 0.0114

0.0194 0.0310

0.0769 0.1060

0.0141 0.0168

0.3398 0.4021

 = 1%

MAUTO MFIX

0.0162 0.0229

0.0158 0.0222

0.0164 0.0226

0.0155 0.0144

0.0135 0.0197

0.0359 0.0551

0.0079 0.0091

0.2678 0.3122

 = 0.1%

MAUTO MFIX

0.0131 0.0203

0.0139 0.0203

0.0106 0.0157

0.0095 0.0110

0.0093 0.0118

0.0150 0.0158

0.0037 0.0041

0.2530 0.2801

The bold value indicates the better performing of the two approaches for the CS-AF area selection (manual or automatic). Table 3 Sparse TFD concentration measure comparison of the newborn EEG seizure signal reconstructed with the various 1 norm based optimization algorithms. IST[15]

TwIST[16]

GPSR[17]

SpaRSA[18]

FISTA[19]

SALSA[20]

YALL1[21]

NESTA[22]

 = 10%

MAUTO MFIX

0.3059 0.2992

0.2144 0.1713

0.3422 0.3492

0.2343 0.1938

0.4678 0.4852

0.3665 0.3679

0.0340 0.0305

0.4829 0.5002

 = 1%

MAUTO MFIX

0.3059 0.2767

0.1694 0.1632

0.2408 0.2079

0.1555 0.1224

0.3822 0.3914

0.3107 0.3002

0.0230 0.0277

0.2041 0.2046

 = 0.1%

MAUTO MFIX

0.1581 0.1383

0.0660 0.0776

0.1319 0.1198

0.1236 0.1121

0.1124 0.1051

0.2693 0.2461

0.0123 0.0161

0.1287 0.1452

The bold value indicates the better performing of the two approaches for the CS-AF area selection (manual or automatic).

Fig. 5. Reconstructed sparse TFD of the bat echolocation signal with NESTA algorithm [22]: (a) Fixed CS-AF area with  = 0.1%, MzS = 0.2801; (b) Fixed CS-AF area with  = 1%, MzS = 0.3122; (c) Fixed CS-AF area with  = 10%, MzS = 0.4021; (d) Automatically selected CS-AF area with  = 0.1%, MzS = 0.2530; (e) Automatically selected CS-AF area with  = 1%, MzS = 0.2678; (f) Automatically selected CS-AF area with  = 10%, MzS = 0.3398.

sparse TFDs reconstructed by the GPSR algorithm [17] are shown in Fig. 2 for  = 10%, 1%, 0.1%. As it can be seen, the sparse TFDs with the CS area selected by the proposed method are significantly better concentrated, and the middle component is clearly distinguishable, which is not the case for the sparse TFDs with manu-

ally selected CS area. Furthermore, the proposed method has resulted in 35 times more concentrated TFD when compared to the WVD, and up to 4 times when compared to the RGK-TFD (shown in Fig. 1(b)).

I. Volaric et al. / Signal Processing 141 (2017) 229–239

235

Fig. 6. Reconstructed sparse TFD of the newborn EEG seizure signal with YALL1 algorithm [21]: (a) Fixed CS-AF area with  = 0.1%, MzS = 0.0161; (b) Fixed CS-AF area with  = 1%, MzS = 0.0277; (c) Fixed CS-AF area with  = 10%, MzS = 0.0305; (d) Automatically selected CS-AF area with  = 0.1%, MzS = 0.0123; (e) Automatically selected CS-AF area with  = 1%, MzS = 0.0230; (f) Automatically selected CS-AF area with  = 10%, MzS = 0.0340.

Fig. 7. Noise impact factor of the sparse TFD reconstructed with the SALSA algorithm with the manually selected CS-AF area (dashed line) and the adaptively detected CS-AF area (solid line): (a)  = 0.1%, (b)  = 1%, (c)  = 10%.

Example 2. We have tested the proposed method on two real-life signals: the bat echolocation signal1 with Nt = 400 samples, and a newborn EEG seizure signal [32] with Nt = 256 samples , shown in Figs. 3 and 4, respectively. The WVD and the RGK-TFD (with α = 3) of the bat echolocation signal are shown in Fig. 3(a) and (b), re-

1 The authors wish to thank Curtis Condon, Ken White, and Al Feng of the Beckman Institute of the University of Illinois for the bat data and for permission to use it in this paper.

spectively, while the WVD and the CKD (with c = 0.01, D = 0.075, and E = 0.075) of the newborn EEG seizure signal are shown in Fig. 4(a) and (b), respectively (note that it has been shown in [32] that CKD with these parameters has the best resolution for the considered newborn EEG seizure signal). The AF of the bat echolocation signal is shown in Fig. 3(c), along with the CS-AF areas of the standard approach (Nτ = Nν = 19) and the proposed CSAF area selection method (Nτ = 71, Nν = 22), while the AF of the newborn EEG seizure signal is shown in Fig. 4(c) with the standard

236

I. Volaric et al. / Signal Processing 141 (2017) 229–239

Fig. 8. Reconstructed noisy sparse TFD of the three LFM component signal with SALSA algorithm and manually selected CS-AF area: (a)  = 0.1%, SNR = 5 dB, NI = 13.12%, (b)  = 1%, SNR = 5 dB, NI = 11.83%, (c)  = 10%, SNR = 5 dB, NI = 11.11%, (d)  = 0.1%, SNR = 10 dB, NI = 3.42%, (e)  = 1%, SNR = 10 dB, NI = 3.26%, (f)  = 10%, SNR = 10 dB, NI = 3.04%, (g)  = 0.1%, SNR = 20 dB, NI = 0.34%, (h)  = 1%, SNR = 20 dB, NI = 0.31%, (i)  = 10%, SNR = 20 dB, NI = 0.29%.

CS-AF area (Nτ = Nν = 15) and the proposed CS-AF area selection method (Nτ = 23 and Nν = 15). The signals TFDs performance of the 1 norm based sparse reconstruction algorithm are shown in Tables 2 and 3. The reconstructed 1 norm based sparse TFD of the bat echolocation signal using NESTA algorithm is shown in Fig. 5, while the reconstructed newborn EEG seizure signal using YALL1 algorithm is shown in Fig. 6. As it can be seen, 1 norm based sparse reconstruction of the bat echolocation signal when CS-AF area is selected by the hereproposed method resulted in a significantly more concentrated sparse TFDs, when compared to the sparse TFDs reconstructed with the fixed CS-AF area. This fact is substantiated by the visual

comparison of the respective TFDs, shown in Fig. 5, and by comparison of the concentration measures listed in Table 2. Moreover, all of the reconstructed sparse TFDs have lower concentration measure than the WVD (shown in Fig. 3(a)), and the RGK-TFD (shown in Fig. 3(b)). Visually, performance of the RGK-TFD is similar to the performance of the sparse TFDs shown in Fig. 6(b) and (f); however, imposing a stricter signal reconstruction stopping criterium has resulted in the visually more concentrated TFDs. The newborn EEG seizure signal, as it can be seen from Figs. 4 and 6 is composed of sinusoids and spikes [32]. The performance of the CS based solutions is better than the classical approach, as observed by visually comparing the respective Figures,

I. Volaric et al. / Signal Processing 141 (2017) 229–239

237

Fig. 9. Reconstructed noisy sparse TFD of the three LFM component signal with SALSA algorithm and adaptively selected CS-AF area: (a)  = 0.1%, SNR = 5 dB, NI = 21.42%, (b)  = 1%, SNR = 5 dB, NI = 19.99%, (c)  = 10%, SNR = 5 dB, NI = 18.71%, (d)  = 0.1%, SNR = 10 dB, NI = 8.98%, (e)  = 1%, SNR = 10 dB, NI = 8.73%, (f)  = 10%, SNR = 10 dB, NI = 8.63%, (g)  = 0.1%, SNR = 20 dB, NI = 4.03%, (h)  = 1%, SNR = 20 dB, NI = 4.52%, (i)  = 10%, SNR = 20 dB, NI = 4.86%.

and confirmed by the lower TFD concentration measures shown in Table 3. The reconstructed TFDs with the fixed CS-AF area have resulted in lower concentration measures for most of the tested algorithms, when compared to the TFDs reconstructed from the adaptively detected CS-AF area. Note, however, that the lower concentration measures of TFDs with the fixed CS-AF area are misleading, as they are caused by the sparse reconstruction algorithms failure to correctly reconstruct the spiky components, which is not the case when the CS-AF area is selected by the here-proposed method, as visually confirmed by the plots shown in Fig. 6.

Example 3. Next, we have embedded the test signal from Example 1 in additive white Gaussian noise with SNR ranging from 0 to 20 dB, and averaged the experiment results over 100 simulations. We have then measured the noise impact factor, NI, calculated as:

 2  ϑzi1 (t, ω ) − ϑzi1REF (t, ω ) · 100%,  2   1  ϑziREF (t, ω )

 NI =

i

i

(24)

238

I. Volaric et al. / Signal Processing 141 (2017) 229–239

Fig. 10. TFD reconstruction with the elliptical CS-AF mask: (a) CS-AF area with N = 1935 samples, (b) reconstructed sparse TFD (SpaRSA) with  = 0.1%, (MzS = 0.1220 ), (c) reconstructed sparse TFD (SpaRSA) with  = 0.1%, from the CS-AF containing 193 (10%) of the maximal AF samples inside the elliptical mask (MzS = 0.0715 ).



where ϑzi1 (t, ω ) is the reconstructed noiseless sparse TFD. The REF

noise impact factor for the SALSA algorithm with  = 0.1%, 1%, 10% as a function of the noise level is shown in Fig. 7. As it can be seen, the NI factor is very similar for all the stopping criteria; moreover, the behavior of the NI factor is similar across all of the tested sparse reconstruction algorithms. The simulations have shown that the proposed method is more susceptible to noise. The reason for this is twofold: firstly, the proposed method results in a larger CS area, and secondly, the noise influences the area selection method itself. A larger CS-AF area means that the reconstruction algorithm input data has overall more noise corrupted samples. The influence of the CS-AF area size on the NI factor is thoroughly investigated in [35]. The proposed method in presence of noise can also result in a smaller CS-AF area, making the resulting TFD less concentrated. In such a case, the NI factor becomes larger since the difference be tween the noiseless (and more concentrated) TFD, ϑzi1 (t, ω ), and 

REF

the noisy (and less concentrated) TFD, ϑzi1 (t, ω ) increases. Figs. 8 and 9 show the noise corrupted TFDs reconstructed with the SALSA algorithm with the manually detected CS-AF area and the adaptively detected CS-AF area, respectively. In both cases, when the noise SNR level is above 5 dB, the signal features are clearly identifiable in the resulting sparse TFD. The added AWGN does not greatly affect the concentration of the resulting sparse TFDs, and the conclusions made in Example 1 are also valid when the considered signal is noise corrupted. Example 4. Finally, we explore non-rectangular possibilities of the CS-AF mask construction using the parameters Nτ and Nν obtained by the proposed algorithm. The first design which naturally arises is an elliptical mask, designed to be tangential to the cross-terms, leading to the following expression for the sensing matrix:

φ (ν , τ ) =

1, 0,

2ν 2 Nν2

+ 2Nτ2 ≤ 1, τ otherwise.

inside the elliptical mask defined by Eq. (25). The reconstructed sparse TFD from a such defined CS-AF area, shown in Fig. 10(c), is significantly better concentrated than the previously considered sparse TFDs (shown in Figs. 2 and 10(b)), achieving the lowest concentration measure value MzS = 0.0715, despite the fact that the respectful mask is the smallest, containing only 194 samples. Indeed, by enlarging the continuous mask, only a small subset of the newly acquired samples provide useful information on the autoterms, while the most of the so acquired samples are close to zero, hence providing little or no useful information. This can be observed in Figs. 1(c), 3(c), and 10(a) where the CS-AF masks contain mostly the blank space, while only a small number of captured samples originate from the auto-terms. These information-less CSAF samples impose additional constraints in signal reconstruction, resulting in a slower and stricter optimization. This is the reason behind the effectiveness of the considered mask; it captures only the samples with the most auto-terms information, while excluding the other less informative samples. 5. Conclusion We have presented a method for an automatic data-driven selection of the CS area for enhancement of the signal TFD. The CS area is generally selected experimentally, with the area size being approximately Nt samples in the ambiguity domain. The method proposed in this paper automates this procedure, resulting in CS areas that contain more useful samples. Furthermore, the proposed method also allows a less strict optimization criterion of the sparse signal reconstruction algorithm, ultimately resulting in highly concentrated TFDs. References

2

(25)

The sensing mask area defined in this way is π /2 times larger than the rectangular mask, allowing more samples to be included from which the TFD is reconstructed. We have tested this design on the three component LFM signal from Example 1. The elliptical CS-AF mask is shown in Fig. 10(a) containing 1935 samples, which is slightly more than the rectangular area containing 73 × 17 = 1241 samples. The reconstructed sparse TFD is shown in Fig. 10(b) achieving MzS = 0.1220, thus outperforming the rectangular CS-AF area by 15% in terms of the resulting TFD concentration. The second tested CS-AF area design is not a continuous area, like the previously discussed CS-AF masks. Instead, the CS-AF area is designed to include only 10% of the maximum values of the AF

[1] B. Boashash, Time frequency signal analysis and processing: a comprehensive reference 2nd edition, Elsevier, 2016. [2] D. Donoho, Compressed sensing, IEEE Trans. Inf. Theory 52 (4) (2006) 1289– 1306, doi:10.1109/TIT.2006.871582. [3] R. Baraniuk, Compressive sensing [lecture notes], IEEE Signal Process. Mag. 4 (24) (2007) 118–121. [4] E. Candes, M. Wakin, An introduction to compressive sampling, Signal Process. Mag., IEEE 25 (2) (2008) 21–30, doi:10.1109/MSP.2007.914731. [5] J. Romberg, Imaging via compressive sampling [introduction to compressive sampling and recovery via convex programming], IEEE Signal Process. Mag. 25 (2) (2008) 14–20. [6] M. Lustig, D. Donoho, J. Santos, J. Pauly, Compressed sensing MRI, Signal Process. Mag., IEEE 25 (2) (2008) 72–82, doi:10.1109/MSP.2007.914728. [7] A. Gramfort, D. Strohmeier, J. Haueisen, M.S. Hämäläinen, M. Kowalski, Time-frequency mixed-norm estimates: sparse m/eeg imaging with non-stationary source activations, NeuroImage 70 (2013) 410–422. [8] A. Gholami, Sparse time-frequency decomposition and some applications, IEEE Trans. Geosci. Remote Sensing 51 (6) (2013) 3598–3604.

I. Volaric et al. / Signal Processing 141 (2017) 229–239 [9] C. Gangodagamage, E. Foufoula-Georgiou, S.P. Brumby, R. Chartrand, A. Koltunov, D. Liu, M. Cai, S.L. Ustin, Wavelet-compressed representation of landscapes for hydrologic and geomorphologic applications, IEEE Geosci. Remote Sensing Lett. 13 (4) (2016) 480–484. [10] J.E. Barcelo-Llado, A. Morell, G. Seco-Granados, Amplify-and-forward compressed sensing as an energy-efficient solution in wireless sensor networks, IEEE Sensors J. 14 (5) (2014) 1710–1719. [11] S. Pudlewski, A. Prasanna, T. Melodia, Compressed-sensing-enabled video streaming for wireless multimedia sensor networks, IEEE Trans. Mobile Comput. 11 (6) (2012) 1060–1072. [12] S. Stankovic, I. Orovic, M. Amin, L-statistics based modification of reconstruction algorithms for compressive sensing in the presence of impulse noise, Signal Process. 93 (11) (2013) 2927–2931. http://dx.doi.org/10.1016/j.sigpro.2013. 04.022. [13] L. Stankovic, I. Orovic, S. Stankovic, M. Amin, Compressive sensing based separation of nonstationary and stationary signals overlapping in time-frequency, IEEE Trans. Signal Process. 61 (18) (2013) 4562–4572, doi:10.1109/TSP.2013. 2271752. [14] L. Stankovic, L. S. Stankovic, L. I. Orovic, M.G. Amin, Robust time-frequency analysis based on the l-estimation and compressive sensing, IEEE Signal Process. Lett. 20 (5) (2013c) 499–502. [15] I. Daubechies, M. Defrise, C. De Mol, An iterative thresholding algorithm for linear inverse problems with a sparsity constraint, Commun. Pure Appl. Math. 57 (2004) 1416–1457. [16] J.M. Bioucas-Dias, M.A. Figueiredo, A new TwIST: two-step iterative shrinkage/thresholding algorithms for image restoration, IEEE Trans. Image Process. 16 (12) (2007) 2992–3004. [17] M.A. Figueiredo, R.D. Nowak, S.J. Wright, Gradient projection for sparse reconstruction: application to compressed sensing and other inverse problems, IEEE J. Sel. Topics Signal Process. 1 (4) (2007) 586–597. [18] S.J. Wright, R.D. Nowak, M.A. Figueiredo, Sparse reconstruction by separable approximation, IEEE Trans. Signal Process. 57 (7) (2009) 2479–2493. [19] A. Beck, M. Teboulle, A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM J. Imaging Sci. 2 (1) (2009) 183–202. [20] M. Afonso, J. Bioucas-Dias, M. Figueiredo, Fast image recovery using variable splitting and constrained optimization, IEEE Trans. Image Process. 19 (9) (2010) 2345–2356, doi:10.1109/TIP.2010.2047910. [21] Y. Zhang, Users guide for YALL1: Your ALgorithms for l1 optimization, Technique report, 2009. [22] S. Becker, J. Bobin, E.J. Candès, NESTA: a fast and accurate first-order method for sparse recovery, SIAM J. Imaging Sci. 4 (1) (2011) 1–39.

239

[23] S. Mallat, Z. Zhang, Matching pursuits with time-frequency dictionaries, IEEE Trans. Signal Process. 41 (12) (1993) 3397–3415, doi:10.1109/78.258082. [24] T. Cai, L. Wang, Orthogonal matching pursuit for sparse signal recovery with noise, IEEE Trans. Inf. Theory 57 (7) (2011) 4680–4688, doi:10.1109/TIT.2011. 2146090. [25] I. Orovic, S. Stankovic, T. Thayaparan, Time frequency-based instantaneous frequency estimation of sparse signals from incomplete set of samples, Signal Process., IET 8 (3) (2014) 239–245, doi:10.1049/iet-spr.2013.0354. [26] P. Borgnat, P. Flandrin, Time-frequency localization from sparsity constraints, in: 2008 IEEE International Conference on Acoustics, Speech and Signal Processing, 2008, pp. 3785–3788. [27] P. Flandrin, P. Borgnat, Time-frequency energy distributions meet compressed sensing, IEEE Trans. Signal Process. 58 (6) (2010) 2974–2982, doi:10.1109/TSP. 2010.2044839. [28] P. Flandrin, N. Pustelnik, P. Borgnat, On wigner-based sparse time-frequency distributions, in: Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), 2015 IEEE 6th International Workshop on, 2015, pp. 65–68. [29] I. Volaric, V. Sucic, Z. Car, A compressive sensing based method for crossterms suppression in the time-frequency plane, in: Bioinformatics and Bioengineering (BIBE), 2015 IEEE 15th International Conference on, 2015, pp. 1–4. doi:10.1109/BIBE.2015.7367703. [30] Y. Wang, J. Xiang, Q. Mo, S. He, Compressed sparse time-frequency feature representation via compressive sensing and its applications in fault diagnosis, Measurement 68 (2015) 70–81. [31] L. Stankovic, A measure of some time-frequency distributions concentration, Signal Process. 81 (3) (2001) 621–631. Special section on Digital Signal Processing for Multimedia. http://dx.doi.org/10.1016/S0165-1684(0 0)0 0236-X. [32] B. Boashash, S. Ouelha, Automatic signal abnormality detection using timefrequency features and machine learning: a newborn EEG seizure case study, Knowl.-Based Syst. 106 (2016) 38–50. http://doi.org/10.1016/j.knosys.2016.05. 027. [33] D.L. Jones, R.G. Baraniuk, An adaptive optimal-kernel time-frequency representation, IEEE Trans. Signal Process. 43 (10) (1995) 2361–2371. [34] N. Levanon, E. Mozeson, Radar Signals, John Wiley & Sons, 2004. [35] I. Volaric, V. Sucic, On the noise impact in the l1 based reconstruction of the sparse time-frequency distributions, in: 2016 International Conference on Broadband Communications for Next Generation Networks and Multimedia Applications (CoBCom), 2016, pp. 1–6.