Time–frequency chirp-Wigner transform for signals with any nonlinear polynomial time varying instantaneous frequency

ARTICLE IN PRESS Mechanical Systems and Signal Processing Mechanical Systems and Signal Processing 21 (2007) 2980–3002 www.elsevier.com/locate/jnlabr/ymssp

Time–frequency chirp-Wigner transform for signals with any nonlinear polynomial time varying instantaneous frequency L. Gelman, J.D. Gould Department of Process and Systems Engineering, School of Engineering, Cranfield University, MK43 0AL, UK Received 13 November 2006; received in revised form 10 May 2007; accepted 10 May 2007 Available online 25 May 2007

Abstract The new technique, the time–frequency chirp-Wigner transform has been proposed recently. This technique is further investigated for the general case of higher order chirps, i.e. non-stationary signals with any nonlinear polynomial variation of the instantaneous frequency in time. Analytical and numerical comparison of the chirp-Wigner transform and the classical Wigner distribution was performed for processing of single-component and multi-component higher order chirps. It is shown for the general case of single component higher order chirps that the chirp-Wigner transform has an essential advantage in comparison with the traditional Wigner distribution: the chirp-Wigner transform ideally follows the nonlinear polynomial frequency variation without amplitude errors. It is shown for multi-component signal where each component is a higher order chirp, that the chirp-Wigner transform adjusted to a single component will follow the instantaneous frequency of the component without amplitude errors. It is also shown that the classical Wigner distribution is unable to estimate component amplitudes of single component and multi-component higher order chirps. r 2007 Elsevier Ltd. All rights reserved. Keywords: Wigner distribution; Time–frequency analysis; Nonlinear variation of the instantaneous frequency; Amplitude estimation

1. Introduction The Wigner distribution has been widely used [1–11] for processing of chirps, i.e. transient signals with constant frequency speeds (chirp rates). These signals have linear variation of the instantaneous frequency in time. The Wigner distribution is suitable for the chirps; the distribution ideally follows the linear variation of the instantaneous frequency without amplitude errors. However, in practical applications, for radar, sonar and mechanical systems, e.g. change of the instantaneous frequency of radar and sonar signals, start-up or shut-down of rotating machinery or change of the shaft instantaneous frequency during machinery operation, one has to process transient signals with any nonlinear polynomial variation of the instantaneous frequency in time. These signals could characterize radar, sonar signals and the shaft and blade passing vibrations of rotating machinery during nonlinear change of shaft frequency. Corresponding author. Tel.: +44 1234750111x5425; fax: +44 1234 754797.

E-mail address: l.gelman@cranfield.ac.uk (L. Gelman). 0888-3270/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2007.05.003

ARTICLE IN PRESS L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002

2981

The Wigner distribution is not suitable for these signals, because the kernel of the distribution is not suitable for nonlinear polynomial frequency variation. The polynomial Wigner–Ville distribution [12–15] has been proposed and investigated for the higher order chirps, i.e. signals with the nonlinear polynomial variation of the instantaneous frequency in time. This very useful transform is a generalization of the Wigner distribution. The disadvantage of the polynomial Wigner–Ville distribution is that the multi-linear kernel of the distribution creates a multiplicity of the cross-terms for multi-component signals. For example [14], if one processes only a two-component signal, the number of the cross-terms is 13 in the polynomial Wigner–Ville distribution of order 4. Application of this distribution requires knowledge (a priori or from independent measurement) of the order of the polynomial variation of the instantaneous frequency. The higher order local polynomial Wigner distribution has been proposed and investigated by Katkovnik [16,17] and investigated by Stankovic [18] for estimation of the unknown variable instantaneous frequency and its derivatives for the higher order chirps. This important and very useful transform is a higher order generalization of the Wigner distribution. It is shown [16] that in the general case the estimation error provided by the transform is not equal to zero. Quite often, for transient vibro-acoustical signals from complex mechanical systems (e.g. turbines, gearboxes, etc.) time variations of the instantaneous frequency and its derivatives are known from independent synchronous measurements (e.g. from tachometer signals). Sometimes, these frequency variations even are known a priori. Therefore, the main problem for these mechanical systems is amplitude estimation for transient signals (on the background of interference) with known nonlinear time variation of the instantaneous frequency. A new time–frequency technique, the chirp-Wigner transform has been proposed and investigated in Ref. [19] for amplitude estimation of the higher order chirps (on the background of interference) with the known nonlinear polynomial time variation of the instantaneous frequency and its derivatives. The transform is compared [19] with the traditional Wigner distribution. It is shown [19] that the main advantage of the proposed transform is that the transform could ideally follow the nonlinear polynomial frequency variation without amplitude errors. However, in spite of the fact that the chirp-Wigner transform is suitable for any nonlinear polynomial frequency variation, the transform was investigated and compared [19] with the traditional Wigner transform only for the particular case of order four of the nonlinear polynomial variation of the instantaneous frequency. Therefore, the purposes of this paper are:

  

To investigate further the chirp-Wigner transform for the general case of known arbitrary nonlinear polynomial variation of the instantaneous frequency To derive a general formula for the time-varying coefficients of the chirp-Wigner transform for transform adjustment To compare the chirp-Wigner transform with the traditional Wigner distribution for the general case of known arbitrary nonlinear polynomial variation of the instantaneous frequency for single-component and multi-component signals.

The topic is not well known within the mechanical engineering community and hence should be brought to its attention. 2. The chirp-Wigner transform Consider a general complex non-stationary signal with constant amplitude A and any known nonlinear polynomial variation of the instantaneous frequency in time [19]:   xðtÞ ¼ A exp j2p½f 0 t þ V ðtÞ , (1) where f0 is the initial frequency, t is time, V ðtÞ ¼ ðc20 t2 =2Þ þ ðc30 t3 =3Þ þ ðc40 t4 =4Þ þ    þ ðcn0 tn =nÞ, c20 is the constant chirp rate (i.e. the frequency speed) of the signal, c30 is the constant frequency acceleration of the signal, c40 ; . . . ; cn0 are higher order constant parameters which determine the nonlinear variation of the

ARTICLE IN PRESS 2982

L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002

instantaneous frequency in time; n is an order of linear or nonlinear frequency variation, nX2, n ¼ 2 is the linear case, nX3 is the nonlinear case. The nonlinear variation of the instantaneous frequency can be written as follows: f i ðtÞ ¼ f 0 þ c20 t þ c30 t2 þ c40 t3 þ    þ cn0 tn1 . This signal is called a higher order chirp [19]. The chirp-Wigner transform is as follows [19]: Z 1  t  t 2 3 4 n W c ðt; f Þ ¼ xn t  x t þ ej2p½f tþðc2 ðtÞ=2Þt þðc3 ðtÞ=3Þt þðc4 ðtÞ=4Þt þþðcn ðtÞ=nÞt  dt, 2 2 1

(2)

(3)

where f is the frequency; t is the time delay, c2 ðtÞ and c3 ðtÞ are the time-varying chirp rate (i.e. frequency speed) and frequency acceleration respectively of the transform, c4 ðtÞ; c5 ðtÞ; . . . ; cn ðtÞ are the time-varying higher order parameters of the transform, * is a symbol for complex conjugation. The main difference between the chirp-Wigner transform and the higher order local polynomial Wigner distribution by Katkovnik [16–18] is that Katkovnik’s transform has been developed for estimation of unknown instantaneous frequency variation and its derivatives, whereas the chirp-Wigner transform is developed for estimation of amplitude of signal (on the background of interference) with known instantaneous frequency variation and its derivatives. The chirp rate, the frequency acceleration and the higher order parameters of the chirp-Wigner transform are time-varying, whereas the appropriate parameters of the higher order local polynomial Wigner distribution are time-invariant. Time variation of the transform parameters provides an essential advantage of the chirp-Wigner transform. The chirp-Wigner transform is bilinear, therefore it has a limited number of cross-terms; this property is similar to the Wigner and Katkovnik’s transforms and in contrast to the polynomial Wigner distribution. The chirp-Wigner transform is a generalization of the higher order local polynomial Wigner distribution [16–18] for the case of the known nonlinear polynomial variation of the instantaneous frequency and its derivatives. 3. Adjustment of the chirp-Wigner transform for any nonlinear polynomial variation of the instantaneous frequency Adjustment of the chirp-Wigner transform is proposed and investigated in Ref. [19]. The purpose of adjustment is that chirp-Wigner transform with adjustment should ideally follow any nonlinear polynomial variation of the instantaneous frequency without amplitude errors. However, in spite of the fact that the chirpWigner transform is suitable for any nonlinear polynomial frequency variation, the transform adjustment was investigated [19] only for the particular case of order four of the nonlinear polynomial variation of the instantaneous frequency. Let us investigate the adjustment [19] of the chirp-Wigner transform for the general case of any order of nonlinear polynomial variation of the instantaneous frequency. In order to adjust the chirp-Wigner transform, one needs:

 



