Signal Processing 92 (2012) 2264–2274

Contents lists available at SciVerse ScienceDirect

Signal Processing
journal homepage: [Link]/locate/sigpro

A generalized synchrosqueezing transform for enhancing signal

time–frequency representation
Chuan Li a,b, Ming Liang b,n
a
Engineering Laboratory for Detection, Control and Integrated System, Chongqing Technology and Business University, Chongqing 400067, China
b
Department of Mechanical Engineering, University of Ottawa, Ottawa, Canada K1N 6N5

a r t i c l e in f o abstract

Article history: High-quality time–frequency representation (TFR) is important for reliable signal analysis.
Received 13 May 2011 The diffusions of the TFR energy along time and/or frequency axes lead to ambiguous TFR
Received in revised form and hence misleading signal analysis results. Synchrosqueezing is an adaptive and
10 October 2011
invertible transform developed to improve the quality or readability of the wavelet-based
Accepted 21 February 2012
TFR by condensing it along the frequency axis. However, the original synchrosqueezing
Available online 1 March 2012
method could be handicapped by time-dimension diffusions of the wavelet coefficients.
Keywords: As such, we propose a generalized synchrosqueezing transform (GST) approach to deal with
Generalized synchrosqueezing transform the diffusions in both time and frequency dimensions. For the signal with a constant
Time–frequency representation
frequency, we have shown that the wavelet diffusion only occurs at frequency dimension.
Wavelet transform
Based on this observation, the original signal with time-varying instantaneous frequency is
Instantaneous frequency
Instantaneous amplitude mapped to another analytical signal with constant frequency to facilitate the synchros-
queezing. A time-scale domain restoration operation is then presented to obtain a TFR with
concentrated wavelet ridge. The performance of the proposed GST for signal TFR enhance-
ment has been demonstrated by our simulation study.
& 2012 Elsevier B.V. All rights reserved.

1. Introduction and wavelet transform (WT) [9] are three popular tools
for TFR. However, there are some inherent limitations
The time or frequency domain description alone may with these TFR tools. For example, the WVD introduces
not provide adequate information for a non-stationary cross-terms between multiple signal components and
signal. One of the most effective tools that helps to as such, it is very difficult to extract the instantaneous
understand such a signal is the time–frequency represen- frequency (IF) from the WVD-based TFR. For STFT, the
tation (TFR) which displays frequency components and invariable window width leads to poor time–frequency
their amplitudes occurring at any given time points [1]. (TF) resolution. Much effort has been taken to improve
The TFR has been successfully applied in a wide range of the TF resolution of different TFRs of a signal [10–12].
fields such as speech, mechanical, biomedicine, seismic However, the TF resolution is limited by the Heisenberg
and radar signals processing [2–6]. Short time Fourier uncertainty principle [13], i.e. Dt Df Z1=4p where Dt
transform (STFT) [7], Wigner–Ville distribution (WVD) [8] represents the time resolution and Df stands for the
frequency resolution. The minimal value of Dt Df is known
as the Heisenberg box. It is obvious that the increase of
n
the time resolution leads to the drop of the frequency
Corresponding author. Tel.: þ 1 613 562 5800x6269;
fax: þ1 613 562 5177.
resolution, and vice versa.
E-mail addresses: chuanli@[Link] (C. Li), It is well known that the wavelet transform can be
liang@[Link] (M. Liang). applied to extract the IF by detecting the ridge of the WT

0165-1684/$ - see front matter & 2012 Elsevier B.V. All rights reserved.
doi:10.1016/[Link].2012.02.019
 C. Li, M. Liang / Signal Processing 92 (2012) 2264–2274 2265

