S-transform based on compact support kernel

Accepted Manuscript S-transform based on compact support kernel

Z. Zidelmal, H. Hamil, A. Moukadem, A. Amirou, D. Ould-Abdeslam

PII: DOI: Reference:

S1051-2004(16)30203-2 http://dx.doi.org/10.1016/j.dsp.2016.11.008 YDSPR 2049

To appear in:

Digital Signal Processing

Please cite this article in press as: Z. Zidelmal et al., S-transform based on compact support kernel, Digit. Signal Process. (2016), http://dx.doi.org/10.1016/j.dsp.2016.11.008

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Highlights • In this study, a Compact Support Kernel (CSK) is used instead of the Gaussian window in the S-transform. • The width of the CSK is controlled by parameters making it more flexible. • These parameters are automatically selected to maximize the energy concentration.

S-transform based on Compact Support Kernel Z.Zidelmal,a, , H.Hamil,a, , A.Moukadem,c, , A.Amirou,b, , D.Ould-Abdeslam,c, b

a D´epartement d’Electronique, Universit´e Mouloud Mammeri, Tizi-Ouzou, Algerie D´epartement de Math´ematiques, Universit´e Mouloud Mammeri, Tizi-Ouzou, Algerie c Laboratoire MIPS,Universit´e de Haute Alsace,France

Abstract In this paper, the use of a Compact Support Kernel (CSK) instead of the Gaussian window in the S-transform is proposed. The CSK is derived from the Gaussian but overcomes its practical drawbacks while preserving a large number of its useful properties. The width of the CSK is controlled by some parameters making it more flexible. These parameters are selected to optimize the energy concentration in the time-frequency domain. Compared to other versions of S-transform, other time-frequency representations and continuous wavelet transform, the achieved results obtained using synthetic and real data show a significant improvement in the time and frequency resolution, energy concentration and instantaneous frequency estimation. Key words: S-Transform, Gausssian window, Compact Support Kernel (CSK), Energy concentration. 1. Introduction Extraction of pertinent information from noisy, multicomponent and nonstationary signals with various and complex backgrounds remains a challenging problem in signal processing. These facts have motivated the use of frequency domain techniques, such as Fourier transform (FT) for the analysis. While the FT does not represent any limitation on the analysis of time-invariant signals, it presents an important handicap when time-variant Email addresses: [email protected] (Z.Zidelmal ), [email protected] (H.Hamil ), [email protected] (A.Moukadem ), [email protected] (A.Amirou ), [email protected] (D.Ould-Abdeslam )

Preprint submitted to Digital Signal Processing

November 29, 2016

signals are studied. To solve these problems, many time-frequency analysis tools have been introduced to process non-stationary signals in the joint time-frequency domain. The short-time Fourier transform (STFT) [1] is one of the most basic tools used for nonstationary signals. This approach introduces a sliding window in the Fourier integral to achieve a better estimation of the time distribution of each frequency component. However, the STFT suffers from the undesirable tradeoff between the time concentration and the frequency concentration due to the limitations of a fixed window width chosen a priori. A Wigner-Ville Distribution is a quadratic transform [2, 3, 4]. It is known to be optimal for monocomponent linear frequency modulation (LFM) signals for any modulation rate [4] since it achieves the best energy concentration around the signal instantaneous frequency law [3]. However, it suffers from cross terms when analyzing multicomponent signals [2, 3]. Over the last few decades, an other tool, the continuous wavelet transform (CWT) has been proposed [5, 6] in order to overcome the fixed window in the STFT. The CWT has been applied to a wide variety of signal processing problems [7, 8]. It can be seen as an extension of the spectrogram in a wide sense, except for representing signals in time-scale space instead of time-frequency space with a frequency variant amplitude response and octave scaling of the frequencies. The CWT is based on scalable and translatable wavelets. The phase of the wavelet transform is relative to the analyzing wavelet center. Thus, as the wavelet translates, the phase reference translates. A more recently developed time-frequency representation method is the S-transform (ST) introduced by R.G.Stockwell [9, 10]. The S-transform is a hybrid of the STFT and the CWT [11] when combining their good features. It uses a variable analyzing window width and preserves the phase information. Furthermore, it maintains a direct relationship with the Fourier spectrum. The S-transform uniquely combines three characteristics: i) progressive resolution. ii) absolutely referenced phase information. iii) frequency invariant amplitude response [11, 12]. In addition, it allows us access to any frequency component or any single frequency without requiring to digital filters. This combination of desirable features has shown promise in a wide variety of applications such as power quality [13, 14], Biomedical engineering [15, 16] and geophysics [17]. The original S-transform uses a Gaussian window chosen for several reasons enumerated in [18]. In addition, the Gaussian functions are the only solution that minimizes the duration bandwidth product posed by Heisenberg uncertainty principle [2]. In the S-transform formulation, the Gaussian 2

standard deviation is the inverse of the frequency whatever the analyzed signal [9, 10, 18]. Nevertheless, some problems considered as limitations of the original ST may be cited: i) the frequency dependence of horizontal and vertical dilatation of the analyzing window. This window has no parameters to allow its width in time or frequency to be adjusted. This leads to degradation of the time resolution at low frequencies and a poor frequency resolution at high frequencies. ii) Because of the infinite support of the Gaussian, its exact implementation is not possible. It must be truncated to a finite window when implemented, leading to information loss. Improving this window attracted much research interest. several modifications were introduced in the original S-transform by proposing new windows. Pinnegar et al. introduced windows having frequency dependence in their shape in addition to their width and height [20]. In [21], the authors introduced a bi-Gaussian window which seems better at resolving the sharp onset of events in a time series. The window proposed in [22] is based on the T-student distribution in order to obtain a more uniform time and frequency resolutions. In [23], a frequency dependent Kaiser window is used for improving the energy concentration. A modified S-transform is also proposed in [24, 25, 26] allowing a better control of the time and frequency resolution by a progressive control of the Gaussian window width. A parameter which maximizes the energy concentration is also selected in [24]. Recently, Moukadem et al. [27] introduced a new version of S-transform where energy concentration enhancement is provided. It also generalizes the standard S-transform and the modified ones presented in [24, 25, 26]. Motivated by these interesting properties, the authors propose in this contribution the use of a compact support kernel (CSK) instead of the Gaussian window in the S-transform. The CSK is derived from the Gaussian kernel and preserves a large number of its useful properties. However, as the Gaussian kernel has an infinite support to be truncated when implemented, the CSK vanishes itself outside a given compact support. The CSK should provide the same time-frequency tiling as the frequency dependent Gaussian window and should maintain the desirable properties of the standard S-transform. The CSK width is controlled by some parameters adjusted by an optimization algorithm where the objective function is maximizing the energy concentration in the time-frequency domain. The paper is organized as follows. In the next section, the construction of the CSK is detailed and then, the expression of the modified S-transform, namely the CSK-ST is given. In this section, the optimization of the scal3