to find relationships between signal parameters (e.g., frequency speed, frequency acceleration, etc.) and time-varying transform parameters (e.g. frequency speed, frequency acceleration, etc.) which provide transform adjustment for any order of nonlinear polynomial variation of the instantaneous frequency to evaluate continuously from known or measured variation of the instantaneous frequency the order n, frequency speed c20, frequency acceleration c30 and the higher order parameters c40 ; . . . ; cn0 of the frequency variation; this estimation could be easily done using standard polynomial approximation techniques and known higher order derivatives of the frequency variation to process the signal by the chirp-Wigner transform with transform parameters defined from the relationships found.

Therefore, let us find firstly the relationships between signal parameters and time-varying transform parameters which provide transform adjustment. Calculating the chirp-Wigner transform W c ðt; f Þ of the

ARTICLE IN PRESS L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002

higher order chirp of order n in Eq. (1) gives the following: ( " #) Z 1 r     X t t c ðtÞ 2kþ1 t2kþ1 W c ðt; f Þ ¼ A2 exp j2p y t þ y t  ft dt, 2 2 2k þ 1 1 k¼1

2983

(4)

where 2pyðtÞ ¼ 2pðf 0 t þ V ðtÞÞ is the instantaneous phase of the signal; 2pfyðt þ ðt=2ÞÞ  yðt  ðt=2ÞÞg is the instantaneous phase of the instantaneous autocorrelation function of the signal. It is shown in Appendix A that the phase of the instantaneous auto-correlation function appearing in Eq. (4) satisfies: r   X t t V ð2kþ1Þ ðtÞ 2kþ1 y t ¼ f i ðtÞt þ y tþ t . (5) 2k 2 2 k¼1 2 ð2k þ 1Þ! The following relationships were found between signal parameters and time-dependent transform parameters which provide transform adjustment: c2k ðtÞ ¼ 0, c2kþ1 ðtÞ ¼ V ð2kþ1Þ ðtÞ=22k ð2kÞ!;

k ¼ 1; . . . ; r,

ð6Þ

ðkÞ

where V ðtÞ denotes the derivative of order k of V(t), r denotes the largest integer less than or equal to (n1)/2 (so that 2r+1 is the largest odd integer less than or equal to n) i.e. r ¼ [(n1)/2], the square brackets denote the integer part of (n1)/2. If we substitute the expression (6) into expression (4), and make use of the relationship (5), we can write the transform as Z 1   W c ðt; f Þ ¼ A2 expfj2p½ðf i ðtÞ  f Þtg dt ¼ A2 d f  f i ðtÞ , (7) 1

where dðf Þ is the Dirac function. The coefficients have been defined in Eq. (6) to ensure precisely that the chirp-Wigner transform evaluates to a Dirac function in Eq. (7). One can see from Eq. (7) that the chirp-Wigner transform with adjustment conditions (6) satisfied is ideally concentrated along the nonlinear polynomial frequency variation of the signal under consideration without amplitude errors. Suppose that the order of the higher-order chirp in Eq. (2) is four, i.e. n ¼ 4. Then, it follows from Eq. (6) that c2 ðtÞ ¼ c4 ðtÞ ¼ 0, and c3 ðtÞ ¼ V ð3Þ ðtÞ=8 ¼ ð2c30 þ 6c40 tÞ=8 ¼ ðc30 þ 3c40 tÞ=4 agreeing with the case given in Ref. [19]. Note that the even coefficients of the chirp-Wigner transform are zero, and so the transform requires about half the number of coefficients that are needed to specify the higher order chirp (1). The fact that the even coefficients are zero also implies that its chirp-Wigner transform is real valued. To see this, note that the chirp-Wigner transform is the Fourier transform of the following function:  t  t 3 5 2rþ1 x t  x t þ ej2p½ðc3 ðtÞ=3Þt þðc5 ðtÞ=5Þt þþðc2rþ1 ðtÞ=2rþ1Þt  , 2 2 which has Hermitian symmetry in t. Using Eq. (5) one can calculate the Wigner distribution of the higher-order chirp of order n to get ( " #) Z 1 r X V ð2kþ1Þ ðtÞ 2kþ1 2 W ðt; f Þ ¼ A t expfj2pðf i ðtÞ  f Þtg exp j2p dt. 2k 1 k¼1 2 ð2k þ 1Þ!

(8)

One can see from Eq. (8), that because of the presence of the second exponential factor in the integrand, the integral will not in general evaluate to the Dirac function expression A2 d½f  f i ðtÞ as in Eq. (7). Therefore, because of this non-unity time-varying multiplier, the traditional Wigner distribution (8) does not ideally follow the nonlinear frequency variation of the higher order chirps. But, note that it is not obvious from Eq. (8) the extent to which the Wigner distribution will differ in its behavior from the chirp-Wigner distribution. It is one of the aims of this paper to investigate this difference.

ARTICLE IN PRESS L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002

2984

4. Theoretical analysis In view of the conditions given in Eq. (6), it makes sense to rewrite the definition of the chirp-Wigner transform given in Eq. (3) by using time-dependent coefficients of the same form as in those appearing in the conditions (6). Thus, the chirp-Wigner transform will be given by the expression: Z 1  t  t W c ðt; f ; l0 ðtÞÞ ¼ xn t  x t þ K c ðt; l0 ðtÞÞej2pf t dt, (9) 2 2 1 where the transform kernel is of the form:

  c3 ðt; l0 Þ 3 c5 ðt; l0 Þ 5 c2rþ1 ðt; l0 Þ 2rþ1 K c ðt; l0 ðtÞÞ ¼ exp j2p t þ t þ  þ t . 3 5 2r þ 1

(10)

In Eq. (9), the dependence of the kernel on a vector of time-dependent transform parameters l0 ðtÞ ¼ ðc3 ðt; l0 Þ; c5 ðt; l0 Þ; . . . ; c2rþ1 ðt; l0 ÞÞ is made explicit. The vector has components c2kþ1 ðt; l0 Þ of the form: c2kþ1 ðt; l0 Þ ¼ V ð2kþ1Þ ðt; l0 Þ=22k ð2kÞ!, ðkÞ

0

(11) 0

0

ðf 00 ; c020 ; c030 ; . . . ; c0n0 Þ

where V ðt; l Þ denotes the derivative of order k of V ðt; l Þ, and l ¼ is a constant vector equivalent via Eq. (11) to l0 (t). Thus, conditions (6) are incorporated into the transform kernel (10), by assuming that the coefficients c2kþ1 ðt; l0 Þ have the same functional form as in conditions (6), rather than being arbitrary functions, and by assuming that the coefficients c2k ðt; l0 Þ are identically zero. In this case, the conditions (6) can be written in time-dependent form as c2kþ1 ðt; l0 Þ ¼ c2kþ1 ðt; lÞ, or in terms of constant form as c0j0 ¼ cj0 ;

j ¼ 2; . . . ; n,

  where the constant vector l ¼ f 0 ; c20 ; c30 ; . . . ; cn0 specifies the signal (1). For polynomial variation of instantaneous phase of degree 6 the customized chirp-Wigner transform is given by  0  Z T  t  t c30 þ 3c040 t þ 6c050 t2 þ 10c060 t3 3 xn t  x t þ expðj2p W c ðt; f ; l0 ðtÞÞ ¼ t 2 2 12 T

0   c þ 5c060 t 5 þ 50 t ej2pf t dt, 80 and the adjustment conditions are c0j0 ¼ cj0 , j ¼ 2; . . . ; 6. There is no analytic expression for the Wigner distribution of a signal with arbitrary polynomial variation of the instantaneous phase. However, it is possible to derive formulae in certain special cases. The cases of general polynomial variation of the instantaneous phase of degrees n ¼ 3, 4, 5, 6 are derived in Appendix A and B and the formulae are given here. Firstly, for polynomial variation of the instantaneous phase of degree 4 of a finite-duration signal for which the duration is in the interval [T, T], the following formula holds for the Wigner distribution: 8 1    n 2 P ð2pc3 ðt;lÞ=3Þ > > Im G 3n þ 1; j2pðf i ðtÞ  f ÞT ; f af i ðtÞ; > 2A ð2pðf i ðtÞf ÞÞ3nþ1 n! < n¼0 (12) W ðt; f Þ ¼ 1 P > ð1Þn ð2pc3 ðt;lÞ=3Þ2n T 6nþ1 2 > ; f ¼ f ðtÞ: 2A > i ð2nÞ! ð6nþ1Þ : n¼0

where

c3 ðt; lÞ ¼

 c30 þ 3c40 t , 4

ARTICLE IN PRESS L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002

2985

R1 Gða; zÞ ¼ z ta1 expðtÞ dt is the incomplete gamma function, l ¼ ðf 0 ; c20 ; c30 ; c40 Þ is the constant parameter vector of the signal. If the instantaneous phase is of degree 3, then c40 is zero. For each time t, the expression for W(t, f) is continuous as a function of frequency f, at all frequencies including the instantaneous frequency of the signal fi(t). The shape of the surface of the Wigner distribution W(t, f) is not immediately apparent from formula (12). For polynomial variation of the instantaneous phase of degree 6 of a finite-duration signal for which the duration is in the interval [T, T], the following expression holds for the Wigner distribution: 8

 1 n n ð1Þk ð2pc5 ðt;lÞ=5Þk ð2pc3 ðt;lÞ=3Þnk > 2 P 1 P > > ImðGð3n þ 2k þ 1; j2pðf i ðtÞ  f ÞTÞÞ; f af i ðtÞ; > ð2pðf i ðtÞf ÞÞ3nþ2kþ1 < 2A n¼0 n! k¼0 k

 W ðt; f Þ ¼ 1 2n 2n > k 2nk T 6nþ2kþ1 2 P ð1Þn P > > f ¼ f i ðtÞ: > ð6nþ2kþ1Þ ; : 2A n¼0 ð2nÞ! k¼0 k ð2pc5 ðt; lÞ=5Þ ð2pc3 ðt; lÞ=3Þ (13) where

 c30 þ 3c40 t þ 6c50 t2 þ 10c60 t3 c3 ðt; lÞ ¼ ; 4

c5 ðt; lÞ ¼

 c50 þ 5c60 t ; 16

  l ¼ f 0 ; c20 ; c30 ; c40 ; c50 ; c60

is the constant parameter vector of the signal. If the phase is of degree 5, then c60 is zero. There is no analytic expression for the chirp-Wigner distribution of a signal with arbitrary polynomial variation of the instantaneous phase. However, as for the Wigner distribution, it is possible to derive formulae in certain special cases. For polynomial variation of the instantaneous phase of degree 4 and a finite-duration signal for which the duration is in the interval [T, T], the following formula holds for the chirp-Wigner transform (see Appendix B): 8 1 ð2pDc3 =3Þn 2 P > > ImðGð3n þ 1; j2pðf i ðtÞ  f ÞTÞÞ; f af i ðtÞ; > ð2pðf ðtÞf ÞÞ3nþ1 n! < 2A W c ðt; f ; l0 ðtÞÞ ¼

> 2 > > : 2A

n¼0 1 P n¼0

i

n

ð1Þ ð2pDc3 =3Þ2n T 6nþ1 ð2nÞ! ð6nþ1Þ ;

(14)

f ¼ f i ðtÞ:

where Dc3 ¼ c3 ðt; lÞ  c3 ðt; l0 Þ ¼

 0  c30 þ 3c40 t c þ 3c040 t  30 ; 4 4

l ¼ ðf 0 ; c20 ; c30 ; c40 Þ

is the constant parameter vector of the signal, l0 ¼ ðf 00 ; c020 ; c030 ; c040 Þ is a constant vector equivalent via Eq. (11) to the time-dependent vector l0 ðtÞ specifying the transform. Notice that expression (14) has the same general form as expression (12). For polynomial variation of the instantaneous phase of degree 6 of a signal of finite duration for which the duration is in the interval [T, T], the following formula holds for the chirp-Wigner transform (see Appendix B): 8

 1 n n ð1Þk ð2pDc5 =5Þk ð2pDc3 =3Þnk > 2 P 1 P > > 2A ImðGð3n þ 2k þ 1; j2pðf i ðtÞ  f ÞTÞÞ; f af i ðtÞ; > n! ð2pðf i ðtÞf ÞÞ3nþ2kþ1 < k n¼0 k¼0 0

 W c ðt; f ; l ðtÞÞ ¼ 1 2n 2n > 2 P ð1Þn P > T 6nþ2kþ1 > ð2pDc5 =5Þk ð2pDc3 =3Þ2nk ð6nþ2kþ1Þ ; f ¼ f i ðtÞ: 2A > ð2nÞ! : k n¼0 k¼0 (15) where Dc3 ¼ c3 ðt; lÞ  c3 ðt; l0 Þ ¼

 0  c30 þ 3c40 t þ 6c50 t2 þ 10c60 t3 c þ 3c040 t þ 6c050 t2 þ 10c060 t3  30 , 4 4

ARTICLE IN PRESS 2986

L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002 0

Dc5 ¼ c5 ðt; lÞ  c5 ðt; l Þ ¼



 0  c50 þ 5c60 t c50 þ 5c060 t  ; 16 16

l ¼ ðf 0 ; c20 ; c30 ; c40 ; c50 ; c60 Þ

is the constant parameter vector of the signal, l0 ¼ ðf 00 ; c020 ; c030 ; . . . ; c060 Þ is a constant vector equivalent via Eq. (11) to the time-dependent vector l0 ðtÞ specifying the transform. Notice that expression (15) has the same general form as expression (13). It is tempting to try to obtain formulae corresponding to Eqs. (12)–(15) for a signal of infinite duration by taking the limit as T-N, in Eqs. (12)–(15). However, this is not permissible mathematically, and will give wrong answers. Therefore, a different approach is required which involves transforming the integral by contour integration (see Appendix C). For polynomial variation of the instantaneous phase of degree 4 and an infinite-duration signal, the following formula holds for the Wigner distribution:

 1 2A2 X Gððn þ 1Þ=3Þ ð2pðf i ðtÞ  f ÞÞn 2ðn þ 1Þp W ðt; f Þ ¼ sin , (16) n! 3 3 n¼0 ð2pc3 ðt; lÞ=3Þðnþ1Þ=3 where c3 ðt; lÞ is the same as in Eq. (12). A special case of formula (16) occurs when 2pc3 ðt; lÞ ¼ 1, for then, the right-hand side equals pAið2pðf i ðtÞ  f ÞÞ where Ai(z) is Airy’s function. The integral defining the Wigner distribution in this case is, up to a constant, Airy’s integral. Notice that formula (16) is a power series in the variable 2pðf i ðtÞ  f Þ=ð2pc3 ðt; lÞ=3Þ1=3 , so that one does not need to consider the case f ¼ fi(t) separately, as we did in the finite duration case. Note also that Eq. (16) is not equal to the limit of Eq. (12). For polynomial variation of the instantaneous phase of degree 6 of a signal of infinite-duration, the following formula holds for the Wigner distribution:



 1 n  n ð2pðf i ðtÞ  f ÞÞk ð2pc3 ðt; lÞ=3Þnk 2A2 X 1X 3n  2k þ 1 ð8n  2k þ 1Þp W ðt; f Þ ¼ G cos , 5 10 5 n¼0 n! k¼0 k ð2pc5 ðt; lÞ=5Þð3n2kþ1Þ=5 (17) where c3 ðt; lÞ, and c5 ðt; lÞ are the same as in Eq. (13). Note also that Eq. (17) is not equal to the limit of Eq. (13). For polynomial variation of the instantaneous phase of degree 4 and an infinite-duration signal, the following formula holds for the chirp-Wigner transform (see Appendix C) when Dc36¼0:

 1 2A2 X Gððn þ 1Þ=3Þ ð2pðf i ðtÞ  f ÞÞn 2ðn þ 1Þp 0 W c ðt; f ; l ðtÞÞ ¼ sin , (18) n! 3 3 n¼0 ð2pDc3 =3Þðnþ1Þ=3 where Dc3 ¼ c3 ðt; lÞ  c3 ðt; l0 Þ is the same as in formula (14). When Dc3 ¼ 0 (i.e. the adjustment conditions are satisfied), we already know that the chirp-Wigner transform is ideally concentrated along the nonlinear polynomial frequency variation of the signal under consideration without amplitude errors. For polynomial variation of the instantaneous phase of degree 6 and a signal of infinite duration, the following formula holds for the chirp-Wigner transform (see Appendix C) provided Dc56¼0:



 1 n  n ð2pðf i ðtÞ  f ÞÞk ð2pDc3 =3Þnk 2A2 X 1X 3n  2k þ 1 ð8n  2k þ 1Þp G W c ðt; f ; l0 ðtÞÞ ¼ cos , 5 10 5 n¼0 n! k¼0 k ð2pDc5 =5Þð3n2kþ1Þ=5 (19) where Dc3 ¼ c3 ðt; lÞ  c3 ðt; l0 Þ, Dc5 ¼ c5 ðt; lÞ  c5 ðt; l0 Þ are the same as in formula (15). When Dc5 ¼ 0 (i.e. the adjustment conditions are partially satisfied), one uses formula (18) for the n ¼ 3, 4 cases, with Dc3 ¼ c3 ðt; mÞ  c3 ðt; l0 Þ given by the definition for the n ¼ 5, 6 cases. When both Dc5 ¼ 0, and Dc3 ¼ 0 (i.e. the adjustment conditions are satisfied) it is already known that the chirp-Wigner transform is ideally concentrated along the nonlinear polynomial frequency variation of the signal under consideration without amplitude errors.

ARTICLE IN PRESS L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002

2987

5. Case study for single-component signal Consider the non-stationary shaft passing vibrations xðtÞ ¼ A expðj2pyðtÞÞ of rotating machinery during machinery start up with known nonlinear change of shaft speed; for this signal the instantaneous phase has the polynomial variation in time of degree 6 given by yðtÞ ¼ 2pt6 =4. Here, the parameter vector is l ¼ ðf 0 ; c20 ; c30 ; c40 ; c50 ; c60 Þ ¼ ð0; 0; 0; 0; 0; 0:25Þ. Graphs of the real part of the signal and of its instantaneous frequency over the time interval 1 sptp2 s, are shown in Fig. 1a and b, respectively. Fig. 2 displays the Wigner distribution for the signal over the time interval 1 sptp2 s, with T ¼ 1 s. Notice in Fig. 2, the clustering of non-zero values near the curve of instantaneous frequency in the time–frequency plane. We can see clearly the oscillation of the transform normal to this curve. It is possible to make out the curve of instantaneous frequency in the picture after looking at Fig. 1, but, as in the exponential case, it is not immediately obvious exactly where it is. It is not clear from Fig. 2, how the maximum magnitude of the Wigner distribution varies over time. To find out, it is necessary to look at time cross-sections of this surface. Fig. 3 shows the graphs of two cross sections of the Wigner distribution of the signal (with amplitude one) at times t ¼ 1.5 s (left), and t ¼ 2 s (right). It can be seen that the maximum magnitude on the left is smaller than maximum magnitude on the right, so the maximum magnitude varies with time. It follows that the Wigner

Fig. 1. Signal with phase 2pt6/4.

Fig. 2. The Wigner distribution of a signal with phase 2pt6/4.

ARTICLE IN PRESS 2988

L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002

Fig. 3. The Wigner distribution with phase above at t ¼ 1.5 s (solid), and t ¼ 2 s (dashed).

Fig. 4. The Wigner distribution (solid) and chirp-Wigner transform (dashed) at t ¼ 2 s.

distribution is not suitable for estimating the amplitude of a signal with polynomial variation of the instantaneous frequency. The study of the Wigner distribution in the polynomial case is continued in Fig. 4, which illustrates by example, how well the distribution tracks the instantaneous frequency of the signal. The solid curve is a crosssection of the Wigner distribution at time t ¼ 2 s, at which time the instantaneous frequency of the signal is fi ¼ 84 Hz. The dotted curve is a cross-section of the chirp-Wigner transform of the same signal with T ¼ 1 s. The peak of the dotted curve lies exactly above the point on the t-axis corresponding to the instantaneous frequency. It can be seen that the Wigner distribution does not track the instantaneous frequency well: using the maximum magnitude as an estimate will give an error of a few percent. The chirp-Wigner transform will be completely adjusted when the kernel parameter vector satisfies the condition l0 ¼ ð0; 0; 0; 0; 0; 1:5Þ, and by Eqs. (9) and (10), this is equivalent to the conditions: c3 ðtÞ ¼ ð15t3 =4Þ, c5 ðtÞ ¼ ð15t=32Þ. Fig. 5 shows the result of using the chirp-Wigner transform with exactly adjusted kernel over the time interval 1 sptp2 s, with T ¼ 1 s. The picture confirms what is already known from the study of the general case, i.e. that the new transform is non-zero on the curve of instantaneous frequency, but it is also close to zero elsewhere. Moreover, the transform is constant with value 2A2T ¼ 2, along this curve. Thus, the transform enables an estimate of the signal amplitude to be made. Also, as was mentioned above, the transform will follow the curve of instantaneous frequency ever more closely as T increases. We continue to investigate the new transform with the signal above to see how the new transform behaves with respect to different degrees of adjustment. We see from formula (11) for the chirp-Wigner transform for a degree 6 polynomial, that there are four parameters c030 , c040 , c050 , c060 , used for adjusting the kernel to the signal. The case of no adjustment corresponds to c030 ¼ c040 ¼ c050 ¼ c060 ¼ 0, and is the Wigner distribution. When adjustment is complete c030 ¼ c040 ¼ c050 ¼ 0, and c060 ¼ 1:5. Partial adjustment corresponds to other values of the parameters.

ARTICLE IN PRESS L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002

2989

Fig. 5. The chirp-Wigner transform with complete adjustment of a signal with phase 2pt6/4.

Fig. 6. The chirp-Wigner at t ¼ 2 s, f ¼ 48 Hz, against the adjustment parameter d with T ¼ 1 s.

Define the adjustment parameter d by d ¼ c030  c30 ¼ c040  c40 ¼ c050  c50 ¼ c060  c60 . Here, c30, c40, c50, c60 are the true values of the parameters for the signal, and c030 , c040 , c050 , c060 are the approximate values used in the kernel of the chirp-Wigner transform. Since we are not assuming that the parameters c30, c40, c50, c60 are non-zero, we cannot define an adjustment ratio in this case. Note that the value d ¼ 0 corresponds to complete adjustment of the transform. At t ¼ 2 s, the instantaneous frequency of the signal is fi ¼ 48 Hz. In Fig. 6, we plot the variation of the chirpWigner transform at t ¼ 2 s, f ¼ 48 Hz, with T ¼ 1 s, against the adjustment parameter d over the interval of values 2pdp2. It can be seen that the maximum value of the chirp-Wigner transform (with value 2) occurs when d ¼ 0, i.e. when the parameter values of the transform are exactly equal to those of the signal. We can also plot the variation of the chirp-Wigner transform at t ¼ 2 s, with T ¼ 1 s, against the adjustment parameter d, for values of frequency f not on the instantaneous frequency curve. Fig. 7 shows the chirp-Wigner transform at t ¼ 2 s, f ¼ 60 Hz (the instantaneous frequency of the signal is fi ¼ 60 Hz at t ¼ 2 s) over the interval of values 5pdp5. It can be seen that the value of the chirp-Wigner transform is close to zero when d ¼ 0, i.e. when the parameter values of the transform are exactly equal to those of the signal.

ARTICLE IN PRESS 2990

L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002

Fig. 7. The chirp-Wigner at t ¼ 2 s, f ¼ 60 Hz, against the adjustment parameter d with T ¼ 1 s.

Fig. 8. The chirp-Wigner transform with d ¼ 0.5.

Fig. 8 shows the chirp-Wigner transform for the same signal (with T ¼ 1 s) with adjustment parameter d ¼ 0.5. Fig. 9 shows the chirp-Wigner transform for the same signal (with T ¼ 1 s) with adjustment parameter d ¼ 0.9. Fig. 10 shows the chirp-Wigner transform for the same signal (with T ¼ 1 s) with adjustment parameter d ¼ 0.9. Fig. 5 shows the case of complete adjustment d ¼ 1, and Fig. 2 shows the Wigner distribution of the signal.

6. Case study for multi-component signal Suppose that the instantaneous phases of the separate complex components of a transient multi-component signal are known a priori or from independent synchronous measurements. Then, in this section it is shown that the chirp-Wigner transform can also be used as an effective amplitude estimator of the different components of the signal. To estimate amplitude of a particular component, take the chirp-Wigner transform of the signal adjusted to that component. The value of the transform along the known instantaneous frequency of the component is related in a simple way to the amplitude. Let’s consider without loss of generality the following non-stationary two-component signal: xðtÞ ¼ x1 ðtÞ þ x2 ðtÞ,

(20)

ARTICLE IN PRESS L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002

2991

Fig. 9. The chirp-Wigner transform with d ¼ 0.9.

Fig. 10. The chirp-Wigner transform with d ¼ 0.9.

where each component xi(t), i ¼ 1, 2, is a constant amplitude complex exponential signal given by xi ðtÞ ¼ Ai expðj2pyi ðtÞÞ. Assume that the instantaneous phase 2pyi ðtÞ of each component has nonlinear polynomial variation in time. Let W i ðt; f Þ½x denote the chirp-Wigner transform of the signal x(t), adjusted to signal component i(i ¼ 1,2). For polynomials Z,x and signal component i (i ¼ 1,2), define Ii(Z,x) by ( " #) Z T r   X t t ci;2kþ1 ðtÞ 2kþ1 t 2 cos 2p Z t þ I i ðZ; xÞ ¼ x t ft dt, (21) 2 2 2k þ 1 0 k¼1 where the functions ci;2kþ1 ðtÞ are given by the adjustment conditions (6) for component i and T6¼N or T ¼ N. Calculating the chirp-Wigner transform, W i ðt; f Þ½x for the two-component signal (20) gives W i ðt; f Þ½x ¼ A21 I i ðy1 ; y1 Þ þ A22 I i ðy2 ; y2 Þ þ A1 A2 I i ðy1 ; y2 Þ þ A1 A2 I i ðy2 ; y1 Þ.

(22)

It follows from Eq. (21) that for i ¼ 1,2, and j ¼ 1,2 the integral Ii(Z,x) is related to the chirp-Wigner transform of each signal component by I i ðyj ; yj Þ ¼ W i ðt; f Þ½xj .

ARTICLE IN PRESS 2992

L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002

Thus, when i ¼ 1, expression (22) can be written as W 1 ðt; f Þ½x ¼ A21 W 1 ðt; f Þ½x1  þ A22 W 1 ðt; f Þ½x2  þ A1 A2 I 1 ðy1 ; y2 Þ þ A1 A2 I 1 ðy2 ; y1 Þ.

(23)

The first term in expression (23) is the chirp-Wigner transform of the first signal component x1(t), completely adjusted to the first component. Therefore, as was shown in Section 3, this term can be expressed in terms of the Dirac function (or sinc function for T6¼N). The second term in expression (23) is the chirpWigner transform of the second component x2(t) adjusted to the first component, and the third and fourth terms are cross terms. Therefore, one can rewrite Eq. (23) as W 1 ðt; f Þ½x ¼ A21 dðf  f 1 ðtÞÞ þ A22 W 1 ðt; f Þ½x2  þ A1 A2 I 1 ðy1 ; y2 Þ þ A1 A2 I 1 ðy2 ; y1 Þ,

(24)

where f1(t) is the instantaneous frequency of the first component. Similarly, the chirp-Wigner transform of the signal adjusted to the second component can be written as W 2 ðt; f Þ½x ¼ A21 W 2 ðt; f Þ½x1  þ A22 dðf  f 2 ðtÞÞ þ A1 A2 I 2 ðy1 ; y2 Þ þ A1 A2 I 2 ðy2 ; y1 Þ,

(25)

where f2(t) is the instantaneous frequency of the second component. Thus, provided that the other three terms in Eqs. (24) and (25) do not interfere with the first term, the chirpWigner transform Wi(t, f)[x] will exactly follow the instantaneous frequency of the component xi(t) without amplitude errors, just as it does in the single-component case. When T6¼N, the Dirac function is replaced by the sinc function so that expressions (24) and (25) becomes, respectively: W 1 ðt; f Þ½x ¼ 2TA21 sincð2pðf  f 1 ðtÞÞÞ þ A22 W 1 ðt; f Þ½x2  þ A1 A2 I 1 ðy1 ; y2 Þ þ A1 A2 I 1 ðy2 ; y1 Þ,

(26)

W 2 ðt; f Þ½x ¼ A21 W 2 ðt; f Þ½x1  þ 2TA22 sincð2pðf  f 2 ðtÞÞÞ þ A1 A2 I 2 ðy1 ; y2 Þ þ A1 A2 I 2 ðy2 ; y1 Þ.

(27)

Provided that there is no interference from the last three terms, the chirp-Wigner transform will look like a sinc function along the curve of instantaneous frequency of the appropriate component. The height of this function above the instantaneous frequency curve will be 2TA2i . This allows an effective estimate to be made of the amplitude of the component i. It is obvious that the above analysis could be easily extended to the multi-component case, i.e. i42. Let us consider an important case from practical applications: the non-stationary two-component shaft passing vibrations generated by rotating shaft during machinery start up with nonlinear polynomial change of the shaft speed. Let us assume that both shaft harmonics (i.e. the fundamental and the second harmonics) have unity amplitude and the polynomial variation in the instantaneous phase of degree 5 given by y1 ðtÞ ¼ 2p 0:1t5 ; y2 ðtÞ ¼ 2p 0:2t5 ; 2ptp3. The graph of the instantaneous frequencies of the two harmonics is shown in Fig. 11. The solid and broken lines are the curves of the instantaneous frequency of the first and second harmonics.

Fig. 11. The instantaneous frequency curves of the harmonics.

ARTICLE IN PRESS L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002

2993

To estimate the transform we have used the adjustment conditions (6) to the first and second harmonics, respectively: c1;3 ðtÞ ¼ 3t2 =4;

c1;5 ðtÞ ¼ 1=32,

c2;3 ðtÞ ¼ 3t2 =2;

c2;5 ðtÞ ¼ 1=16.

We take T ¼ 2 for the duration of the interval of integration in the transform. Fig. 12 displays the chirp-Wigner transform for the signal completely adjusted to the first harmonic. The pronounced front ridge in the time–frequency plane corresponds to the first term in Eq. (26), i.e. to the instantaneous frequency of the first harmonic. There is no obvious ridge corresponding to the second term in Eq. (26), i.e. to the instantaneous frequency of the second harmonic. Fig. 13 displays the chirp-Wigner transform for the signal completely adjusted to the second harmonic. The pronounced rear ridge in the time–frequency plane corresponds to the second term in Eq. (27), i.e. to the instantaneous frequency of the second harmonic. There is no obvious ridge corresponding to the first term in Eq. (27), i.e. to the instantaneous frequency of the first harmonic.

Fig. 12. The chirp-Wigner transform adjusted to the first harmonic.

Fig. 13. The chirp-Wigner transform adjusted to the second harmonic.

ARTICLE IN PRESS 2994

L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002

Fig. 14. The Wigner distribution of the signal.

Fig. 15. The time cross-section of the chirp-Wigner transform at t ¼ 2.5 s adjusted to the first harmonic.

These results show that for two-component signal where each component is a higher order chirp the chirpWigner transform adjusted to a single component will follow the instantaneous frequency of the component without amplitude errors. Fig. 14 shows the classical Wigner distribution of the signal (20). It can be seen that it is not possible to make out the instantaneous frequency curves of either harmonic and estimate harmonic amplitudes without errors. Thus, the chirp-Wigner transform has a major advantage over the classical Wigner distribution for amplitude estimation of the higher order multi-component chirps. Let us further investigate the chirp-Wigner transforms adjusted to the appropriate harmonics. For the first harmonic the front ridge is of interest in Fig. 12. To determine to what extent this ridge corresponds to the sinc function along the instantaneous frequency, time cross-sections of the transform are examined and compared with the sinc function in Eq. (26). The solid curves in Figs. 15 and 16 are time cross-sections of the chirp-Wigner transform adjusted to the first harmonic at t ¼ 2.5 and 3 s respectively, plotted against frequency. The broken curves are the sinc functions centered at the instantaneous frequency of the first harmonic at t ¼ 2.5 and 3 s respectively. The close match between the solid and broken curves in Figs. 15 and 16 shows that the chirp-Wigner transform is providing an effective amplitude estimate of the first harmonic along the instantaneous frequency. The time cross-section of the classical Wigner distribution is displayed in Fig. 17 at t ¼ 3 s. Notice that it is unsuitable for determining the amplitude of either harmonic. The broken curves are sinc functions corresponding to the instantaneous frequencies of each harmonic.

ARTICLE IN PRESS L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002

2995

Fig. 16. The time cross-section of the chirp-Wigner transform at t ¼ 3 s adjusted to the first harmonic.

Fig. 17. The time cross-section of the Wigner distribution at t ¼ 3 s.

Fig. 18. The time cross-section of the chirp-Wigner transform at t ¼ 3 s adjusted to the second harmonic.

Next, consider the chirp-Wigner transform adjusted to the second harmonic. For the second harmonic the rear ridge is of interest in Fig. 13. Fig. 18 displays the time cross-section of the chirp-Wigner transform at t ¼ 3 s adjusted to the second harmonic. The broken curve is a sinc function corresponding to the instantaneous frequency of the second harmonic, which is f ¼ 81 Hz at this time. The figure indicates the close correspondence between the chirp-Wigner

ARTICLE IN PRESS 2996

L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002

transform and the sinc function. It follows that the chirp-Wigner transform adjusted to the second harmonic allows amplitude estimation of this harmonic without errors. 7. Final remarks Taking into account the results presented in Sections 2–6, the proposed chirp-Wigner transform could be used for processing transient signals with any non-polynomial time variation of instantaneous frequency. This processing could be done by the following two-step adjustment procedure: (i) evaluate from independent measurement (e.g. from frequency sensors) the time variation of the instantaneous frequency of the signal and approximate this variation using the polynomial approximation; estimate an order and parameter values of the approximation; (ii) develop an appropriate version of the chirp-Wigner transform (3) with adjustment conditions (6) according to the obtained approximation of the instantaneous frequency and employ it for processing. The proposed adaptation could be done on-line (i.e. in real time) or off-line. Thus, to use the transform (3), one needs to know a priori or from independent measurement the instantaneous frequency variation of the signal in time. The empirical mode decomposition (EMD) method is proposed in Ref. [20] for processing of non-stationary nonlinear signals. Novel comparison of the chirp-Wigner transform and EMD method for higher order chirps was brought to authors’ attention by the paper reviewer. Future investigation will be performed for solution of this task. 8. Conclusions 1. The new technique, the time–frequency chirp-Wigner transform is proposed and investigated in Ref. [19]. The chirp-Wigner transform is a generalization of the higher order local polynomial Wigner distribution proposed by V. Katkovnik [16,17] for the case of the known nonlinear polynomial variation of the instantaneous frequency and its derivatives. This technique is further investigated in this paper for the general case of higher order chirps, i.e. non-stationary signals with any nonlinear polynomial variation of the instantaneous frequency in time. 2. Analytical and numerical comparison of the chirp-Wigner transform and the classical Wigner distribution was performed for processing of single-component and multi-component higher order chirps. 3. It is shown that the chirp-Wigner transform uses approximately half the number of coefficients that are required to specify the higher order chirps, and that it is real valued and bilinear. 4. The adjustment of the chirp-Wigner transform is investigated for the general case of the higher order chirps. The new analytical relationships between time-varying parameters of the chirp-Wigner transform and signal parameters were found for this general case. These relationships provide transform adjustment. 5. It is shown for the general case of single component higher order chirps that the chirp-Wigner transform with adjustment has an essential advantage in comparison with the traditional Wigner distribution: the chirp-Wigner transform ideally follows the nonlinear polynomial frequency variation without amplitude errors. 6. It is shown for multi-component signal where each component is a higher order chirp that the chirp-Wigner transform adjusted to a single component will follow the instantaneous frequency of the component without amplitude errors. It is also shown that the classical Wigner distribution is unable to estimate component amplitudes of multi-component higher order chirps. Acknowledgement The authors are very thankful to paper reviewer for the comments, propositions and reference. Appendix A Note that since y(t) is assumed to be a polynomial of degree n in the variable t, it follows that the expression fðtÞ ¼ yðt þ ðt=2ÞÞ  yðt  ðt=2ÞÞ is also a polynomial of degree at most n in the variable t. Using Taylor’s

ARTICLE IN PRESS L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002

2997

formula to expand f(t) as a function of t about zero, we can write fðtÞ ¼ fð0Þ þ fð1Þ ð0Þt þ fð2Þ ð0Þt2 þ fð3Þ ð0Þt3 þ    þ fðnÞ ð0Þtn . We note that Eq. (28) is exact since f(t) is a polynomial of degree n. Then, if yðtÞ ¼ straightforward to show that the derivative of order m of f(t) is given by

 n  X kðk  1Þ . . . ðk  m þ 1Þak  tkm t km m fðmÞ ðtÞ ¼ t þ  ð1Þ t  . 2 2 2m k¼m

Pn

(28) k

k¼1 ak t

, it is

(29)

Putting t ¼ 0 in Eq. (29) gives for m ¼ 1; . . . ; n fðmÞ ð0Þ ¼

n X kðk  1Þ . . . ðk  m þ 1Þak k¼m

¼

8 n P > < > :

k¼m

2m kðk1Þ...ðkmþ1Þak km t 2m1

ð1  ð1Þm Þtkm ðmÞ

¼ y2m1ðtÞ ;

m odd; m even:

0

ðkÞ

Note that y ðtÞ ¼ V ðkÞ ðtÞ; k ¼ 2; . . . ; n: Substituting the above into Eq. (28), gives the following: r r   X t t X yð2kþ1Þ ðtÞ 2kþ1 V ð2kþ1Þ ðtÞ 2kþ1 y tþ t t ¼ f ðtÞt þ . y t ¼ i 2k 2k 2 2 k¼0 2 ð2k þ 1Þ! k¼1 2 ð2k þ 1Þ! Appendix B The first part of this appendix presents the derivation of formula (12). Suppose that the signal has polynomial variation of the instantaneous phase of degree 4 and that the duration is in the finite interval [T, T]. Calculating the Wigner distribution for (1), using the notation in formula (11) gives



 Z T Z T c3 ðt; lÞ 3 2 2 W ðt; f Þ ¼ A t exp j 2pðf i ðtÞ  f Þt þ 2p expðjðat3 þ btÞÞ dt, dt ¼ A 3 T T where l ¼ ðf 0 ; c20 ; c30 ; c40 Þ is the constant parameter vector of the signal, c3 ðt; lÞ ¼ ðc30 þ 3c40 t=4Þ, a ¼ 2pðc3 ðt; lÞ=3Þ, b ¼ 2pðf i ðtÞ  f Þ. Since the integrand is Hermitian symmetric one has that Z T Z T expðjðat3 þ btÞÞ dt ¼ 2Re expðjðat3 þ btÞÞ dt. 0

T 3

Replacing the factor expðjat Þ by its power series expansion at zero, and interchanging the order of integration and summation (permissible because of uniform convergence) one has when f af i ðtÞ that ) Z T Z T (X Z T 1 3 n ðjat Þ expðjðat3 þ btÞÞ dt ¼ 2Re expðjat3 Þ expðjbtÞ dt ¼ 2Re expðjbtÞ dt n! T 0 0 n¼0 ¼ 2Re

1 X ðjaÞn n¼0

¼2

1 ðGð3n þ 1Þ  Gð3n þ 1; jbTÞÞ n! ðjÞ3nþ1 b3nþ1

1 n X a n¼0

1 ImðGð3n þ 1; jbTÞÞ. n! b3nþ1

In the last step we used the fact that the function Gð3n þ 1Þ is real valued. When f ¼ f i ðtÞ, a different integral must be evaluated. Since in this case b ¼ 0, one has Z T Z 1 1 1 X X X ðjaÞn T 3n ðjaÞn T 3nþ1 ð1Þn a2n T 6nþ1 3 ¼2 . expðjðat þ btÞÞ dt ¼ 2Re t dt ¼ 2Re n! 0 n! ð3n þ 1Þ ð2nÞ! ð6n þ 1Þ T n¼0 n¼0 n¼0

ARTICLE IN PRESS 2998

L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002

This result could also be obtained from the result for the f af i ðtÞ case by using l’Hoˆpital’s rule to evaluate the limit of the expression ð1=b3nþ1 ÞImðGð3n þ 1; jbTÞÞ, as b ! 0. Next we present a derivation of formula (13). Suppose that the signal has polynomial variation of the instantaneous phase of degree 6 and that the duration is in the finite interval [T, T]. Then calculating the Wigner distribution for (1), using the notation in formula (11) gives



 Z T c3 ðt; lÞ 3 c5 ðt; lÞ 5 2 t þ 2p t W ðt; f Þ ¼ A exp j 2pðf i ðtÞ  f Þt þ 2p dt 3 5 T Z T expðjðat5 þ ct3 þ btÞÞ dt, ¼ A2 T

where l ¼ ðf 0 ; c20 ; c30 ; c40 ; c50 ; c60 Þ is the constant parameter vector of the signal,



 c30 þ 3c40 t þ 6c50 t2 þ 10c60 t3 c50 þ 5c60 t c3 ðt; lÞ ¼ ; c5 ðt; lÞ ¼ 16 4 c5 ðt; lÞ c3 ðt; lÞ ; b ¼ 2pðf i ðtÞ  f Þ; c ¼ 2p . a ¼ 2p 5 3 If the phase is of degree 5, then c60 is zero. The power series expansion of the exponential function about zero, together with the Binomial theorem gives the following: 1 1 n  X n k nk ðu þ vÞn X 1X ¼ u v . expðu þ vÞ ¼ (30) n! k¼0 k n! n¼0 n¼0 This formula is used to expand the integrand factor expðjðat5 þ ct3 ÞÞ as an infinite series in the following calculation. The interchange of the order of integration and summation is permissible because of uniform convergence of the series. When f af i ðtÞ, one has the following after transformations Z T Z T expðjðat5 þ ct3 þ btÞÞ dt ¼ 2Re expðjðat5 þ ct3 ÞÞ expðjbtÞ dt T 0 ! 1 n X n 1X ak cnk ¼2 ð1Þk 3nþ2kþ1 ImðGð3n þ 2k þ 1; jbTÞÞ. n! k¼0 k b n¼0 In the last step we used the fact that Gð3n þ 2k þ 1Þ is real valued. When f ¼ f i ðtÞ, a different integral must be evaluated. Since in this case b ¼ 0, one has the following calculation: ! Z T Z T 1 nX X n j n 5 3 expðjðat þ ct þ btÞÞ dt ¼ 2Re t3nþ2k dt ak cnk n! k T 0 n¼0 k¼0 ! 1 2n nX 6nþ2kþ1 X 2n ð1Þ T . ¼2 ak c2nk ð2nÞ! ð6n þ 2k þ 1Þ k n¼0 k¼0 This result could also be obtained from the result for the f af i ðtÞ case by using l’Hoˆpital’s rule to evaluate the limit of the expression ð1=b3nþ2kþ1 ÞImðGð3n þ 2k þ 1; jbTÞÞ, as b ! 0. Formula (14) for the chirp-Wigner transform is obtained using formula (12) for the Wigner distribution and the following formula



 Z 1 Dc3 3 t exp j 2pðf i ðtÞ  f Þt þ 2p W c ðt; f ; l0 ðtÞÞ ¼ A2 dt, 3 1 where Dc3 ¼ c3 ðt; lÞ  c3 ðt; l0 Þ, since the integral on the right-hand side is the same one that is evaluated to get the Wigner distribution (12).

ARTICLE IN PRESS L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002

2999

Formula (15), the chirp-Wigner transform, is obtained using formula (13) for the Wigner distribution and the following formula:



 Z 1 Dc3 3 Dc5 5 t þ 2p t exp j 2pðf i ðtÞ  f Þt þ 2p W c ðt; f ; l0 ðtÞÞ ¼ A2 dt, 3 5 1 where Dc3 ¼ c3 ðt; lÞ  c3 ðt; l0 Þ, since the integral on the right hand side is the same one that is evaluated to get the Wigner distribution (13). Appendix C This appendix presents the derivation of formulae (16)–(19). Suppose that the signal has polynomial variation of the instantaneous phase of degree 4 and the signal duration is infinite. Then in view of the calculation of the Wigner distribution given in Appendix B one must evaluate the integral: Z 1 expðjðat3 þ btÞÞ dt, (31) 0

where a ¼ 2p

c3 ðt; lÞ b ¼ 2pðf i ðtÞ  f Þ. 3

We cannot proceed as before, by expanding the integrand as an infinite series, because when we interchange summation and integration, the resulting integrals do not converge. Therefore, we transform the integral by a contour integration. Let C be the contour in the complex plane given by C ¼ C1 þ C2 ¼ C3, shown in Fig. 19, where C1 is the line segment from the origin to the point R on the real axis, C2 is the arc of the circle centered at the origin from R to R expðjp=2nÞ, n is a positive integer, C3 is the line segment from R expðjp=2nÞ to the origin.ref SHAPE Let C be the contour above with n ¼ 3. Define the function of a complex variable f(z) by f ðzÞ ¼ expðjðaz3 þ bzÞÞ. Then since f(z) is analytic everywhere in the complex plane, we can apply Cauchy’s theorem to conclude that Z Z Z Z f ðzÞ dz ¼ f ðzÞ dz þ f ðzÞ dz þ f ðzÞ dz ¼ 0. (32) C

C1

C2

C3

Im

C3 C2 π/2n C1

R

R e Re

Fig. 19. Complex plane for integration.

ARTICLE IN PRESS L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002

3000

R RR The integral along C1 is given by C1 f ðzÞ dz ¼ 0 expðjðat3 þ btÞÞ dt, and the integral along C2 is given by Z Z 0 f ðzÞ dz ¼ expðjp=6Þ expðjðaðexpðjp=6ÞtÞ3 þ b expðjp=6ÞtÞÞ dt C3 R Z R expðat3 þ jb expðjp=6ÞtÞ dt. ¼  expðjp=6Þ 0

R One can show using Jordan’s lemma that C2 f ðzÞ dz ! 0, as R ! 1. Therefore, taking the limit in Eq. (24) as R ! 1 gives the following result: Z 1 Z 1 3 expðjðat þ btÞÞ dt ¼ expðjp=6Þ expðat3 þ jb expðjp=6ÞtÞ dt. 0

0

Now we can evaluate (23) by developing the integrand as a power series as in Appendix B. Note that the resulting integrals will converge because of the expðat3 Þ factor: ( ) Z 1 Z 1 1 X ðjb expðjp=6ÞtÞn 3 3 expðjðat þ btÞÞ dt ¼ expðjp=6Þ expðat Þ dt n! 0 0 n¼0 Z 1 X ðjb expðjp=6ÞÞn 1 n t expðat3 Þ dt ¼ expðjp=6Þ n! 0 n¼0

 1 n X b j expðj2ðn þ 1Þp=3Þ 1 nþ1 ¼  G . n! 3 3aðnþ1Þ=3 n¼0 It follows from Hermitian symmetry of the integrand that the integral over the whole real line is given by

 Z 1 Z 1 1 n 2X b sinð2ðn þ 1Þ=3Þ 1 nþ1 expðjðat3 þ btÞÞ dt ¼ 2Re expðjðat3 þ btÞÞ dt ¼ G . 3 n¼0 n! 3 aðnþ1Þ=3 1 0 Finally, the formula for the Wigner distribution in expression (16) follows from:



 Z 1 c3 ðt; lÞ 3 2 t exp j 2pðf i ðtÞ  f Þt þ 2p dt W ðt; f Þ ¼ A 3 1 Z 1 expðjðat3 þ btÞÞ dt ¼ A2 1 Z 1 2 ¼ 2A Re expðjðat3 þ btÞÞ dt. 0

Suppose next, that the signal has polynomial variation of the instantaneous phase of degree 6 and that the signal duration is infinite. Then in view of the calculation of the Wigner distribution given in Appendix B one must evaluate the integral: Z 1 expðjðat5 þ ct3 þ btÞÞ dt, (33) 1

where a ¼ 2p

c5 ðt; lÞ ; 5

b ¼ 2pðf i ðtÞ  f Þ;

and

c ¼ 2p

c3 ðt; lÞ . 3

As in the n ¼ 6 case above we cannot expand the integrand as a power series, because the resulting integrals will not converge. Therefore, we transform the integral by a contour integration. Let C be the contour above with n ¼ 5. Define the function of a complex variable f(z) by f ðzÞ ¼ expðjðaz5 þ cz3 þ bzÞÞ.

(34)

ARTICLE IN PRESS L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002

3001

Then since f(z) is analytic everywhere in the complex plane, we can apply Cauchy’s theorem to conclude that Z Z Z Z f ðzÞ dz ¼ f ðzÞ dz þ f ðzÞ dz þ f ðzÞ dz ¼ 0. (35) C

C1

C2

C3

The integral along C1 is given by Z Z R f ðzÞ dz ¼ expðjðat5 þ ct3 þ btÞÞ dt C1

0

and the integral along C2 is given by: Z Z 0   f ðzÞ dz ¼ expðjp=10Þ exp jðaðexpðjp=10ÞtÞ5 þ cðexpðjp=10ÞtÞ3 þ b expðjp=10ÞtÞ dt C3 R Z R expðat5 þ c expðj8p=10Þt3 þ b expðj6p=10ÞtÞ dt. ¼  expðjp=10Þ 0

R One can show using Jordan’s lemma that C2 f ðzÞ dz ! 0, as R ! 1. Therefore, taking the limit as R ! 1 gives the following result: Z Z 1     1 5 3 exp jðat þ ct þ btÞ dt ¼ exp jp=10 expðat5 þ c expðj8p=10Þt3 þ b expðj6p=10ÞtÞ dt. 0

0

Now we can evaluate the integrand as a power series. Note that the resulting integrals will converge because of the exp(3at35) factor: Z 1 expðjðat5 þ ct3 þ btÞÞ dt 0 Z 1 ¼ expðjp=10Þ expðat5 Þ expðc expðj8p=10Þt3 þ b expðj6p=10ÞtÞ dt 0 ! ( ) Z 1 1 n X n  k 1X 5 3 nk ¼ expðjp=10Þ expðat Þ b expðj6p=10Þt ðc expðj8p=10Þt Þ dt n! k¼0 k 0 n¼0 !

 1 n X n k 1X 1 3n  2k þ 1 nk expðjð8n  2k þ 1Þp=10Þ ð3n2kþ1Þ=5 G . b c ¼ n! k¼0 k 5 5a n¼0 Finally, the formula for the Wigner distribution in Eq. (17) follows from:



 Z 1 c3 ðt; lÞ 3 c5 ðt; lÞ 5 t þ 2p t exp j 2pðf i ðtÞ  f Þt þ 2p dt W ðt; f Þ ¼ A2 3 5 1 Z 1 ¼ 2A2 Re expðjðat5 þ ct3 þ btÞÞ dt. 0

Formula (18), for the chirp-Wigner transform, is obtained using the formula for the Wigner distribution (16) and the following formula:



 Z 1 Dc3 3 2 0 W c ðt; f ; l ðtÞÞ ¼ A t exp j 2pðf i ðtÞ  f Þt þ 2p dt, 3 1 where Dc3 ¼ c3 ðt; lÞ  c3 ðt; l0 Þ, since the integral on the right-hand side is the same one that is evaluated to get the Wigner distribution (16). Formula (19), for the chirp-Wigner transform, is obtained using the formula for the Wigner distribution (17) and the following formula:



 Z 1 Dc3 3 Dc5 5 2 0 t þ 2p t W c ðt; f ; l ðtÞÞ ¼ A exp j 2pðf i ðtÞ  f Þt þ 2p dt; 3 5 1 where Dc3 ¼ c3 ðt; lÞ  c3 ðt; l0 Þ, since the integral on the right hand side is the same one that is evaluated to get the Wigner distribution (17).

ARTICLE IN PRESS 3002

L. Gelman, J.D. Gould / Mechanical Systems and Signal Processing 21 (2007) 2980–3002

Appendix D. Supplementary materials Supplementary data associated with this article can be found in the online version at doi:10.1016/ j.ymssp.2007.05.003.

References [1] E. Wigner, On the quantum correction for thermodynamic equilibrium, Physical Review 40 (5) (1932) 749–759. [2] W. Mecklenbrauker, F. Hlawatsch, The Wigner distribution-theory and application in signal processing, Elsevier Science, Amsterdam, 1997. [3] T.J. Wahl, J.S. Bolton, The application of the Wigner distribution to the identification of structure-borne noise components, Journal of Sound and Vibration 163 (1) (1993) 101–122. [4] M. Poletti, Linearly swept frequency measurements and the Wigner–Ville distribution, chapter 19 in book B. Boashash, Time–Frequency Signal Analysis (1992) 424–444. [5] M.S. Safizadeh, A. Lakis, Time frequency distributions and their application to machinery fault diagnosis, International Journal of Comadem 5 (2) (2002) 41–56. [6] J. Moss, J. Hammond, A comparison between the modified spectrogram and the pseudo Wigner–Ville distribution with and without modification, Mechanical Systems and Signal Processing 8 (3) (1994) 243–258. [7] J. Zou, J. Chen, A comparative study on time–frequency feature of cracked rotor by Wigner–Ville distribution and wavelet transform, Journal of Sound and Vibration 276 (2004) 1–11. [8] P. Flandrin, P. Rioul, Affine smoothing of the Wigner–Ville distribution, IEEE International Conference on Acoustics, Speech, Signal Processing, USA, New Mexico (1990) 2455–2458. [9] B. Forrester, Use of the Wigner–Ville distribution in helicopter fault detection, in: Proceedings, International Symposium Signal Processing Algorithms 89, USA, 1989, pp. 78–82. [10] L. Gelman, D. Adamenko, Usage of the Wigner–Ville distribution for vibro-acoustcal diagnostics of cracks, Proceedings of the Meeting of the Acoustical Society of America, Journal of the Acoustical Society of America 108 (5 Part 2) (2000). [11] B. Boashash, P. Black, An efficient real-time implementation of the Wigner–Ville distribution, IEEE Transactions on Acoustics, Speech and Signal Processing 35 (1987) 1611–1618. [12] B. Boashash, P. O’Shea, Polynomial Wigner–Ville distributions and their relationship to time-varying higher order spectra, IEEE Transactions on Signal Processing 42 (1) (1995) 216–220. [13] B. Boashash, B. Ristic, Polynomial WVD’s and time-varying polyspectra, in: B. Boashash, et al. (Eds.), Higher Order Statistical Processing, Longman Cheshire, 1993. [14] L. Stankovic, On the realization of the polynomial Wigner–Ville distribution for multicomponent signals, IEEE Signal Processing Letters 5 (7) (1998) 157–159. [15] B. Ristic, B. Boashash, Relationship between the polynomial and higher order Wigner–Ville distribution, IEEE Signal Processing Letters 2 (1996) 227–229. [16] V. Katkovnik, Local polynomial approximation of the instantaneous frequency: asymptotic accuracy, Signal Processing 52 (3) (1996) 343–356. [17] V. Katkovnik, Erratum to ‘‘Local polynomial approximation of the instantaneous frequency: asymptotic accuracy,’’ Signal Processing 59 (2) (1997) 251–252. [18] L. Stankovic, Local polynomial Wigner distribution, Signal Processing 59 (1997) 123–128. [19] L. Gelman, Adaptive time–frequency transform for non-stationary signals with nonlinear polynomial frequency variation, Mechanical Systems and Signal Processing 21 (2007) 2684–2687. [20] N. Huang, Z. Shen, S. Long. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proceedings of the Royal Society, UK, 1998, vol. 454, pp. 903–995.