[14]. However, the WT also yields blurred TFR and thus The observations made in this process inspire us to
reduces the readability of the instantaneous amplitude develop the GST approach.
(IA) from the TFR. Daubechies et al. [15] proposed a
wavelet-based TFR enhancement method which was 2.1. The complex continuous wavelet transform for TFR
called ‘‘synchrosqueezing’’. This method has been applied
to TF analysis of electrocardiography and paleoclimatol- A monocomponent asymptotic AM-FM signal is given by
ogy signals [16,17].
sðtÞ ¼ AðtÞ cosðfðtÞÞ ð1Þ
The synchrosqueezing approach focuses the TFR, Rðt,f Þ,
of the signal s(t) by squeezing its value along the frequency f where Ainst(t)¼A(t) is the IA and fðtÞ the instantaneous
(scale a) axis. Compared to other TFR enhancement methods phase (IP) of the signal. The IF of the signal is accordingly
such as reassignment [18], synchrosqueezing offers better defined as
adaptability and an exact reconstruction formula for con-
dfðtÞ
stituent components. As pointed out by Daubechies and f inst ðtÞ ¼ ð2Þ
2p dt
Maes [19], however, this transform is limited to scale
variable of the continuous wavelet transform (CWT). In fact, The complex CWT of the signal s(t) is given by
Z  
there is also a time offset b (corresponding to time t)   1 tb
variable in addition to the scale for the CWT [20]. Though W s ða,bÞ ¼ sðtÞ, ca,b ðtÞ ¼ pffiffiffi sðtÞcn dt ð3Þ
a a
the a-dimension spread is dominant for an ideal signal,
b-dimension diffusion cannot be overlooked when dealing where symbol ‘‘n’’ denotes the complex conjugate opera-
with noisy, amplitude modulated-frequency modulated tion, c is a properly selected complex wavelet that is
(AM-FM) signals. Even if the IF trajectory of an AM-FM concentrated on the positive frequency axis, i.e. j
^ ðxÞ ¼ 0
signal can be extracted by the synchrosqueezing algorithm, for x o 0. cðtÞ can be expressed as
the related IA remains questionable due to the b-dimension cðtÞ ¼ gðtÞ expðio0 tÞ ð4Þ
diffusion. In addition, the TF resolution of the synchros-
queezing, as a post-processing approach for TFR enhance- where g(t) is the window function and o0 the central
ment, is still limited as stated by the Heisenberg uncertainty angular frequency of the wavelet.
principle. Along the scale direction, a maximum value of the
As such, we propose a simple yet flexible approach modulus of W s ða,bÞ exists at each time point. A curve
which may be called generalized synchrosqueezing trans- composed of these maximum values is known as wavelet
form (GST) to resolve the bi-dimensional (a and b) smear ridge, or simply, ridge, that can be expressed by a set as
problem. Instead of a direct bi-dimensional squeezing P ¼ fðar ,bÞ 2 R2 ; M s ðar ,bÞ ¼ maxð9W s ðai ,bÞ9Þg ð5Þ
transform, however, the proposed GST paves the way
towards much more condensed TFR by a combination of where (ar, b) is the ridge point at any time point b,
multiple TF plane manipulations (such as translation, 9W s ða,bÞ9 denotes the modulus of the wavelet coefficient,
rotation, and twisting) in one operation. We will demon- i.e.
strate that the diffusion only occurs at a-dimension for   qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
W s ða,bÞ ¼ ðReðW s ða,bÞÞÞ2 þ ðlmðW s ða,bÞÞÞ2 ð6Þ
the TFR of the signal with a constant frequency. This
observation inspires us to map a time-varying frequency where ReðW s ða,bÞÞ and lmðW s ða,bÞÞ respectively represent
signal to another corresponding signal with constant the real part and the imaginary part of the WT of signal s(t).
frequency to facilitate the synchrosqueezing. A time-scale The analytical expression of the signal s(t) as shown in
domain restoration operation is presented to obtain a TFR (1) is often unavailable in real applications. Through
with concentrated wavelet ridge from where both IF and performing the WT, fortunately, the IA and the IF of
IA can be extracted more accurately. Since the proposed the signal can be calculated from the extracted ridge (ar, b)
algorithm is no longer a simple post-processing techni- as [21]
que, the TF resolution dilemma is also alleviated during o0
the mapping and the restoration operations. f inst ðtÞ ¼ ð7Þ
2par ðtÞ
This paper is organized as follows. In Section 2, we
provide a brief introduction to the synchrosqueezed WT. 29W s ðar ðtÞ,tÞ9
Some observations of synchrosqueezed wavelet coeffi- Ainst ðtÞ  pffiffiffiffiffiffiffiffiffiffi ð8Þ
ar ðtÞ9gð0Þ9
^
cient modulus and chirp rate are also made in this section.
The GST approach is then detailed in Section 3. Section 4 ^
where gð0Þ ^ oÞ9o ¼ 0 , gð
¼ gð ^ oÞ denotes the Fourier trans-
presents a simulation case study and discussions. The form of g(t).
conclusions are given in Section 5. The TF information of a signal can be accordingly
extracted using the WT as described by (7) and (8).
2. TFR using synchrosqueezed WT However, the WT often produces a blurred TF picture in
terms of finst(t) and/or Ainst(t). This may cause misinter-
In this section, we first briefly introduce the complex pretation of the signal. To mitigate the TFR blurring effect
continuous wavelet transform used for TFR of the signal. stemming from the WT, Daubechies et al. [15] proposed
The synchrosqueezing method is then presented to the synchrosqueezing algorithm which will be briefly
enhance the WT-based TFR. The relationship between introduced with our observations in the following
the synchrosqueezing effect and the chirp rate is analyzed. subsection.
2266 C. Li, M. Liang / Signal Processing 92 (2012) 2264–2274

2.2. The synchrosqueezing algorithm and observations Such deficiency leads to substantially distorted IF as com-
pared to the true IF shown in Fig. 1(b). Although the IF of
Suppose that the signal has a constant IF and a region B may be vaguely observed from the a-dimensional
constant IA, i.e. synchrosqueezed TFR (Fig. 1(d)), the extracted IA trajectory
f inst ðtÞ ¼ f , Ainst ðtÞ ¼ A ð9Þ of region B is unacceptable (see Fig. 1(e)). It is clear that the
higher the 9c(t)9, the worse the synchrosqueezing effect is.
Substituting (9) into (3) yields Comparing regions A and B of Fig. 1(d) and (e) reveals that,
Z
1 A for a horizontal wavelet ridge (c(t)¼0), a frequency-axis
^ ðaxÞ eibx dx ¼ pffiffiffi j
W s ða,bÞ ¼ pffiffiffi s^ ðxÞj ^ ðaf Þ eibf ð10Þ
a 2 a synchrosqueezing step should be sufficient but for a curved
ridge (c(t)a0) there will be a substantial amount of smear
W s ða,bÞ often spreads out and leads to a smear projection
‘‘residual’’ after such uni-axial synchrosqueezing.
in the time-scale plane. For a given time offset b, this
To overcome the aforementioned shortcomings, one
diffusion occurs mainly in a-dimension. If the smear
may subsequently suggest a bi-axial synchrosqueezing
behavior in b-dimension is neglected, at any (a, b) location
approach by taking b into account in (11) and (12).
for which W s ða,bÞa0, its IF f inst ða,bÞ can be computed by
However, it may not be feasible for such a reassignment-
i @ðW s ða,bÞÞ like synchrosqueezing to find a solution. Even if it can find
f inst ða,bÞ ¼ ¼ ð11Þ
W s ða,bÞ @b a solution, it would be computationally burdensome.
The above observations motivate us to map a signal
In actual computations, a, b and f are all discrete.
whose c(t)a0 to another corresponding signal with
Suppose ðDaÞk ¼ ak ak1 for any ak. When mapping
c(t)¼ 0 to facilitate synchrosqueezing. As such, we pro-
from the time-scale plane to the TF plane with ðb,aÞ-
pose a generalized synchrosqueezing approach that is
ðb,f inst ða,bÞÞ, the synchrosqueezing transform Sðf ,bÞ is per-
simple to implement, flexible in signal mapping and
formed at the center fl of the frequency range ½f l Df =2,
robust to noise with fine resolution.
f l þ Df =2 (where Df ¼ f l f l1 ):
1 X 3=2
Sðf l ,bÞ ¼ W s ðak ,bÞak ðDaÞk ð12Þ 3. The generalized synchrosqueezing transform
Df Df
ak :9f s ðak ,bÞf l 9 r 2