ing parameter is also presented. In section 3, the achieved results obtained with the CSK-ST are presented and compared to those obtained with others approaches. An application on pregnant women ECG is also given to reveal the temporal resolution quality obtained with this proposal. Finally, Section 4 concludes the paper. 2. Method 2.1. The standard S-Transform The S-transform originates from two advanced signal processing tools, the short time Fourier transform (STFT) and the continuous wavelet transform (CWT). Derived from the STFT, the standard S-transform of a time varying signal x(t) is given by:[9, 10]  +∞ S(τ, f ) = x(t)w(t − τ )e−i2πf t dt. (1) −∞

where w(t) is a time window centered in t = 0 and used to extract a segment of x(t). The S-transform can be found by defining a particular window function w(t), a normalized gaussian −t2 1 w(t) = √ e 2σ2 . (2) σ 2π and allowing w(t) a translation τ and a dilatation (a variable width σ). A constraint is added to restrict the window width σ to be a function of the frequency 1 σ(f ) = . (3) |f |

The window width σ varying inversely with frequency makes the ST performing a multi-resolution analysis on the signal. The ST (1) is rephrased as  +∞ −(t−τ )2 f 2 |f | S(τ, f ) = √ x(t)e 2 e−i2πf t dt. (4) 2π −∞ where τ and f denote respectively the time of the spectral localization and the Fourier frequency. There are two vital terms in the S-transform definition [18]: i) the phase function e−i2πf t localizing the phase spectrum as well as the amplitude spectrum. This is referred to as absolutely referenced phase information meaning that the phase information given by the S-transform 4

is always referenced to time t = 0. ii) the normalization parameter √|f2π| normalizes the localizing window to have unit area. Therefore, the amplitude of the S-transform has the same meaning as the amplitude of the Fourier transform. In equation (4), the S-transform can be seen as a continuous wavelet transform  +∞ 1 t−τ W (τ, a) = )dt. (5) x(t) √ ψ ∗ ( a a −∞ with a specific mother wavelet ψ(.) multiplied by a phase factor correction [9] and a constraint that the wavelet dilatation factor a = σ(f ). Hence, the mother wavelet is |f | −t2 f 2 ψ(t, f ) = √ e 2 e−i2πf t . (6) 2π and S(τ, f ) = W (τ, f )e−i2πf τ . (7) Note that the wavelet (6) does not satisfy the admissibility condition of having a zero mean, so (7) is not strictly a CWT. The ST returns correct amplitude response for all frequencies and produces a time frequency representation instead of the time scale representation developed by the CWT [18]. Even though advantages of the S-transform becoming a valuable tool for signal analysis in many applications, it presents some drawbacks such as: • the fact that the window width is always defined as a reciprocal of the frequency when some signals would benefit from different window widths. • the window has no parameters to make it more flexible and more adaptive to the analyzed signal. • the only considered window is a Gaussian, chosen for several reasons [18]. However, because of its infinite support, exact implementation of this kernel is not possible so, approximating the infinite support by a finite support becomes inevitable. 2.2. S-transform with Compact Support Kernel 2.2.1. Compact Support Kernel construction To remedy the standard ST limitations already mentioned above, an improvement can be obtained by replacing the Gaussian window with a frequency dependent compact support kernel which must satisfy the following conditions: 5

• having a finite support while retaining the advantage of the Gaussian. • allowing the same time-frequency tiling as that of the standard ST. • getting better progressive control of the window width in order to maximize the energy concentration. The CSK construction is inspired from some proposals [28, 29, 30] where compact support kernels were proposed. The compact support kernel proposed in [28] was applied to scale-space image processing and extracting handwritten data. In [29, 30], the authors introduced a particular quadratic TFD, a kernel-based transform [31]. This quadratic TFD may be considered as a smoothed Wigner-Ville distributions (WVDs), where the ”smoothing” is described in the time-frequency domain by convolution with a two-dimensional compact support kernel. In this paper, the CSK is derived from the Gaussian window by deforming the real axe into a finite segment. There are several transformations allowing this. Let us consider a C 1 diffeomorphism, h−1

h : R+ −→ [0, +1[, a map between manifolds which is differentiable and has a differentiable inverse. It was also shown how similar are the derivatives of the obtained kernel and those of the Gaussian [28]. The map h is defined as  t → h(t) = −ln(1 − t2 ) with 0 ≤ t < 1. (8) After applying this variable change to a normalized Gaussian kernel defined as −t2 1 gσ (t) = √ e 2σ2 , (9) σ 2π the following kernel is obtained  1 γln(1−t2 ) e = C1γ (1 − t2 )γ , if t2 < 1 , Cγ φγ (t) = (10) 0, elsewhere . where γ = 2σ1 2 . Although γ could be a real number, it is considered to be a positive integer. Hence, the proposed kernel has a polynomial form. Cγ is a normalization parameter which must satisfy  +1  +1 1 2 γ (1 − t ) dt = 1 ⇒ Cγ = (1 − t2 )γ dt. (11) C γ −1 −1 Unfortunately, there is no analytical expression for Cγ , it must be computed by a numerical method. If N is the number of points in the discretization 6

 2 2 γ process, Cγ = Δt N i=1 (1 − ti ) where ti = −1 + iΔt and Δt = N . For a variable kernel support λ, the one dimensional family of CSK is defined as:  1 (λ2 − t2 )γ , if t2 < λ2 , Cγ Dγ φλ (t) = (12) 0, elsewhere . where the normalization parameter Cγ1Dγ is introduced to insure the invertibility of the CSK-ST. Let’s put x = t/λ and dt = λdx,  +1 2γ+1  +λ 1 λ 2 2 γ (λ − t ) dt = (1 − x2 )γ dx = 1. (13) −λ Cγ Dγ −1 Cγ Dγ From equation (11) and (13), Dγ = λ2γ+1 . A shifted and scaled CSK can then be expressed as:  1 (λ2 − (t − τ )2 )γ , if (t − τ )2 < λ2 , Cγ λ(2γ+1) φλ (t − τ ) = (14) 0, elsewhere . The kernel width is controlled through λ and its peak is adjusted through the parameter γ allowing a tradeoff between a good autoterm resolution and a sufficient cross-term suppression. This parameter is empirically set to 2. As in the standard S-transform where the window width is controlled by σ(f ) = |f1 | , the scale parameter λ is also assumed to be a function of frequency and defined as . λ(f ) =

m . p + fr

(15)

The introduced parameters m, p and r in λ aim to make the CSK, more flexible and more adaptive to the analyzed signal. Hence, the resolution in time and in frequency will be tuned depending on these parameters. Figure.1 illustrates the variation rule of the CSK width against frequency, compared to that of the frequency dependent Gaussian. In Fig.2, one can compare the Gaussian window and the CSK. At low frequencies, the CSK is more compact than the Gaussian window. It should provide better time resolution. At high frequencies, the Gaussian is very narrow, inevitably with consequent loss of resolution in the frequency direction. Let’s now examine and compare the frequency domain behavior of the CSK and the Gaussian window at low and at high frequencies. As one important property of the Gaussian window is that its Fourier transform is also 7

1 0.9

Proposed CSK−ST Stadard ST

Window width variation

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0

200

400 600 Frequency(Hz)

800

1000

Figure 1: Variation law of the normalized window width over frequency: CSK with m=1.02, p=0.5, r=0.8 (red) and the Gaussian used in standard ST (black).

0.015

Amplitude(V)

Amplitude(V)

0.14

0.05

0.005

0 −0.1

−0.05

0

0.05

0.1

Time(s)

0 0.05

0 Time(s)

0.05

Figure 2: Comparison of the windows in the time domain, for normalized frequency f = 0.04Hz (left) and f = 0.35HZ (right). Frequency dependent Gaussian (dashed) and the CSK (solid) with γ = 2 (black) and γ = 5 (red).

8

Normalized PSD

theoretically Gaussian, ie positive and unimodal. This property is interesting in many applications but it is lost since the Gaussian is truncated when implemented. Therefore, neither the CSK nor the Gaussian are unimodal. The graphs in Figs.3 and 4 represent power spectra of the two windows respectively on a linear scale and on a logarithmic scale. These illustrations 1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

Normalized PSD

0 −0.2 1

−0.1

0

0.1

0.2 γ=2 γ=5

0.8

0 −0.2 1

0.6

0.4

0.4

0.2

0.2

−0.1

0

0.1

0.2

Normalized frequency

0

0.1

0.2 γ=2 γ=5

0.8

0.6

0 −0.2

−0.1

0 −0.2

−0.1

0

0.1

0.2

Normalized frequency

Figure 3: Power spectrum on a linear scale of the frequency dependant Gaussian window (top) and that of the frequency dependant CSK (bottom) for f = 0.04 (left) and f = 0.36 (right).

depict an interesting result of such analysis. As shown in Figs.3 and 4, the CSK provides much narrower mainlobe than that of the Gaussian, with significantly reduced side-lobes. It is clear that, not only the CSK provides the same time-frequency tiling as the frequency dependent Gaussian window, but it should improve the resolution in the two axis notably for γ = 2. Considering equations (14, 15) respectively of the CSK and its scale parameter λ, the CSK-ST can be expressed as

S

m,p,r

1 p + f r 2γ+1 (τ, f ) = ( ) Cγ m



+∞

x(t)[ −∞

m2 −(t−τ )2 ]γ e−i2πf t dt. (16) (p + f r )2

2.2.2. Optimization of the scaling parameter Equation (15) shows that the kernel support depends on some parameters (m, p and r). Obviously, the crucial question is how to choose these parameters. Select their values empirically may not be adequate for some types of 9

0

−100

−100

PSD

dB

0

−200

−200

−300

−300

−400 −0.5 0

−400 0.5 −0.5 0

0 γ=2 γ=5

−100

0

0.5 γ=2

−100

γ=5

PSD

dB

−200 0

−200 −100 −300 −0.5

−200 0.5 −0.5

0

0

Normalized frequency

Normalized frequency

0.5

Figure 4: Power spectrum on a logarithmic scale of the Gaussian window (top) and that of the CSK (bottom) for a normalized frequency f = 0.04 (left) and f = 0.36 (right).

signals. To provide automatic selection of these parameters with respect to the analyzed signal, the methodology proposed in [24, 27] is adopted. The objective function is maximizing the energy concentration: CM (m, p, r) =  +∞  +∞ −∞

−∞

1 |S m,p,r (τ, f )|dτ df

.

(17)

The module of the S-transform is normalized as: S m,p,r (τ, f ) . S m,p,r (τ, f ) =  +∞  +∞ m,p,r (τ, f )|2 dτ df |S −∞ −∞

(18)

To maximize the energy concentration (17), an optimization problem is formulated under some constraints related to the width λ(f ) of the kernel. This later should not be very narrow to alter the frequency resolution and not very large to destroy the time resolution KTs ≤

m ≤ LTs . p + fr

(19)

where Ts is the sampling period, f ∈ [fmin , fmax ]. fmin = 1, is assumed to be the first voice while fmax depends on the analyzed signal. As in [27], K and L are respectively fixed to 10 and 1000. Developing the double inequality 10

(19) leads to the constraints: 

r KTs fmax + pKTs − m ≤ 0 , m/(p + 1) − LTs ≤ 0.

(20)

The last constraint should limit the intervals that contain respectively the parameters m, p and r where m, p ∈ [0, 2]. The lower bound is set to 0 to ensure a positive window width. The upper limit is set to 2 according to tests we performed using a set of synthetic test signals that we will define in the following paragraphs. However, r is limited to the range [0, 1]. Any negative value of r would make the window wider as frequency increases. Similarly, values greater than 1 provide a window which may be too narrow in the time domain. Hence, the optimization problem can be set as: ⎧ max N fmax 1 m,p,r , ⎪ ⎪ (τ,f )| m,p,r ⎪ 1 fmin |S ⎪ ⎪ r ⎨ Sc : KTs fmax + pKTs − m ≤ 0 , m/(p + 1) − LTs ≤ 0 , ⎪ ⎪ ⎪ 0 ≤ m, p ≤ 2 , ⎪ ⎪ ⎩ 0≤r≤1 .

(21)

In fact, since the parameter γ in (16) controls the shape of the window in the time domain and in the frequency domain, it could be signal dependent. In this case, it must be introduced in the optimization problem. This later would then become heavier. The constrained optimization problem (21) is resolved using an active-set based algorithm [33] as used in [27]. This strategy aims to transform the problem into a basis easier one to be used in an iterative process. The implementation is performed under Matlab environment. 3. Performance analysis In this section, the performance of the proposed scheme is examined using a set of synthetic and real test signals. The goal is to examine how the CSK-ST performs compared to the standard S-transform, other modified Stransforms proposed in the literature and to other classical time-frequency representations. All used synthetic signals are with 1 second length and sampled at 1kHz. 3.1. Illustration on simulated signals As a first illustrative test, the CSK-ST is applied to a Linear Frequency Modulated (LFM) signal with transients s1 (t) including low (f1 = 100Hz), 11

high (f2 = 350Hz) frequency and two chirp components of frequency ranges [250Hz–200Hz] and [200Hz–250Hz]. The achieved result is presented in Fig.5 (left). To check the frequency resolution at low and at high frequencies, the signal s2 (t) that consists of two pairs of closely spaced parallel chirps of normalized frequency ranges [0.05–0.15], [0.1–0.2] and [0.3–0.35], [0.35–0.4] respectively is considered. Figure 5 (right) shows clean plots with notable frequency resolution.

Figure 5: Time-frequency representation of s1 (t) and s2 (t) using CSK-ST.

3.2. Comparison with other versions of S-transform In the paragraphs above, some existing versions of modified S-transform have been cited, notably that presented in [27] generalizing the standard ST and the versions introduced in [24, 25, 26]. In order to evaluate the performance of the CSK-ST compared with other versions of S-transform, the proposed approach is applied to the same synthetic signals used in [27] where a rigorous comparison with [24, 25, 26] has been made. Firstly, we consider a multi-component transient signal s3 (t) including lower (f1 = 100Hz), medium (f2 = 200Hz), and high (f3 = 400Hz) frequency components and s4 (t), a sinusoidal FM signal with a crossing linear chirp, defined as follow: s4 (t) = cos(aπt − bπt2 ) + cos[4πsin(πf t) + 80πt].

(22)

where a, b and f control the chirp rate and the sinusoidal modulated component. The results illustrated in Fig.6 (left) show that the CSK-ST with optimized parameters (m=1.0709, p=0.9985, r=0.3944), leads to a better frequency resolution than the standard S-transform and a better time resolution than Moukadem’s proposed approach. For s4 (t), the results are shown 12

in Fig.6 (right). The standard S-transform provide poor energy concentration and worse resolution, especially at low frequencies. The CSK-ST with optimized parameters (m=1.0797, p=0.9969, r=0.6239) shows comparable results to those obtained with Moukadem’s ST for the linear and the non linear components; except that Moukadem’s ST shows a slight loss in frequency resolution at very low frequencies (close to zero) for the sinusoidal FM component. 2 4

x (t)

x3(t)

1 0

−1

−2

(a)

1

500 400 300 200 100

(b)

0.5

Time(s)

100 0.5

0

100

1

(c)

Time(s)

0.5

1

0.5

1

80 60 40 20 0

1

Time(s)

0.5 Time(s)

20

100 Frequency(Hz)

(c)

(b)

40

Frequency(Hz)

200

1

60

100

300

0.5

(a)

80

1

400

500

0

100

500 Frequency(Hz)

0.5

Frequency(Hz)

Frequency(Hz)

0

Frequency(Hz)

0

400 300 200

80 60 40 20

100 0.5

0

1

Time(s)

0

Time(s)

Figure 6: Comparison between different versions of S-transform using s3 (t) (left) and s4 (t) (right) with a = 100, b = 20 and f = 5 (a) Standard ST, (b) Moukadem’s ST, (c) CSK-ST.

This comparison is also made using other kind of signals. s5 (t) is composed of four Gaussian modulated kernels. This transient signal in both time and frequency is given by: 2

2

s5 (t) = e−35π(t−t1 ) cos(πf1 t) + e−35π(t−t2 ) cos(πf1 t) 2 2 +e−55π(t−t3 ) cos(πf2 t) + e−45π(t−t3 ) cos(πf3 t). 13

(23)

where t1 , t2 , t3 , f1 , f2 and f3 control the time and frequency positions of the components. Figure 7 (left) shows that the standard S-transform suffers from poor frequency resolution at high frequencies and poor temporal resolution at low frequencies. Although Moukadem’s ST presents a good time-frequency resolution throughout over the whole image, some loss of temporal resolution compared to the CSK-ST (with m=1.0602, p=0.9981, r=0.4918) can be seen. A slightly more complicated signal s6 (t) is defined as: s6 (t) = cos[20π ln(at + 1)] + cos(bπt + cπt2 ).

(24)

where the parameters a, b and c control the two crossing components. To perform a comparison with [27] these parameters are fixed as a = 30, b = 40 and c = 150. This example is considered to be complicated: Firstly, the components are crossing. Secondly, as the frequency of the hyperbolic component decreases, the frequency of the linear chirp increases. Often, it is difficult to provide good concentration for both. The different time-frequency representations of s6 (t) are depicted in Fig.7 (right). For the hyperbolic component, the CSK-ST (with m=1.0818, p=0.9980, r=0.5893) leads to a better time resolution at high frequencies than Moukadem’s ST. At low frequencies, the CSK-ST shows a comparable result to that obtained with the standard ST. This later is very favorable for this component following the strategy of the Gaussian window: higher temporal resolution at higher frequencies and higher frequency resolution at lower frequencies. For the LMF component, the standard ST shows a poor resolution at high frequencies. However, the CSK-ST leads to comparable frequency resolution than Moukadem’s ST. 3.2.1. Effect of additive noise To check the effect of additive noise, the signals s5 (t) and s6 (t) embedded in medium (SN R = 5dB) and high (SN R = 0dB) levels of an Additive White Gaussian Noise (AWGN) have been used. The time-frequency plots are shown in Figs.8 and 9. Here, one can see that the CSK-ST is more robust to noise than the other considered versions of S-transform overall the time-frequency plane, for medium and high noise levels. Instead of relying solely on visual inspection of the presented time-frequency plots, a quantitative comparison is also made using energy concentration. Table 1 summarizes the energy concentration results obtained in the same measurement conditions with the standard ST, Moukadem’s ST and the CSK-ST when dealing with some already mentioned signals and a monocomponent LFM signal. The results show that for the analyzed signals, the CSK-ST 14

2

0

−2 100

0

(a)

0.5

80

0

−2

1

150 Frequency(Hz)

Frequency(Hz)

x6(t)

5

x (t)

2

0

(a)

0.5

1

(b)

Time(s)

0.5

1

(c)

Time(s)

0.5

1

0.5

1

100

60 40

50

20

(b)

0.5 Time(s)

1

150 Frequency(Hz)

Frequency(Hz)

100 80

100

60 40

50

20

(c)

0.5 Time(s)

1 150 Frequency(Hz)

Frequency(Hz)

100 80

100

60 40

50

20 0

0.5

1

Time(s)

Time(s)

Figure 7: Comparison between different versions of S-transform using s5 (t) (with t1 = 0.3, t2 = 0.7, t3 = 0.5, f 1 = 45, f 2 = 80 and f 3 = 15)(left) and s6 (t) (with a = 30, b = 40 and c = 150) (right). (a) Standard ST, (b) Moukadem’s ST, (c) CSK-ST.

leads to an energy concentration better than that obtained with the standard ST and comparable to that obtained with Moukadem’s ST. However, it should be noted that for the LFM monocomponent signal, the CSK-ST chows a better energy concentration than Moukadem’s ST, knowing that for a monocomponent LFM signal, performance analysis is usually defined in terms of the energy concentration [3]. Another important aspect in signal analysis is instantaneous frequency (IF) estimation. The LFM signal with transients s1 (t) is used. The IF estimation is based on the peaks values of local spectra. Figure 10 shows that the CSK-ST provides comparable IF estimation than Moukadem’ST. Compared to the standard ST, the CSK-ST shows notable improvement in IF estimation especially in noisy environment.

15

3

s5(t)+wn

s5(t)+wn

3

0

0

Frequency(Hz)

0.5

1

100

50

100

(b)

0.5

Time(s)

(c)

0.5 Time(s)

100

0.5

1

(b)

0.5 Time(s)

1

1

Frequency(Hz)

50

0.5 Time(s)

(a)

50

100

0

0

50

1

50

100 Frequency−Hz)

(a)

Frequency(Hz)

Frequency(Hz)

−3

0

Frequency(Hz)

−3

100

1

(c)

0.5

1

0.5

1

Time(s)

50

Time(s)

Figure 8: Comparison between different versions of S-transform using s5 (t) corrupted by an AWGN (5dB: left and 0dB: right). (a) Standard ST, (b) Moukadem’s ST, (c) CSK-ST.

3.3. Comparison with continuous wavelet transform In this test, the proposed approach is compared to the CWT when applied to s1 (t), s4 (t) and to a monocomponent LFM signal. For this purpose, the Morlet wavelet consistently used in the literature [18, 26], is used. It consists of a complex exponential modulated by a Gaussian. Its Fourier transform is a shifted Gaussian for each scale leading to a Bandpass filter bank analysis. The behavior of each method with a LFM signal is illustrated in Fig.11 where the CSK-ST is shown to have a frequency invariant amplitude response as the standard ST, in contrast to the CWT response where the amplitude is large for lower frequencies, and diminishes as the frequency of the chirp increases. In terms of resolution, the CSK-ST shows better frequency resolution all over 16

3

s (t)+wn 0.5

1

50

(b)

0.5

Time(s)

1

0.5

1

0.5

1

0.5

1

0.5 Time(s)

1

50

(b)

Time(s)

(c)

Time(s)

100

50

150

Frequency(Hz)

Frequency(Hz)

0.5

Time(s)

(a)

100

150

50

(c)

0

1

100

150

−3

150

Frequency(Hz)

Frequency(Hz)

(a)

100

150 Frequency(Hz)

0

Frequency(Hz)

−3

150

0

6

0

6

s (t)+wn

3

100

100

50

0.5 Time(s)

1

50

Figure 9: Comparison between different versions of S-transform using s6 (t) corrupted by an AWGN (5dB: left and 0dB: right). (a) Standard ST,(b) Moukadem’s ST, (c) CSK-ST.

the time-frequency plane. When the CWT is applied to multicomponent signals, the results indicated in Fig.12 show that the use of CWT for this kind of signals does not seem to be relevant compared to the CSK-ST. 3.4. Comparison with other time frequency representations 3.4.1. Test with linear and non linear FM signals To compare our approach to the Smoothed Pseudo WignerVille Distribution (SPWVD) [2, 4, 32] known to be optimal for LFM monocomponent signals [3, 4], and to the STFT, an application to s1 (t), s4 (t), s6 (t) and to crossing chirps s7 (t) is considered. Implemented in the function tfrspaw of the TF toolbox [32], the SPWVD locates in a perfect way all the components of the s1 (t) seen separately. Unfortunately, it still suffers from some interference problems at the transition points as shown in Fig.13 (right) and 17

Table 1: Energy concentration (CM) obtained with different versions of S-transform.

CM Standard ST Moukadem’s ST CSK-ST

s1 (t) s2 (t) s3 (t) s4 (t) s5 (t) s6 (t) LFM 0.0026 0.0018 0.0028 0.0043 0.0062 0.0043 0.0049 0.0044 0.0033 0.0051 0.0049 0.0078 0.0050 0.0050 0.0043 0.0032 0.0051 0.0048 0.0078 0.0049 0.0052

many cross terms which appear midway between true signal components of s7 (t) in Fig.13 (left). For the same realizations, the CSK-ST shows better time-frequency resolution. The components are well localized with no cross-terms. In Fig.14 the SPWVD shows some interference terms and bad time-frequency resolution for the sinusoidal FM component in s4 (t) and for the logarithmic FM component in s6 (t) especially at high frequencies. However, it still gives butter resolution for the linear chirps. For these examples, the CSK-ST gives a good time-frequency resolution for the linear components as for the non-linear components. Among these test, the STFT leads to results having poorer definition. Since the SPWVD is known for its best energy concentration, especially for LFM signals, a quantitative comparison is also made. Table 2 summarizing the energy concentration results CM shows that the CSK-ST provides a good compromise between components separability and concentration. Table 2: Energy concentration (CM) obtained with different time frequency representations.

signal s1 (t) s4 (t) s6 (t) s7 (t)

STFT 0.0037 0.0039 0.0043 0.0027

SPWVD 0.0044 0.0050 0.52 0.35

CSK-ST 0.0043 0.0048 0.0049 0.0032

3.4.2. Test with real ECG of pregnant women In this part of tests, a composite maternal ECG with 4 second length and sampled at 250Hz is considered. It is a low frequency signal where fetal ECG can also be seen. The goal is the fetal heartbeat detection. In such 18

Frequency(Hz)

Frequency(Hz)

400

(a)

400

300

300

200

200

100 400

100

(b)

0.5

1

300

200

200

100

Frequency(Hz)

400

300

400

(a)

(b)

0.5

1

(c)

0.5

1

0.5 Time(s)

1

100

(c)

0.5

1

400

300

300

200

200

100

100 0.5 Time(s)

1

Figure 10: Estimation of the instantaneous frequency of s1 (t). Ideal time-frequency representation (green), Estimated IF (Black). Free noise signal:(left), signal with high noise level, SN R = 0dB:(right). (a) Standard ST, (b) Moukadem’s ST, (c) CSK-ST.

application, a good time resolution is required. The CSK-ST is applied to the signal and compared to the STFT, the SPWVD and to the Standard ST. One can see in Fig.15 that the STFT gives a bad time-frequency resolution. The SPWVD gives a good time resolution for the R peaks of mother’s heartbeats, but the fetal heartbeats and the low frequency components of the ECG such as P and T waves are masked by the cross terms. The CSK-ST shows better time-frequency resolution than the standard ST. This can be especially seen at the first and the forth mothers heartbeats which are almost superimposed with the fetal ones. This test, present a new attempt to automatic detection of fetal heartbeats in the time-frequency plane. For this, the frequency range where fetal heartbeats are as relevant as mother’s ones, namely above 40Hz must be used. Accordingly, the already proposed method in [15] to heartbeat detection is used. This method is based on Shannon energy of normalized local spectra. n1 SSE(xj ) = − [Snorm (j, n)]2 log[Snorm (j, n)]2 . (25) n=no

19

(a)

2.5

(b)

400

2

300

1.5 200 1

Frequency(Hz)

Frequency(Hz)

400

300

200

100

100 0.5

0

0

0.5 Time(s)

1

0

0

0.5 Time(s)

Figure 11: Illustration of the CWT amplitude response (a) vs the CSK-ST response (b). Because of octave sampling of the scale axis in CWT, approximation is made to get the frequency axis.

where Snorm are the normalized local spectra and [no → n1 ] corresponds to the frequency range [40Hz, .., 80Hz]. Figure 16 shows that all fetal heartbeats have been detected although some of them are almost superimposed with the mother’s ones (the encircled ones). 4. Conclusion In this paper, a frequency dependent Compact Support Kernel (CSK) is introduced as an analysis window in the S-transform instead of the Gaussian analysis window used in the standard S-transform. The CSK is derived from the Gaussian window and retains its important properties. As in the standard S-transform, the width of the CSK is frequency dependent. It turns out that the control parameters of the kernel width make the CSK more flexible and more adaptive to the analyzed signal. These parameters are selected using an active set algorithm where the objective function is maximizing the energy concentration in time-frequency domain. The performance of the proposed scheme is examined using a set of synthetic test signals with: transients, linear and nonlinear FM laws, multi20

1

(a)

400

2.5

350

1.5

200 1 100

0.5

0 100

(a)

0.5 Time(s)

1

250 200

100

0

(b)

7

100

0.5 Time(s)

1

6

80 Frequency(Hz)

Frequency(Hz)

250

350

5 60

4

40

3

Frequency(Hz)

Frequency(Hz)

2

(b)

400

80 60 40

2 20

1 0

0.5

1

Time(s)

0

20 0

0

0.5 Time(s)

Figure 12: Representation of the CWT amplitude response (a) vs the CSK-ST response (b) using s1 (t) (Top) and s4 (t) (Bottom).

components and with white noise effect. The results of numerical analysis have shown that the proposed window enhance the time-frequency resolution of the standard ST and some existing modified ST. The comparisons made are not based solely on visual inspection of the time-frequency domain plots but also by examining the energy concentration especially for LFM signals. These results are mainly confirmed in presence of noise and when estimating the instantaneous frequency. Similar comparison is also made with other standard methods such as STFT, SPWVD and CWT. The achieved results show better resolution especially for signals with transitients, and with non linear FM laws. However, no method is superior to another, since they all have benefits for specific types of signals and applications. The performance of our method is also examined through a test performed with an ECG of 21

1

1

1

0

0.5

0

450

0.1

(b)

0.2

0.3

Normalized frequency

0

0.45 0.5

0

(b) 0.1

300

0.5

200

200

100

100

(c)

0.5

Time(s)

(c)

1

0.5

1

0.5 Time(s)

1

0.5

1

Time(s)

350

450

Frequency(Hz)

Frequency(Hz)

0.2 0.3 0.35 0.4 Normalized frequency

350 Frequency(Hz)

Frequency(Hz)

(a)

PSD

DSP

(a)

0.5

300 200

200 100

100

(d)

0.5

Time(s)

(d)

1 350

450 Frequency(Hz)

Frequency(Hz)

0

300 200 100 0.5 Time(s)

1

200 100

Time(s)

Figure 13: Time frequency representations using multicomponent LFM signals, s1 (t) (right) and s7 (t) (left). (a) Analyzed signal PSD, (b) SPWVD, (c) STFT and (d) CSK-ST.

a pregnant woman for detecting fetal heartbeats. The KSC-ST may be useful in ECG segmentation since it presents good time resolution even at low frequencies. References [1] D.Gabor. Theory of communications. J.Inst.Elec.Eng. 93, pp.429-457, (1946). [2] L. Cohen. Time-Frequency Analysis. Prentice-Hall PTR, New York, (1995). [3] B. Boashash and V. Sucic. Resolution measure criteria for the objective 22

(a)

Frequency(Hz)

Frequency(Hz)

80

100

60 40 20

(b)

0.5

1

Time(s)

150 Frequency(Hz)

Frequency(Hz)

100

50

0

0

80

0.5

1

(b)

Time(s)

(c)

0.5 Time(s)

1

0.5

1

100

60 40

50

20

(c)

0.5 Time(s)

0

1

150 Frequency(Hz)

100 Frequency(Hz)

(a)

150

100

80

100

60 40

50

20 0

0

0.5 Time(s)

1

Time(s)

Figure 14: Time frequency representations using s4 (t) (left) and s6 (t)(right). (a)SPWVD, (b) STFT, (c) CSK-ST.

assessment of the performance of quadratic time-frequency distributions. IEEE.Trans.Signal.Process, 51(5), pp. 1253-1263, (2003). [4] P.Flandrin, F. Auger and E. Chassande-Mottin. Applications in TimeFrequency signal proceessing. Time-frequency reassignment: from principles to algorithms. Taylor and Francis Group, (2003). [5] P.Goupillaud, A.Grossmann and J.Morlet. Cycle-octave and related transforms in seismic analysis. Geoexploration 23, pp.84–102, (1984). [6] I.Daubechies. The wavelet transform, time-frequency localization and signal analysis. IEEE.Trans.Inf.Theory, 36(5), pp.961-1005, (1990).

23

(a)

1 ECG

C

M F

−1

0

Frequency(Hz)

125

Frequency(Hz)

F

M

F

F M

F

M

F

F

(b)

1

(c)

1

2

3

4

2

3

4

40 0

125

Frequency(Hz)

F

40

0 0 125

(d)

1

2

3

4

(e)

1

2

3

4

2

3

40

125

Frequency(Hz)

M

F

0

40 0

1

Time(s)

4

Figure 15: Comparison between different time-frequency representations applied to a pregnant woman ECG. (a) ECG composite (M:Mother’s beat, F:fetal beat), (b) STFT, (c) SPWVD, (d) Standard ST, (e)Proposed CSK-ST.

[7] Z.Zidelmal, A.Amirou, M.Adnane and A.Belouchrani. QRS detection using wavelet coefficients. Computer Methods and Programs in Biomedicine, 107(3), pp.490–496, (2012). [8] E.Castillo, D.P.Morales, G.Botella, A.Garcia, L.Parrilla and A.J.Palma. Efficient wavelet-based ECG processing for single-lead FHR extractionOriginal. Digital Signal Processing, 23(6), pp.1897–1909, (2013). [9] R.G.Stockwell, L.Mansinha and R.P.Lowe. Localisation of the complex spectrum: The S-Transform. IEEE Trans.signal.process, 44 (4), 998– 1001, (1996).

24

1

(a) ECGC

0.5 0

−0.5 −1

0

1

2

3

4

1

(b) SSEn

0.8 0.6 0.4 0.2 0

1

2

Time(s)

3

4

Figure 16: Shannon energy showing the candidate fetal heartbeats (b) and the detected fetal heartbeats on Composite ECG (a).

[10] R.G.Stockwell. A basic for efficient representation of the S-transform. Digital signal processing, 17(1), pp.371–393, (2007). [11] S.Ventosa, C.Simon, M.Schimmel, J.J.Danobeitia, and A.Manuel. The S-Transform From a Wavelet Point of View. IEEE.Trans.Signal.process, 56(7), pp.2771–2780, (2008). [12] A.Moukadem, D.Ould Abdeslam and A.Dieterlen. Time-Frequency Domain for Segmentation and Classification of Non-Stationary Signals. Wiley-ISTE, ISBN 978-1-84821-613-6, (2014). [13] A.Amirou, D.Ould Abdeslam, Z.Zidelmal, M.Aiden and J.Merckle. STransform and Shannon energy for electrical disturbances detection. IECON’2014, the 40th Annual Conference of the IEEE Industrial Electronics Society, Dallas, TX-USA Oct 29–Nov 1, (2014). [14] M.J.Bharata.Reddya, R.Krishnan.Raghupathya, K.P.Venkatesha and D.K.Mohanta. Power quality analysis using Discrete Orthogonal Stransform (DOST). Digital Signal Processing, 23(2), pp.616-626,(2013). [15] Z.Zidelmal, A.Amirou, D.Ould Abdeslam, A.Moukadem and A.Dieterlin. QRS detection using S-Transform and Shannon Energy. Computer Methods and Programs in Biomedicine, 116(1), pp.1–9, (2014).

25

[16] R.A.Brown, M.Louis.Lauzon and R.Frayne. A General Description of Linear Time-Frequency Transforms and formulation of a Fast, Invertible Transform That Samples the Continuous S-Transform Spectrum Nonredundantly. IEEE.Tran.Signal.Process, 58(1), pp.281–289, (2010). [17] R.G.Stockwell. S-Transform Analysis of Gravity Wave Activity. Ph.D. Dissertation, Dept. of Physics and Astronomy, The University of Western Ontario, Canada, (1999). [18] R.G.Stockwell Why use the S-Transform ?. AMS Pseudodifferential operators: partial differential equations and timefrequency analysis, 52, 279–309, (2007). [19] P.Flandrin. Time-frequency/time-scale analysis. San Diego, Academic Press (1999) [20] C.R.Pinnegar and L.Mansinha. The S-transform with windows of arbitrary and varying shape. Geophysics, 68(1), 381-385, 2003. [21] C.R.Pinnegar and L.Mansinha. The Bi-Gaussian S-transform. SIAM J.Sci.Comput 24(5), 1678–1692, (2003). [22] K.Kazemi, M.Amirian and M.J.Dehghani. The S-Transform using a new window to improve frequency and time resolution. Signal, Image and Video Processing (SIVP), 8(3), pp.533–541, (2013). [23] E.Sejdic, I.Djurovic and J.Jiang. S-Transform with frequency dependent kaiser window. IEEE Explore, Acoustics, Speech and Signal Processing, (ICASSP), (2007). [24] E.Sejdic, I.Djurovic, J.Jiang. A window widh optimized S-transform. EURASIP J.Adv.Signal.Process, (2008). [25] I.Djurovic, E.Sejdic, J.Jiang. Frequency-based window widh optimization for S-transform. Int.J.Electron.Commun, 62(4), pp.245–250, (2008). [26] S.Assous and B.Boashash. Evaluation of the modified S-transform for time-frequency synchrony analysis and source localisation. EURASIP J.Adv.Signal.Process, (49), pp.1–18, (2012).

26

[27] A.Moukadem, Z.Bouguila, D.Ould-Abdeslam, A.Dieterlen. A new optimized Steckwell transform applied on synthetic and non stationary signals. Digital Signal processing, 46, pp.226–238, (2015). [28] S. Saryazdi and M. Cheriet. PKCS: A polynomial kernel family with compact support for scale-space image processing. IEEE Trans.Image.Process. 16(9), pp.2299–2307, (2007). [29] A. Belouchrani and M. Cheriet. On the use of a new compact support kernel in time-frequency analysis. in Proc. 11th IEEE Workshop on Statist. Signal Process, Singapore, pp.333-336, May, (2001). [30] M.Abed, A.Belouchrani, M.Cheriet and B.Boashash. Time-Frequency Distributions Based on Compact Support Kernels: Properties and Performance Evaluation. IEEE Tans.Signal.Process, 60(6), pp.2814–2827, (2012). [31] B.Boashash. Time-Frequency Signal Analysis and Processing: A Comprehensive Reference. Ed. Amsterdam, The Netherlands: Elsevier, (2003). [32] F.Auger, P.Flandrin, P.Gon¸calvˇcs and O.Lemoine. Time-frequency toolbox. CNRS France, Rice University, (1996) [33] J.Nocedal and S.Wright. Numerical Optimization. (2nd ed.), New York: Springer-Verlag, ISBN 978-0-387-30303-1, (2006).

27

Zahia Zidelmal received the State Engineering degree in 1990 and the M.S. degree in signal processing and biomedical engineering in 1996 both from University Mouloud Mammeri of Tizi-Ouzou (UMMTO), Algeria. Since 1993, she has been with the Electrical Engineering Department of UMMTO first as Teaching Assistant and then as Teaching Assistant Fellow since 1996. She is currently and since 2012 associate professor at UMMTO. Her research interests are in statistical signal processing with applications in biomedical and communications, time-frequency representations and pattern recognition with Support Vector Machines (SVMs). Hocine Hamil received the M.Sc degree in computer engineering in 2014 from University Mouloud Mammeri of Tizi-Ouzou (UMMTO). Currently he is working towards her Ph.D. Thesis. His research interests are time-frequency array signal processing with applications in biomedical and FPGA implementation of signal processing algorithms. Ali Moukadem received the M.Sc. Degree in Software Engineering from the University of Technology of Belfort-Montbéliard, France in 2007, the M.Sc. and the Ph.D. degrees in signal processing form the University of Haute Alsace, France in 2008 and 2011, respectively. The main research interests of Dr. Moukadem are: time–frequency analysis and methods, nonstationary signal processing, applied mathematics and pattern recognition. He is currently an Associate Professor in the University of Haute Alsace. Ahmed Amirou received the State Engineering degree in computer engineering in 1986 and the M.Sc. degree in Operational Research in 2006 both from University Mouloud Mammeri of Tizi-Ouzou (UMMTO), Algeria. Since 1990, he has been with the Electrical Engineering department of UMMTO first, as Teaching Assistant and then as Teaching Assistant Fellow since 2006. He is currently and since 2015 associate professor at UMMTO. His research interests are in Operational Research and Global Optimization. Djaffar Ould Abdeslam received the M.Sc. degree in electrical engineering from the University of Franche-Comté, Besançon, France, in 2002, the Ph.D. degree from the University of Haute-Alsace, Mulhouse, France, in 2005. From 2007, he is an Associate Professor at the University of Haute-Alsace. He received the Habilitation in Electrical Engineering from the University of Haute-Alsace in 2014. His research interest includes Artificial Neural Networks applied to Power Systems, Signal Processing for Power Quality Improvement, Smart Metering, Smart Building and Smart Grids.