The above equation shows that the TFR of the signal To facilitate the synchrosqueezing of an FM signal, the
s(t) is synchrosqueezed along the frequency (or scale a) IF trajectory mapping is first introduced in this section.
axis only. Let us employ the chirp rate c(t) to represent The GST approach is then presented, followed by compar-
the rate of change of the IF, i.e., [22] isons with several existing TFR tools to validate the
effectiveness of the proposed GST approach.
cðtÞ ¼ df inst ðtÞ=dt ð13Þ
If c(t)a0, the diffusion occurs in both time and frequency
3.1. IF trajectory mapping and restoration for generalized
dimensions. As a result, a uni-axial synchrosqueezing as
synchrosqueezing
described by (12) is inadequate for TFR enhancement.
This can be illustrated by an example as shown in Fig. 1
Fig. 1 shows that there is no blur at the time-dimen-
where a signal s(t) is defined by
( sion if c(t) ¼0. It is therefore necessary to squeeze the
0:8 cosðpt=15Þ 0 r t r600 horizontal ridge only in the frequency dimension. This
sðtÞ ¼
cosðpt=10 þ 2p sinðpt=150ÞÞ 600 o t r 1024 motivates us to map the signal from c(t)a0 to c(t) ¼0 (or
cðtÞ-0) such that the time-dimension spread can be
ð14Þ
avoided or at least alleviated. The generalized Fourier
Fig. 1(a) shows the temporal waveform of s(t) whose transform (GFT) [23] is employed to illustrate our idea.
ideal IF is plotted in Fig. 1(b). Using the Morlet CWT, a The GFT of the signal s(t) is given by
sketchy frequency trajectory is obtained as shown in Z þ1
Fig. 1(c). From this figure, we can see that the WT-based SG ðf Þ ¼ FG ðsðtÞÞ ¼ sðtÞ ei2pðf t þ x0 ðtÞÞ dt ð15Þ
1
TFR is too vague for us to comprehend the TF information of
the signal. To get a clear projection, the synchrosqueezing where x0(t) is a real-valued transform function specifying
algorithm is employed to sharpen Fig. 1(c) and the shar- the phase transform behavior of the signal. The signal
pened result is shown in Fig. 1(d). Along the wavelet ridge in may be also computed from the inverse GFT, i.e.
the TFR of Fig. 1(d) and (e) shows the modulus maxima Z þ1
i2px0 ðtÞ
MOD(t) corresponding to the IA as described by (8). sðtÞ ¼ F1
G ðSG ðf ÞÞ ¼ e SG ðf Þ ei2pf 0 df ð16Þ
1
Comparing Fig. 1(b) and (c) indicates that the WT
method generates heavily blurred TFR, even for the
horizontal portion. Hence, synchrosqueezing is applied If SG ðf Þ  dðf f 0 Þ, then sðtÞ ¼ ei2pðf 0 t þ x0 ðtÞÞ where f0 is the
to get much better result shown in Fig. 1(d). Region A is desired horizontal GFT frequency trajectory expressed as
part of a horizontal wavelet ridge (c(t)¼0) where both the
dx0 ðtÞ
IF and IA trajectories can be accurately extracted from the f 0 ¼ f ðtÞ ¼ f ðtÞx00 ðtÞ ð17Þ
dt
synchrosqueezed WT-based TFR. On the contrary, region B
is located in a curved ridge (c(t)a0) for which synchros- Since both the GFT and the CWT produce the same
queezing is less effective than for the horizontal ridge. (in an ideal situation) IF trajectory, the horizontal GFT
 C. Li, M. Liang / Signal Processing 92 (2012) 2264–2274 2267

1 0.1

0.08

Frequency (Hz)
0.5
Amplitude
0.06
0
0.04
-0.5
0.02

-1 0
0 200 400 600 800 1000 0 200 400 600 800 1000
Time (s) Time (s)

0.1 0.1 250

70
0.08 60 0.08 200
Frequency (Hz)

Frequency (Hz)
50
0.06 0.06 150
40
B
0.04 30 0.04 100

0.02 20 0.02 A
50
10
0 0
200 400 600 800 1000 200 400 600 800 1000
Time (s) Time (s)

260
240
Maximal coefficients

220 A
200
180
160
B
140
120
0 200 400 600 800 1000
Time (s)

Fig. 1. Observations of the synchrosqueezing algorithm illustrated using an example: (a) temporal waveform of s(t) given by (14) with sampling
frequency of 1 Hz, (b) real IF trajectory of the signal, (c) TFR of s(t) produced by the Morlet CWT, (d) synchrosqueezed result of (c), and (e) modulus
maxima MOD(t) (corresponding to the IA) of the wavelet ridge from the result of (d). (Note: The TFR coefficients shown in the subplots, i.e., the vertical
bars, of (c) and (d) are not in the same scale. Such values from different methods are not comparable.)

frequency ridge is associated with a horizontal IF for the Hilbert transform is applied to create a new analytic
same signal. According to the GFT theory, the energy of signal that is free from the above influence
the signal s(t) with IF f(t) can converge to a horizontal
wðtÞ ¼ vðtÞ þ iHðvðtÞÞ ð21Þ
trajectory f0 (i.e. c(t)¼0) using a mapping function ei2px0 ðtÞ
(i.e., sðtÞ-sðtÞei2px0 ðtÞ ). Similarly, applying this mapping This analytic signal will replace the original s(t) to
function to the frequency representation (i.e., SG ðf 0 Þ- produce CWT Ww. This time-scale result is then trans-
SG ðf 0 Þei2px0 ðtÞ ) will restore the time-domain signal. This formed as
process is partly motivated by the generalized demodula-
tion approach proposed by Olhede and Walden [24]. W y ða,bÞ ¼ W w ða,bÞei2px0 ðtÞ ð22Þ
To apply the IF trajectory mapping in the complex It should be noted that the modulus of Wy is equal to
domain, s(t) is first transformed to an analytic signal that of Ww, with respectively mapped real and imaginary
uðtÞ ¼ sðtÞ þ iHðsðtÞÞ ð18Þ parts. This is the main reason why the GST approach can
enhance TFR of the signal. According to (11), the joint
where HðsðtÞÞ represents the Hilbert transform of s(t), i.e.,
real part and imaginary part mapping restores the IF
HðsðtÞÞ ¼ sðtÞnð1=ðptÞÞ ð19Þ f yinst ða,bÞ ¼ f sinst ða,bÞ. However, the synchrosqueezing
The analytic signal u(t) is then transformed with the transform defined by (12) shows that the TF convergence
mapping function ei2px0 ðtÞ as feature of analytic signal w(t) that has frequency f0 can be
preserved by the IF trajectory mapping, i.e.
vðtÞ ¼ uðtÞei2px0 ðtÞ ð20Þ  
To eliminate the influence of negative frequency com- f ðtÞx00 ðtÞ
Ss ðf s ,bÞ3Sy ðf y ,bÞ ¼ Sw f w ,b ð23Þ
ponents in v(t) during the IF trajectory mapping, another f ðtÞ
2268 C. Li, M. Liang / Signal Processing 92 (2012) 2264–2274

The above equation indicates that the IF of s(t) has two purposes of the IF trajectory mapping and restora-
been restored and the TFR has been further condensed by tion: (1) to facilitate the synchrosqueezing for a time-
the proposed GST. In this way, an enhanced TFR may be varying frequency signal regardless of the influence of the
obtained by eliminating the diffusion in both time and time-dimension TFR diffusion; and (2) to alleviate the TF
frequency dimensions. Hence, both the IF and the IA resolution dilemma. The first purpose has been described
trajectories can be accordingly extracted from the GST in Section 3.1. Here, we elaborate the second purpose.
enhanced TFR. As illustrated in Section 1, the TF resolution of the CWT
According to the above introduction, the GST approach is limited by the Heisenberg box. Within the lower scales,
can be implemented following the flowchart illustrated the CWT features higher frequency resolution and lower
in Fig. 2. time resolution. On the contrary, the increase of the scale
With the TFR obtained by the GST, the signal can be results in the improvement of the time resolution and the
reconstructed in a way analogous to that used by the deterioration of the frequency resolution. When imple-
original synchrosqueezing method. The original synchros- menting the GST, fortunately, two separate steps are
queezing reconstructs a signal based on the following applied to improve the frequency resolution without
equation [15]: reducing the time resolution: (1) choose a small constant
! f0 as the target frequency so as to obtain higher frequency
X
sðbÞ ¼ Re C 1
j S ðf
s l ,bÞðDf Þ ð24Þ resolution; and (2) perform the IF restoration operation to
l restore the original time resolution. In the aforementioned
R þ1 first step, the rise of the frequency resolution may lead to
where C j ¼ 0:5 0 ðj ^ ðxÞ=xÞ dx. Replacing Ss ðf l ,bÞ with
the drop of the time resolution according to the Heisenberg
Sy ðf y ,bÞ in the above equation leads to the exact signal
box. Thankfully, the horizontal ridge is a special case. Since
reconstruction formula of the GST:
X an ideal horizontal ridge is time-independent (f0(t)f0),
sðbÞ ¼ ReðC 1j Sy ðf y ,bÞðDf ÞÞ ð25Þ the drop of time resolution do not affect this IF trajectory
mapping. In other words, the original time resolution is
preserved with increased frequency resolution for the GST.
3.2. GST performance in comparison with existing TFR This will be further demonstrated using the simulation
methods case study as shown in Section 4.
According to the analysis above, usually a low fre-
It should be noted that desired horizontal GFT fre- quency horizontal ridge (constant frequency) can be
quency trajectory (i.e. constant frequency f0) is an impor- chosen as f0. If a true horizontal ridge is not available in
tant parameter in our TFR approach. There are, in general, real applications, a quasi-horizontal ridge is also accep-
table for the proposed GST method.
To examine the effectiveness of the proposed method,
the signal s(t) given by (14) is mapped to the analytic
signal w(t) with target frequency f0 ¼ 1/30 Hz. According
to (17), we have
(
1 0 r t r600
x0 ðtÞ ¼ ð26Þ
5t þsinð2ptÞ 600 o t r1024

The TFR of the associated analytic signal w(t) is then

restored by ei2px0 ðtÞ for GST. Fig. 3(a) shows the enhanced
TFR produced by the GST approach, which is very close to
the true TFR shown in Fig. 1(b). Again, the modulus
maxima MOD(t) along the wavelet ridge are employed
to reveal the sharpness of the TFR. Fig. 3(b) plots the
modulus maxima MOD(t) of the TFR using the GST
approach. The reference regions A and B are also high-
lighted in Fig. 3(a) and (b). The enhancement can be
clearly observed by comparing regions A and B in Fig. 3
with their counterparts obtained using CWT (Fig. 1(c))
and synchrosqueezed CWT (Fig. 1(d)).
According to (8), the IA Ainst ðtÞ can be extracted from
the MOD(t), i.e.

MODðtÞ2L þ 1
Ainst ðtÞ ¼ ð27Þ
2L þ 1
where L stands for the number of the WT scales. In
our demonstration, L¼6. Fig. 3(c) shows the comparison
between the real and the extracted IA trajectories (calcu-
Fig. 2. Flowchart of the proposed GST approach for signal TFR. lated using (27)). It is clear that, in addition to the IF
 C. Li, M. Liang / Signal Processing 92 (2012) 2264–2274 2269

0.1 250 260

0.08 B
200 240

Maximal coefficients
Frequency (Hz)

0.06 150 220 A

B
0.04 100 200

A
0.02 50 180

0 160
200 400 600 800 1000 0 200 400 600 800 1000
Time (s) Time (s)

1.5
1 Original signal
1 Reconstructed signal
Instantueous amplitude

0.8
0.5

Amplitude
0.6
0

0.4 -0.5
Extracted IA
0.2 Real IA -1

0 -1.5
0 200 400 600 800 1000 0 100 200 300 400 500 600 700 800 900 1000
Time (s) Time (s)

0.2
0.1
0
-0.1
Error

-0.2
-0.3
-0.4
-0.5
-0.6
0 100 200 300 400 500 600 700 800 900 1000
Time (s)

Fig. 3. Demonstration of the GST performance: (a) TFR obtained by the GST, (b) modulus maxima of the IF trajectory, (c) IA trajectory estimated from
MOD(t), (d) signal reconstruction from the GST-based TFR as shown in (a), and (e) the reconstruction error.

trajectory enhancement, the IA trajectory can also be smoothed pseudo Wigner–Ville distribution (SPWVD) [25],
estimated by the proposed GST approach (terminal effects reassigned SPWVD and reassigned Morlet CWT as follows.
are not taken into consideration). One may better appreci- Fig. 4(a)–(c) show the TFRs of the signal s(t) given by (14)
ate the effectiveness of the GST by comparing Fig. 3(c) using the SPWVD, the reassigned SPWVD and the reas-
with Fig. 1(e). Fig. 3(d) shows the comparison between signed Morlet CWT, respectively. It is evident that the
the original signal and the signal reconstructed from the assignment method is also effective to condense the hor-
GST-based TFR shown in Fig. 3(a). The reconstruction izontal ridge (i.e. c(t)¼0). When dealing with a time-varying
error (the difference between the original and the recon- IF trajectory (c(t)a0), however, there are apparent trajec-
structed signals) is plotted in Fig. 3(e) which is quite close tory distortions in the TFRs obtained by the bi-axial assign-
to zero indicating good reconstruction accuracy. ment approaches, i.e., the reassigned SPWVD and the
The comparison between Figs. 1 and 3 indicates that the reassigned CWT. Comparing Fig. 3 with Fig. 4 reveals that
proposed GST yields better TFR than that resulting from the the proposed GST performs much better than the existing
original synchrosqueezing in terms of both IA and IF. Besides TFR tools. More importantly, the present GST does not suffer
synchrosqueezing, reassignment is another popular TFR from the heavy computation burden present in the bi-axial
enhancement approach. To further assess the performance reassignment methods. In comparison to the original syn-
of the proposed GST method, we will then compare it with chrosqueezing, the additional computation required by the
2270 C. Li, M. Liang / Signal Processing 92 (2012) 2264–2274

TFR. When the signal becomes longer, however, more

0.1 22
computing time is required by the synchrosqueezing or
the GST or any other methods. In addition, in the absence
0.08 18
of prior knowledge, one should calculate the mapping
Frequency (Hz)

function from (17) by estimating f(t) from the standard

0.06 14
synchrosqueezing TFR. This estimation step hinders the
real time application of the proposed GST algorithm.
0.04 10
Hence, the proposed GST algorithm can be applicable in
real time situations if the sampling frequency is relatively
6
0.02 low and the mapping function is known beforehand. In
summary, the GST is faster than other popular reassign-
2
ment methods. Even so, it may not be fast enough for real
0 200 400 600 800 1000 time applications mainly due to a separate step to esti-
Time (s) mate the mapping function (if the prior knowledge is
available to all the methods, the GST is almost as fast as
0.1
the standard synchrosqueezing transform). Moreover, a
180
long signal may also slow down the GST but it does the
0.08 same (or even more so) to all most other methods.
140 In the above, we have shown that the GST could
Frequency (Hz)

0.06 produce better TFR for extracting both the IF and IA

100 trajectories. Furthermore, we wish to reiterate that the
0.04 TF resolution is always limited by the Heisenberg uncer-
tainty principle [13] regardless which post-processing
60
0.02 enhancement approaches, e.g., synchrosqueezing or reas-
signment, is adopted. However, it is important to point
20 out that the proposed GST is no longer a simple post-
0
0 200 400 600 800 1000 processing enhancement step. The operations of the GST
Time (s) include TF plane manipulation (i.e. IF trajectory mapping),
TF decomposition (i.e. CWT), and post-processing (syn-
0.1 chrosqueezing). Hence, the TF resolution dilemma may be
1200 alleviated through refining frequency resolution and time
0.08 resolution respectively in different steps. This will be
1000
illustrated in the following simulation study.
Frequency (Hz)

0.06 800
4. Simulation case study
0.04 600

400 A simulation study is presented to demonstrate the

0.02 GST’s capability in TF resolution refinement even for noisy
200 signals. We simulate the meshing frequency of a gearbox
0 modulated by a faulty gear’s rotation frequency and its
0 200 400 600 800 1000 harmonics. The gearbox fault can be reflected by the
Time (s) existence of sidebands [26]. A simulated vibration signal
Fig. 4. TFRs of s(t) obtained respectively by: (a) SPWVD, (b) reassigned
x(t) of a faulty gearbox with variable speed is given by
SPWVD, and (c) reassigned Morlet CWT. (Note: The values of the TFR
8
> xðtÞ ¼ sðtÞ þ dðtÞ
coefficients in the three subplots are not in the same scale and hence >
>
>
> sðtÞ ¼ s1 ðtÞ þ s2 ðtÞ þ s3 ðtÞ
such values from different methods are not comparable.) >
<
s1 ðtÞ ¼ cosð50pð3t þ 5t 2 2t 3 Þ þ 2:5 sinð0:6ptÞÞ ð28Þ
>
>
>
> s2 ðtÞ ¼ 0:5 cosð60pð3t þ 5t 2 2t 3 Þ þ 2:5 sinð0:6ptÞÞ
>
>
GST is for Hilbert transforms which is negligible comparing : s ðtÞ ¼ 0:5cosð40pð3t þ5t 2 2t 3 Þ þ 2:5 sinð0:6ptÞÞ
3
to the synchrosqueezing itself. The bi-axial assignment, on
the contrary, requires much more computations than the where t 2 ½0,1, s1(t), s2(t) and s3(t) respectively denote the
proposed GST technique. gear meshing component and two sideband components
However, the proposed algorithm may still not be (voltage V), and dðtÞ is a white Gaussian noise with signal-
applicable in real time situations, though it is much faster to-noise ratio (SNR) SNRðsðtÞÞ ¼ 10 lgðPðsðtÞÞ=PðdðtÞÞÞ ¼ 3 dB,
than some popular reassignment methods [18]. When where P(.) denotes the average power of a signal. It should
running on a laptop (Intel i3-370 CPU, 2G DDR3 memory, be noted that the above SNR(s(t)) is defined using the total
Matlab version 2007b), it takes more than 10 s to compute power of the signal s(t). Since the simulated signal given
the TFR for s(t) given by (14) using either the reassigned by (28) is of multi-component, we calculate hereafter
Morlet CWT or the reassigned smoothed pseudo-Wigner– different SNRs in terms of three constituent components
Ville distribution [18]. On the contrary, it only takes the s1(t), s2(t) and s3(t) respectively. In this study, the sam-
synchrosqueezing or the GST less than 1 s to generate the pling frequency is 4 kHz.
 C. Li, M. Liang / Signal Processing 92 (2012) 2264–2274 2271

It is calculated that average powers of the simulated 250t150t 2 þ 0:65 cosð0:6ptÞ that is the analytic form
signal and its constituent components are respectively of the meshing frequency acquired from the assumed
P(s(t))¼R(rms(s(t)))2 E0.7497RW, P(s1(t))E0.4981RW, ‘‘tachometer’’.
P(s2(t))¼P(s3(t))E0.1258RW, and P(dðtÞ)E0.3757RW, According to the analysis in Section 3.2, we choose
where R stands for the system impedance. One may f0 ¼50 Hz, which is a lower frequency horizontal ridge (to
accordingly calculate the signal-component-to-noise ratios improve the frequency resolution), in attempt to ensure
of s1(t), s2(t) and s3(t), i.e., SNR(s1(t))¼1.22 dB, and that all components be in the positive frequency territory
SNR(s2(t))¼SNR(s3(t))¼  4.75 dB. The time-domain wave- during the IF trajectory mapping. According to the chosen
form and the synchrosqueezed Morlet CWT of s(t) are f0 and the f(t) acquired from the assumed ‘‘tachometer’’,
respectively plotted in Fig. 5(a) and (b). Due to the higher one can calculate x0(t) using (17) as well as the mapping
energy level of s1(t) and the low resolution in the function ei2px0 ðtÞ . Hence the assumed ‘‘tachometer’’ is
frequency window, none of the two sidebands can be considered the prior knowledge of the mapping function.
observed in Fig. 5(b). Mapping s(t) to the target signal with constant fre-
Usually a tachometer is installed for a gearbox. Accord- quency f0 ¼50 Hz (Fig. 5(c)) reveals a clear lower sideband
ing to (1), (2) and (28), the meshing frequency of the and a blurred and distorted upper sideband in s(t). This
simulated gearbox signal can be specified as f ðtÞ ¼ 75 þ validates that the GST can indeed mitigate the TF

5 250 250
Observed meshing frequency

200 200
Frequency (Hz)
Amplitude

150 150
0
100 100

50 50

-5 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
Time (s) Time (s)

250 250
250 250

200 200 200

200
Frequency (Hz)

Frequency (Hz)

150 150 150 150

100 Upper sideband 100 Upper sideband

100 100

Meshing frequency
50 50 50 Meshing frequency 50
Lower sideband
Lower sideband
0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
Time (s) Time (s)

250 250

200 200
Frequency (Hz)

150 150

100 100
Segmented upper sideband

50 50

0
0 0.2 0.4 0.6 0.8 1
Time (s)

Fig. 5. TFRs of simulated signal using the proposed GST approach: (a) time-domain waveform of s(t) defined by (28), (b) synchrosqueezing-based TFR of
s(t), (c) GST-based TFR using real f1(t), (d) GST-based TFR using estimated f1(t) given by (29), and (e) GST result of s2(t).
2272 C. Li, M. Liang / Signal Processing 92 (2012) 2264–2274

resolution dilemma (i.e., it shows a higher frequency Table 1

resolution without compromising the time resolution). SNRs of different white Gaussian noises.
To further challenge this method, we assume that the
Power of dðtÞ SNR(s(t)) SNR(s1(t)) SNR(s2(t)) SNR(s3(t))
tachometer of the gearbox is not available. This means (dB) (dB) (dB) (dB)
(RW)
that we no longer have the prior knowledge of the
mapping function. Fortunately, one can still estimate the 0.7497 0  1.78  7.75  7.75
meshing frequency from the synchrosqueezing-based 2.3710 5  6.78  12.75  12.75
7.4970  10  11.78  17.75  17.75
TFR, although observation error may exist. Therefore, we
estimate the meshing frequency from Fig. 5(b) approxi-
mated by linear segments other third-party meters available). This provides prior
( knowledge of the mapping function ei2px0 ðtÞ (and the
75 þ100t=0:6 0 rt r 0:6 mapped s1(t) can be guaranteed to be true horizontal).
f 1 ðtÞ ¼ ð29Þ
175 0:6o t r1 Fig. 6(d)–(f) respectively illustrate the TFRs generated by
the proposed GST approach for the three simulated signals.
The observed meshing frequency using the linear seg- We wish to point out that a mapping function will not
ments is also shown in Fig. 5(b). The target frequency is benefit the standard synchrosqueezing method because it
chosen again as f0 ¼50 Hz (c(t)¼0). Owing to the difference does not need such information. It should also be noted that
between the real and the estimated meshing frequencies, if the meshing speed is unavailable, one can still estimate the
however, the target frequency using the estimated meshing meshing frequency from the standard synchsqueezing trans-
frequency is not really constant (horizontal in the TFR plot). form. However, if the meshing frequency trajectory is com-
The GST-based TFR using the estimated f1(t) expressed by pletely masked by the noise (e.g. SNR(s(t))¼ 10 dB) and the
(29) is shown in Fig. 5(d) that is similar to Fig. 5(c). Compar- speed information cannot be provided by a third-party
ing Fig. 5(c) and (d) suggests that informative results can still device, the GST will not be able to produce a clear TFR.
be obtained even in the presence of the mapping function Comparing with the standard synchrosqueezing-based TFR, it
error (caused by replacing a curve with two linear segments). is shown that the proposed GST-based TFR features better
From either Fig. 5(c) or Fig. 5(d) the sidebands of the gearbox noise immunity in addition to alleviating the time–frequency
can be identified and thus the gearbox fault is detected in our resolution dilemma. This more clearly demonstrates the
simulations. effectiveness of the proposed GST approach in enhancing
The ability for partial signal segmentation is another the signal TFR.
added benefit of the GST. With the mapping of ei2px0 ðtÞ , part
of the signal mixture may be transformed to the negative 5. Conclusions
frequency band if f ðtÞx0’ ðtÞ o 0 (see (17)). According to (23),
the negative frequency components are removed directly A generalized synchrosqueezing transform has been
from the TFR. In this way, any unwanted portion of a signal developed to enhance signal TFR. This method can per-
can be removed from the original signal TFR. form a wide range of TF plane transforms such as transla-
To illustrate, we apply the partial signal segmentation tion, rotation and twisting within one framework. We
feature of the GST to make the upper sideband in Fig. 5(c) have observed that, for the signal with a constant fre-
more observable. In doing so, we map the upper sideband quency, the wavelet diffusion only occurs at frequency
frequency f2(t)¼0:75cosð0:6ptÞ180t 2 þ 300t þ90 to f0 ¼ dimension. This observation has led us to introduce an IF
10 Hz. After this GST operation, both s1(t) and s3(t) are trajectory mapping operation to eliminate the time-
segmented and removed from s(t). Fig. 5(e) illustrates the dimensional spread whereas synchrosqueezing is
GST result of s2(t). employed to suppress the frequency-dimensional diffu-
In the above simulation case, we illustrated the proposed sion. With the combination of signal transform and TF
method at a relatively high SNR (3 dB). We now examine the enhancement, the GST is able to create a sharpened TFR in
capability of the proposed GST method of processing signals both time and frequency dimensions.
with lower SNRs. For this purpose, we fix the constituent The advantages of the proposed approach include: (a) it
components s1(t), s2(t) and s3(t) and increase the power of the can produce a clearer TFR without excessive computations
added white Gaussian noise dðtÞ shown in (28) to 0.7497RW, and trajectory distortion comparing to the several existing
2.3710RW and 7.4970RW respectively. The corresponding TFR tools, and (b) it helps to alleviate the TF resolution
SNRs of s(t), s1(t), s2(t) and s3(t) different components are dilemma reflected by the Heisenberg uncertainty principle.
calculated and listed in Table 1. This novel approach has been used to extract a simulated
The synchrosqueezing-based TFRs for the signals with gearbox fault signature from a noisy signal. The result shows
SNR ¼0, –5 and –10 dB are plotted in Fig. 6(a)–(c) respec- that the proposed approach is effective for enhancing the
tively. With the rise of the noise strength, it becomes signal TFR.
increasingly difficult to estimate the meshing frequency
using the standard synchrosqueezing. When SNR(s(t)) ¼
 10 dB, especially, the meshing frequency trajectory is Acknowledgments
almost completely masked by the noise as shown in
Fig. 6(c). This work is supported in part by the Natural Sciences
Suppose that the meshing speed information of the and Engineering Research Council of Canada, the Ontario
simulated gearbox signal is available (e.g. speedometer or Centers of Excellence, the National Natural Science
 C. Li, M. Liang / Signal Processing 92 (2012) 2264–2274 2273

250 250 250 250

200 200 200 200

Frequency (Hz)

Frequency (Hz)
150 150 150 150

100 100 100 100

50 50 50 50

0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
Time (s) Time (s)

250 250 250 250

200 200 200 200

Frequency (Hz)

Frequency (Hz)
150 150 150 150

100 100 100 100

50 50 50 50

0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
Time (s) Time (s)

250 250 250 250

200 200 200 200

Frequency (Hz)

Frequency (Hz)

150 150 150 150

100 100 100 100

50 50 50 50

0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
Time (s) Time (s)

Fig. 6. TFRs of simulated signals with different noise levels: (a), (b) and (c) are temporal waveforms of s(t) defined by (28) with SNR(s(t)) ¼0,  5, and
 10 dB respectively; (d)–(f) are GST-based TFRs for the three simulated signals shown in (a)–(c) respectively.

Foundation of China (50905193) and the Natural [5] W. Abdulla, L. Wong, Neonatal EEG signal characteristics using time
frequency analysis, Physica A-Statistical Mechanics and Its Applica-
Science Foundation Project of CQ CSTC (2010BB3312, tions 390 (2011) 1096–1110.
2010BB4261). [6] V.N. Ivanovic, S. Jovanovski, Signal adaptive system for time–
frequency analysis, Electronics Letters 44 (2008) 1279–1280.
[7] D. Gabor, Theory of communication, Journal of Institution of
References Electrical Engineering 93 (1946) 429–457.
[8] E.P. Wigner, On the quantum correction for thermodynamic equili-
[1] E. Sejdic, E. Djurovic, J. Jiang, Time–frequency feature representa- brium, Physical Review 40 (1946) 749–759.
tion using energy concentration: an overview of recent advances, [9] A. Lefebvre, T. Corpetti, L.H. Moy, Estimation of the orientation of
Digital Signal Processing 19 (2009) 153–183. textured patterns via wavelet analysis, Pattern Recognition Letters
[2] C. Tantibundhit, F. Pernkopf, G. Kubin, Joint time–frequency seg- 32 (2011) 190–196.
mentation algorithm for transient speech decomposition and [10] J.G. Zhong, Y. Huang, Time–frequency representation based on an
speech enhancement, IEEE Transactions on Audio Speech and adaptive short-time Fourier transform, IEEE Transactions on Signal
Language Processing 18 (2010) 1417–1428. Processing 58 (2010) 5118–5128.
[3] L.R. Padovese, Hybrid time–frequency methods for non-stationary [11] O. Talakoub, J. Cui, W. Wong, Approximating the time-frequency
mechanical signal analysis, Mechanical Systems and Signal Proces- representation of biosignals with chirplets, EURASIP Journal on
sing 18 (2004) 1047–1064. Advances in Signal Processing 2010 (2010) 857685.
[4] C. Li, M. Liang, Separation of the vibration-induced signal of oil [12] L. Zhang, G.L. Xiong, H.S. Liu, W.Z. Guo, Time–frequency representation
debris for vibration monitoring, Smart Structures & Structures 20 based on time-varying autoregressive model with applications to non-
(2011) 045016. stationary rotor vibration analysis, Sadhana 35 (2010) 215–232.
2274 C. Li, M. Liang / Signal Processing 92 (2012) 2264–2274

[13] S.G. Mallat, A Wavelet Tour of Signal Processing, Academic, [20] I. Daubechies, Ten Lectures on Wavelets, CBMS-NSF Regional
New York, 1998. Conference Series in Applied Mathematics, 61, SIAM, Philadelphia,
[14] Y. Dai, Q. Ma, W.Y. Tang, Efficient wavelet ridge extraction method PA, 1992.
for asymptotic signal analysis, Review of Scientific Instruments 79 [21] J.M. Lilly, J.C. Gascard, Wavelet ridge diagnosis of time-varying
(2008) 124703. elliptical signals with application to an oceanic eddy, Nonlinear
[15] I. Daubechies, J.F. Lu, H.T. Wu, Synchrosqueezed wavelet trans- Processes in Geophysics 13 (2006) 467–483.
forms: an empirical mode decomposition-like tool, Applied and [22] I. Djurovic, Viterbi algorithm for chirp-rate and instantaneous
Computational Harmonic Analysis 30 (2011) 243–261. frequency estimation, Signal Processing 91 (2011) 1308–1314.
[16] G. Thakur, H.T. Wu, Synchrosqueezing-based recovery of instanta- [23] J.S. Cheng, Y. Yang, D.J. Yu, The envelope order spectrum based on
neous frequency from nonuniform samples, SIAM Journal on generalized demodulation time-frequency analysis and its applica-
Mathematical Analysis (2011). tion to gear fault diagnosis, Mechanical Systems and Signal Proces-
[17] E. Brevdo, H.T. Wu, G. Thakur, N.S. Fukar, Synchrosqueezing and its sing 24 (2010) 508–521.
applications in the analysis of signals with time-varying spectrum, [24] S. Olhede, A.T. Walden, A generalized demodulation approach to
Proceedings of the National Academy of Sciences of the United time–frequency projections for multicomponent signals, Proceed-
States of America (2011). ings of the Royal Society A-Mathematical Physical and Engineering
[18] F. Auger, P. Flandrin, Improving the readability of time-frequency Sciences 461 (2005) 2159–2179.
and time-scale representations by the reassignment method, IEEE [25] P.L. Shui, H.Y. Shang, Y.B. Zhao, Instantaneous frequency estimation
Transaction on Signal Processing 43 (1995) 1068–1089. based on directionally smoothed pseudo-Wigner–Ville distribution
[19] I. Daubechies, S. Maes, A nonlinear squeezing of the continuous bank, IET Radar Sonar and Navigation 1 (2007) 317–325.
wavelet transform based on auditory nerve models, in: A. Aldroubi, [26] W.D. Mark, Stationary transducer response to planetary-gear vibra-
M. Unser (Eds.), Wavelets in Medicine and Biology, CRC Press, 1996, tion excitation II: effects of torque modulations, Mechanical Systems
pp. 527–546. and Signal Processing 23 (2009) 2253–2259.