Effect of interpolation on PWVD computation and instantaneous frequency estimation

Signal Processing 84 (2004) 107 – 116

www.elsevier.com/locate/sigpro

Eect of interpolation on PWVD computation and instantaneous frequency estimation S. Chandra Sekhar, T.V. Sreenivas∗ Department of Electrical Communication Engineering, Indian Institute of Science, Bangalore 560 012, Karnataka, India Received 15 February 2002; received in revised form 25 July 2003

Abstract The peak of the polynomial Wigner–Ville distribution (PWVD) can be used for estimating the instantaneous frequency (IF) of monocomponent polynomial phase signals. However, the PWVD kernel, optimized to yield a time–frequency distribution (TFD) localized along the IF, comprises of fractional-time-sampled signals. When implemented in a discrete-time scenario, this calls for signal interpolation. We study three interpolation schemes—linear, cubic polynomial and sinc and derive expressions for the variance of the interpolated samples in the presence of noise. In representing nonstationary signals using the PWVD, the instantaneous energy content of noise auto-terms and signal-noise cross-terms is found to be the least for linear interpolation scheme. For polynomial IF estimation using the peak of the PWVD, it was found that linear interpolation is a computationally e:cient way of obtaining reasonably good estimates at low signal-to-noise ratios (SNRs). For high SNRs, sinc interpolation outperforms the other two schemes. Similar results were found when the experiment was extended to sinusoidal IF signals also. ? 2003 Elsevier B.V. All rights reserved. Keywords: Interpolation; Polynomial Wigner–Ville distribution; Instantaneous frequency

1. Introduction Most information carrying signals such as speech, radar and sonar signals have time-varying spectral characteristics and cannot be eectively described by Fourier analysis. Joint time–frequency analysis has been used as a tool to characterize such nonstationary signals by describing their spectral characteristics as a function of time [6]. Of particular interest is the problem of instantaneous frequency (IF) estimation, which attempts to provide a high-resolution

∗

Corresponding author. Tel.: +91-80-293-2285; fax: +9180-360-0683. E-mail address: [email protected] (T.V. Sreenivas).

description of the time-varying spectral information. Such a description is often useful in several applications as in [8]. Several methods have been devised to estimate the IF of monocomponent phase signals [3–5]. The Wigner–Ville distribution (WVD), a member of the bilinear time–frequency (TF) distributions [6] has been shown to be eective in localizing monocomponent phase signals with linear frequency modulation. Hence, the IF can be obtained by picking the peak of the WVD. However, peak-picking the WVD of polynomial phase signals will yield biased estimates of the IF. To circumvent the problem, polynomial Wigner–Ville distribution (PWVD) has been introduced which is a generalization of the WVD to higher-order polynomial phase signals [1]. A pth

0165-1684/$ - see front matter ? 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2003.07.015

108

S. Chandra Sekhar, T.V. Sreenivas / Signal Processing 84 (2004) 107 – 116

order PWVD best localizes, in the joint time– frequency plane, a polynomial phase signal of order less than or equal to p. As a result, the IF of such signals can be easily obtained by picking the peak of the PWVD. However, the theory of PWVD was originally developed and the optimal kernel was designed for the continuous-time case. Maintaining the optimality of the distribution in discrete-time implementation requires interpolation to Hnd the signal amplitude at nonintegral instants of time. There are several methods for interpolation such as linear, cubic polynomial, bandlimited sinc interpolation [7] etc. The eectiveness of interpolation will also depend on additive noise present in the signal. The focus of this paper is to study the eect of dierent interpolation schemes on the PWVD of noisy signals.

0:675; and c3 = −c3 ≈ 0:675. The design procedure for higher-order PWVD is illustrated in [2]. Consider a discrete-time pth order polynomial phase signal z[n] over a Hnite duration [0; N − 1] given by z[n] = Ae{j


p

i=0


p

i=0

ai t i }

;

 q=2   m

(1)

where the coe:cients {ai ; i = 0; 1; : : : ; p} are real valued. The IF of such a signal is a polynomial of order p − 1 and is given by

(4)

P[n; f)

2. Polynomial Wigner–Ville distribution

z(t) = Ae{j

:

The sampling period is normalized to unity. We would like to estimate the IF of this signal using the peak of its PWVD. This is clearly a nonlinear estimator. The discrete-time equivalent of the PWVD for this signal, denoted by P[n; f) 1 (maintaining the optimality of the distribution) can be expressed as

=

The PWVD was designed to obtain an optimum time–frequency representation for a polynomial phase signal of the form

ai ni }

 ∗

z[n + ci m]z [n + c−i m] e−j2 fm : (5)

i=1

Clearly, the distribution depends on the value of the signal obtained at nonintegral instants of time. Since the signal samples are available only at discrete instants of time, we need to interpolate to obtain the signal at the fractional time instants. There are many schemes for interpolation and the computed time– frequency distribution and the resulting IF estimate depends on the scheme of interpolation employed.

p

1  f(t) = iai t i−1 : 2

(2)

i=1

The PWVD of the signal z(t), denoted by P(t; f) is deHned as P(t; f)  =

+∞

−∞

 q=2   z(t + ci )z ∗ (t + c−i ) e−j2 f d; i=1

(3)

where q is an even integer that indicates the order of the nonlinearity of the PWVD kernel. The coe:cients {ci ; c−i ; i = 1; 2; : : : ; q=2}, are obtained such that the PWVD is real and is localized around the IF for a polynomial phase signal. The WVD is a member of this class with parameters q = 2 and c1 = −c−1 = 0:5. For the sixth-order PWVD (q = 6), the optimal kernel coe:cients are c1 = −c−1 ≈ −0:85 and c2 = −c−2 ≈

3. Variance of error in interpolation We consider three interpolation schemes, namely, linear interpolation, cubic polynomial interpolation and bandlimited sinc interpolation. We discuss the relative performance of these methods in performing interpolation in the presence of noise by deriving the expressions for bias and variance of the estimates as a function of the noise variance. 3.1. Linear interpolation In this scheme, we Ht a straight line through two points and use it to obtain the values at instants in between the two points.

1

The notation means that n is discrete and f is continuous.

S. Chandra Sekhar, T.V. Sreenivas / Signal Processing 84 (2004) 107 – 116

Let n1 and n2 be two time instants and it is desired to interpolate the signal z[n] at a time t (referred to as the query point in this paper) between n1 and n2 , i.e., n1 ¡ t ¡ n2 . The estimate of the signal at time t is given by   t − n1 {z[n2 ] − z[n1 ]} + z[n1 ] z(t) = (6) n2 − n 1 with the boundary conditions satisHed. Assume that the samples z[n] are corrupted with white noise samples, w[n], of zero mean and variance w2 . The observed signal, y[n] is given by y[n]=z[n]+ w[n]. Using these samples for performing interpolation, we have   t − n1 {y[n2 ] − y[n1 ]} + y[n1 ]: y(t) = (7) n2 − n1 The contribution of noise in the above estimate, denoted by e(t), is given by   t − n1 {w[n2 ] − w[n1 ]} + w[n1 ]: (8) e(t) = n2 − n1 Clearly, the estimator is unbiased, i.e., E{y(t)} = z(t);

(9)

where E{ } denotes the expectation operator. 2 , is The variance of the estimate, denoted by y(t) given by 2 = w2 {s2 (t) + (1 − s(t))2 }; y(t)

point. If t is any query point, then the samples in the neighbourhood are chosen at instants n1 = t − 1; n2 = t; n3 = t; n4 = t + 1. We Ht a smooth cubic polynomial z(t) = a3 t 3 +  2 a2 t + a1 t + a0 through these points. 2 The coe:cients of the polynomial can be found uniquely by solving a set of linear equations using the information at the time instants n1 ; n2 ; n3 and n4 . For the noisy signal y(t), let the interpolation coe:cients be a3 ; a2 ; a1 and a0 , respectively. The coe:cient estimates are aected by noise. It can be easily shown that the coe:cient estimates are unbiased and therefore, E{y(t)} = z(t):

3.2. Cubic polynomial interpolation

(11)

The variance of the estimate is given by 2 y(t) = w2

6 

vi t i ;

(12)

i=0

where the coe:cients vi are given by the following equations: v0 =

4 

4;2 j ;

(13)

j=1

v1 = 2

4 

3; j 4; j ;

(14)

j=1

(10)

where s(t) = (t − n1 )=(n2 − n1 ). This equation is easily interpreted as follows: At the given instants n1 and n2 , the variance of the estimate is only due to additive noise. At other instants, in between n1 and n2 , the variance is a quadratic function of the distance between the point of interpolation and the known points. The time dependent slope factor s(t) is symmetric half-way through the points n1 and n2 . The estimate obtained by linear interpolation has least variance at this point. This is true of any interval between the sampling instants. In the present work, n1 and n2 are chosen as successive time instants.

109

v2 =

4 

{3;2 j + 22; j 4; j };

(15)

j=1

v3 = 2

4 

{2; j 3; j + 1; j 4; j };

(16)

j=1

v4 =

4 

{2;2 j + 21; j 3; j };

(17)

j=1

v5 = 2

4 

1; j 4; j ;

(18)

j=1

v6 =

4 

1;2 j :

(19)

j=1

In this scheme, we consider four samples in the neighbourhood of the query point and Ht a cubic polynomial which is used to estimate the signal at the query

2

The prime on the coe:cients indicates that noise is absent.

110

S. Chandra Sekhar, T.V. Sreenivas / Signal Processing 84 (2004) 107 – 116

Here, i; j given by  3 n1  3  n2 =  n3  3 n34

is the (i; j)th element of the matrix, , n21

n1

n22 n23 n24

n2 n3 n4

1

The contribution of noise to the error in the estimate is given by

−1

 1  1  1

e(t) = ;

(20)

sin( (t − n)) :

(t − n)

n=0

In all schemes discussed so far, we have found that the estimators are unbiased. The variance of the interpolated value depends on the noise variance and in general, has the form 2 y(t) = (t)w2 ;

0.5

0.5

0

0

-0.5

-0.5

-1 -1.5 -2 -2.5

-1

-2 -2.5

-3

-3 -3.5

(a)

58.4

58.6

time

58.8

-4 58

59

(b)

x 10

-3

SINC

0

-1.5

-3.5 58.2

5

Noise Factor (dB)

1

Noise Factor (dB)

Noise Factor (dB)

CUBIC

1

-4 58

(26)

where (t) is referred to as the noise factor and is determined by the interpolation technique. It is desired that (t) be as small as possible. Fig. 1 shows 2 a plot of (t) = (y(t) =w2 ) which is nothing but the variance of the estimate normalized with respect to the noise variance. (t) is a function of the time instant at which interpolation is performed and gives a measure of the amount of error involved. The plots are identical for every sampling interval except for sinc which over the entire observation interval has an envelope which varies as 1=t 2 . It is clear from

(22)

LINEAR

(24)

The variance of the estimator can be obtained as 2 N −1   sin( (t − n)) 2 = w2 : (25) y(t)

(t − n)

where the sampling period has been normalized to unity. A similar relation holds for y(t). In practice, we have discrete-time signals of Hnite duration only, in which case, we can write,

n=0

(23)

E{y(t)} = z(t):

For bandlimited signals, sinc interpolation is the optimum method of signal reconstruction from its sample sequence. Thus, the signals z(t) and its sample sequence, z[n] are related by ∞  sin( (t − n)) z(t) = z[n] ; (21)

(t − n) n=−∞

y[n]

sin( (t − n)) :

(t − n)

It is easy to show that the estimator is unbiased, i.e.,

3.3. Sinc interpolation

N −1 

w[n]

n=0

which arises in the solution for the coe:cients which are unknowns in a set of linear equations. The interrelations between n1 ; n2 ; n3 and n4 can be used to yield simpler expressions for the coe:cients; this is very straightforward and is omitted. For large values of n1 ; n2 ; n3 and n4 , singularity problems may arise in computing the matrix inverse.

y(t) =

N −1 

-5

-10

-15

58.2

58.4

58.6

time

58.8

59

-20 58

(c)

58.2

58.4

58.6

58.8

time

Fig. 1. Noise factor for dierent interpolation schemes: (a) linear, (b) cubic polynomial and (c) sinc.

59

S. Chandra Sekhar, T.V. Sreenivas / Signal Processing 84 (2004) 107 – 116

the Hgure that (t) = 1 at the sample instants. In between the samples, (t) ¡ 1, indicating decreased variance. The least variance is at the centre of the interval for all the techniques. Within an interval, the least variance is for linear interpolation, whereas sinc interpolation has most variance. Compared to linear and cubic, sinc has almost Lat variance within a sampling interval.

111

where Ne [n; m], under high-SNR assumption, is approximated by n  1  k1 ∗k1 Ne [n; m] ≈ z [n + ci m]z [n − ci m] i=1 n1  × ki [z −1 [n + ci m]e[n + ci m] i=1

+ z ∗−1 [n − ci m]e∗ [n − ci m]]: 4. PWVD computation and IF estimation In this section, we study the eects of interpolation in computing the PWVD in the presence of noise and its eect on IF estimation. 4.1. Power of the noise-terms in the PWVD

q=2  {z[n + ci m] + e[n + ci m]} i=1

×{z[n − ci m] + e[n − ci m]}∗ ;

(27)

where e[ · ] denotes the error due to interpolation and is a linear combination of the additive noise samples for the methods discussed. Rewriting, we get

m=0 = E{Ne [n; 0]Ne∗ [n; 0]} n  1  2(q−1) 2 ki w2 : = 2A

(31) (32)

When the lag is zero, no interpolation is required and the noise present is just the additive noise. For nonzero lags, this will also include error due to interpolation which depends on the type of interpolation scheme used and is nonstationary in nature. The power of the noise terms for nonzero lags is given by m=0 = E{Ne [n; m]Ne∗ [n; m]}  n 1  = A2q E ki [z −1 [n + ci m]e[n + ci m] i=1

+z

∗−1

[n − ci m]e∗ [n − ci m]]

n1  −1 × kj [z ∗ [n + ci m]e∗ [n + ci m]

Ky [n; m] = {z[n + ci m] + e[n + ci m]}k1 · · ·

j=1

×{z[n + cn1 m] + e[n + cn1 m]}kn1

+z

×{z ∗ [n − ci m] + e∗ [n − ci m]}k1 · · · ×{z ∗ [n − cn1 m] + e∗ [n − cn1 m]}kn1 ; (28) where n1 is the number of dierent coe:cients ci and ki ; i = 1; 2; 3; : : : ; n1 , is the multiplicity of each coe:cient ci . Using the binomial theorem and approximating the kernel by neglecting the product terms of more than one noise term, we get Ky [n; m] = Kz [n; m] + Ne [n; m];

Assuming i.i.d zero-mean complex Gaussian noise of variance w2 , it can be readily shown that the power of the noise terms (denoted by ) for a zero lag, i.e., m = 0, is given by

i=1

To compute the power of the noise-terms in the PWVD domain, we follow the analysis presented in [1]. The kernel for the noisy signal is given by Ky [n; m] =

(30)

(29)

−1

 [n − ci m]e[n − ci m]]

:

(33)

This involves expectations of products, in general, of the form e[n + ci m]e[n + cj m]. For the case of sinc interpolation, this can be very easily calculated as E{e[n + ci m]e[n + cj m]}   sin( (n + ci m − l))  2 =w

(n + ci m − l) l   sin( (n + cj m − l)) : ×

(n + cj m − l)

(34)

112

S. Chandra Sekhar, T.V. Sreenivas / Signal Processing 84 (2004) 107 – 116

50 45 40 35 30 25 20

50

(a)

100

150

200

50 45 40 35 30 25 20

250

55

Instantaneous Noise Energy (dB)

55

Instantaneous Noise Energy (dB)

Instantaneous Noise Energy (dB)

55

50

(b)

sample index

100

150

200

50 45 40 35 30 25 20

250

50

(c)

sample index

100

150

200

250

sample index

Fig. 2. Instantaneous noise energy for dierent methods of interpolation, for a phase signal with fourth-order polynomial IF, dotted: linear, dash-dot: cubic, solid: sinc. (a) SNR=5 dB, (b) SNR=10 dB, (c) SNR=20 dB.

45 40 35 30 25

50

100

150

200

50 45 40 35 30 25 20

250

sample index

55

Instantaneous Noise Energy (dB)

50

20

(a)

55

Instantaneous Noise Energy (dB)

Instantaneous Noise Energy (dB)

55

50

(b)

100

150

200

sample index

50 45 40 35 30 25 20

250

(c)

50

100

150

200

250

sample index

Fig. 3. Instantaneous noise energy for dierent methods of interpolation, for a phase signal with a sinusoidal IF, dotted: linear, dash-dot: cubic, solid: sinc. (a) SNR=5 dB; (b)=10 dB; (c) SNR=20 dB.

This can be used to compute the power of the noise terms for nonzero lags. For the case of linear and cubic interpolation, computation of the expectation of such product terms is by no means di:cult, but more tedious and is omitted. The purpose of this calculation is only to show that the noise introduced due to interpolation is nonstationary. To illustrate the nonstationarity of noise, we deHne a quantity called the instantaneous noise energy (INE), denoted by I(t) and deHned as  I(t) =

0

∞

|Py (t; f) − Pz (t; f)|2 df;

(35)

where Py (t; f) and Pz (t; f) are the PWVDs of y(t) and z(t), respectively. For purposes of implementation, a discrete-time, discrete-frequency equivalent of the above equation is used in the simulations.

4.2. IF estimation The IF of the signal z[n] is estimated using y[n] as ˆ = arg max Py [n; f): f[n] f

(36)

When analyzing signals with polynomial IF, if the order of the PWVD is matched to the polynomial order, then PWVDs localized along the IF can be obtained. Picking the peak of the PWVD can be used to perform IF estimation. However, if there is a mismatch between the signal and the order of the PWVD, peak-picking will result in biased IF estimates; but, the focus of the current work is to study the eect of dierent interpolation schemes on IF estimation. 4.3. Simulation experiment In this section, we report experimental results on the eect of interpolation in computing the PWVD and

S. Chandra Sekhar, T.V. Sreenivas / Signal Processing 84 (2004) 107 – 116 PWVD

I F estimate

0.3 0.2 0.1 50

100

150

200

0.4 0.3 0.2 0.1 0

250

0

squared error (dB)

frequency (Hz)

frequency (Hz)

0.4

0

50

sample index

250

0.3 0.2 0.1 50

100

150

200

0.2 0.1 50

100

150

200

50

0.2 0.1 200

250

200

250

200

250

200

250

-40 -60 -80 0

50

100

150

0

0.4 0.3 0.2 0.1 0

150

sample index

squared error (dB)

frequency (Hz)

0.3

100

-20

-100

250

0.5

150

0

sample index

0.4

100

-80

0

0.3

0

250

0.5

50

-60

sample index

squared error (dB)

frequency (Hz)

frequency (Hz)

200

0.4

sample index

frequency (Hz)

150

-40

-100

0.5

0.4

0

100

-20

sample index

0.5

0

I F error

0.5

0.5

113

50

sample index

100

150

200

-20 -40 -60 -80 -100

250

0

50

sample index

100

150

sample index

Fig. 4. Fourth-order polynomial IF estimation using the peak of PWVD (noise-free phase signal). The Hrst, second and third rows correspond to linear, cubic polynomial and sinc interpolation schemes, respectively. The Hrst, second and third columns correspond to the PWVD, the IF estimate obtained by peak-picking and the instantaneous squared error (dB) in the estimate, respectively. PWVD

I F estimate

0.3 0.2 0.1 50

100

150

200

0.4 0.3 0.2 0.1 0

250

0

squared error (dB)

frequency (Hz)

frequency (Hz)

0.4

0

50

sample index

250

0.3 0.2 0.1 50

100

150

200

0.4 0.3 0.2 0.1 0

250

50

100

150

200

250

0.2 0.1 150

200

250

150

200

250

200

250

200

250

-40 -60 -80 0

50

100

150

0

0.4 0.3 0.2 0.1 0

100

sample index

squared error (dB)

frequency (Hz)

0.3

100

50

-20

-100

0.5

0.4

sample index

0

sample index

0.5

50

-80

sample index

squared error (dB)

frequency (Hz)

frequency (Hz)

200

-60

0

sample index

frequency (Hz)

150

-40

sample index

0.4

0

100

-20

-100

0.5

0.5

0

I F error

0.5

0.5

50

100

150

sample index

200

250

-20 -40 -60 -80 -100

0

50

100

150

sample index

Fig. 5. Sinusoidal IF estimation using the peak of PWVD (noise-free phase signal). The Hrst, second and third rows correspond to linear, cubic polynomial and sinc interpolation schemes, respectively. The Hrst, second and third columns correspond to the PWVD, the IF estimate obtained by peak-picking and the instantaneous squared error (dB) in the estimate, respectively.

114

S. Chandra Sekhar, T.V. Sreenivas / Signal Processing 84 (2004) 107 – 116 PWVD

I F estimate

0.3 0.2 0.1 50

100

150

200

0.4 0.3 0.2 0.1 0

250

0

squared error (dB)

frequency (Hz)

frequency (Hz)

0.4

0

50

sample index

250

0.3 0.2 0.1 50

100

150

200

0.2 0.1 50

100

150

200

50

0.2 0.1 200

250

200

250

200

250

200

250

-40 -60 -80 0

50

100

150

0

0.4 0.3 0.2 0.1 0

150

sample index

squared error (dB)

frequency (Hz)

0.3

100

-20

-100

250

0.5

150

0

sample index

0.4

100

-80

0

0.3

0

250

0.5

50

-60

sample index

squared error (dB)

frequency (Hz)

frequency (Hz)

200

0.4

sample index

frequency (Hz)

150

-40

-100

0.5

0.4

0

100

-20

sample index

0.5

0

I F error

0.5

0.5

50

sample index

100

150

200

-20 -40 -60 -80 -100

250

0

50

sample index

100

150

sample index

Fig. 6. Fourth-order polynomial IF estimation using the peak of the PWVD for SNR=8 dB. The Hrst, second and third rows correspond to linear, cubic polynomial and sinc interpolation schemes, respectively. The Hrst, second and third columns correspond to the PWVD, the IF estimate obtained by peak-picking and the instantaneous squared error (dB) in the estimate, respectively. PWVD

0.3 0.2 0.1 50

100

150

200

0.4 0.3 0.2 0.1 0

250

0

squared error (dB)

frequency (Hz)

frequency (Hz)

0.4

0

50

sample index

250

0.3 0.2 0.1 50

100

150

200

250

0.2 0.1 50

100

150

200

250

0.2 0.1 200

250

150

200

250

200

250

200

250

-40 -60 -80 0

50

100

150

0

0.4 0.3 0.2 0.1 0

100

sample index

squared error (dB)

frequency (Hz)

0.3

150

50

-20

-100

0.5

100

0

sample index

0.4

sample index

-80

0

0.3

0

0.5

50

-60

sample index

squared error (dB)

frequency (Hz)

frequency (Hz)

200

0.4

sample index

frequency (Hz)

150

-40

-100

0.5

0.4

0

100

-20

sample index

0.5

0

I F error

I F estimate 0.5

0.5

50

100

150

sample index

200

250

-20 -40 -60 -80 -100

0

50

100

150

sample index

Fig. 7. Sinusoidal IF estimation using the peak of the PWVD for SNR=8 dB. The Hrst, second and third rows correspond to linear, cubic polynomial and sinc interpolation schemes, respectively. The Hrst, second and third columns correspond to the PWVD, the IF estimate obtained by peak-picking and the instantaneous squared error (dB) in the estimate, respectively.

S. Chandra Sekhar, T.V. Sreenivas / Signal Processing 84 (2004) 107 – 116 -16

-15 linear cubic sinc

Cumulative Mean Square Error (dB)

-18

-20

-22

-24

-26

0

5

10

15

-19

-23

-27

-31

-35

20

0

5

10

SNR (dB)

,

Linear best

15

20

SNR (dB)

,

Sinc best

(b)

,

,

(a)

,

,

Linear best

, Sinc best

,

Cumulative Mean Square Error (dB)

linear cubic sinc

-28

115

Fig. 8. Performance of the PWVD peak based IF estimator with dierent interpolation schemes: (a) fourth-order polynomial IF, (b) sinusoidal IF.

IF estimation by peak-picking the PWVD. Consider a real-valued phase signal with a fourth-order IF given by 4  n − (N=2) f[n] = 0:098 + 0:2940 ; (37) N=2 0 6 n 6 N −1. The coe:cients are chosen arbitrarily. For the purpose of simulation, N is chosen as 256. To illustrate the nonstationarity of noise in the PWVD domain, white Gaussian noise is added to the signal to achieve a desired SNR. The PWVD was computed using the analytic signal of the resulting noisy signal. The INE is computed for the three interpolation schemes for SNR=5,10 and 20 dB. Fig. 2 shows the computed value (in dB) as a function of time. As an extension, we consider sinusoidal IF, consisting of two cycles in the observation interval with frequency swinging between 0.0196 and 0:2706 Hz (chosen arbitrarily). The computed INE is shown in Fig. 3. From Figs. 2 and 3, it is clear that in the presence of noise, on an average, linear interpolation gives rise to a PWVD with least INE. In Figs. 4 and 5, we show the IF estimates obtained using the PWVD of the noise-free phase signal and the instantaneous squared error (dB) for fourth-order polynomial and sinusoidal IF, respectively, for dierent interpolation schemes. The Hrst column shows the

joint time–frequency PWVD plots. The plots in the second column show the IF estimated by peak-picking the PWVD. The frequency shown on the y-axis is normalized to 0:5 Hz. The third column shows plots of the instantaneous squared error in the IF estimate, in dB. The rows correspond to the interpolation schemes—linear, cubic polynomial and sinc interpolation, respectively. From the plots, it is clear that for the noise-free phase signal, the IF estimates for the polynomial IF are more or less identical whereas for sinusoidal IF, sinc interpolation outperforms the other two schemes. Similar plots for SNR = 8 dB are shown in Figs. 6 and 7. The desired SNR was obtained by adding white Gaussian noise to the phase signal. From these Hgures, we can conclude that in the presence of noise, all three interpolation schemes yield almost identical results. In such cases, the use of linear interpolation can result in considerable computational saving. In Fig. 5, where noise-free sinusoidal FM signals are considered, sinc interpolation yields a better PWVD and IF estimate than linear interpolation. Additional errors for the noisy signal (Fig. 7) are therefore due to interpolation in the presence of noise. To study the performance of the IF estimator as a function of SNR, we use a measure called the

116

S. Chandra Sekhar, T.V. Sreenivas / Signal Processing 84 (2004) 107 – 116

cumulative mean square error, denoted by Hned as =

N −1 M 1  (f[n] − fˆ l [n])2 ; MN

and de(38)

n=0 l=1

where M is the total number of realizations and fˆ l [n] is the IF estimate for the lth realization and f[n] is the actual IF. This measure is essentially an average of squared instantaneous estimation errors. Fig. 8 shows the performance plots corresponding to the same as a function of SNR. The results are obtained by averaging over 50 Monte-Carlo trials. The performance characteristics show that for low SNRs, linear interpolation oers least error. For very high SNRs, sinc interpolation turns out to be the best. This is expected because sinc interpolation is ideally suited for bandlimited signals. For high SNRs, the performance of cubic interpolation scheme is in between that of linear and sinc interpolation schemes. 5. Conclusion We have addressed the problem of signal interpolation in PWVD in the presence of noise; in particular, the eect of interpolation on IF estimation is studied. From the expressions of sample variance, linear interpolation is seen to be the best. In the computation of PWVD, it is found that linear interpolation results in noise auto-terms and signal-noise cross-terms

with least energy compared to the other two schemes. In estimating IF, at low SNRs, linear interpolation is the best. Its simplicity also provides a clear advantage. At high SNRs, however, sinc interpolation turns out to be the best. References [1] B. Barkat, B. Boashash, Instantaneous frequency estimation of polynomial FM signals using the peak of the PWVD: statistical performance in the presence of additive Gaussian noise, IEEE Trans. Signal Process. 47 (9) (September 1999) 2480–2490. [2] B. Barkat, B. Boashash, Design of higher order polynomial Wigner–Ville distributions, IEEE Trans. Signal Process. 47 (9) (September 1999) 2608–2611. [3] B. Boashash, Estimating and interpreting the instantaneous frequency of a signal—Part 1: Fundamentals, Proc. IEEE 80 (4) (April 1992) 519–538. [4] B. Boashash, Estimating and interpreting the instantaneous frequency of a signal—Part 2: Algorithms and applications, Proc. IEEE 80 (4) (April 1992) 539–568. [5] S. Chandra Sekhar, T.V. Sreenivas, Polynomial instantaneous frequency estimation using level-crossing information, Proceedings of the IEEE-EURASIP, Workshop on Nonlinear Signal and Image Processing (NSIP), Grado, Italy, 2003. [6] L. Cohen, Time-Frequency Analysis, Prentice-Hall, Englewood Clis, NJ, 1995. [7] A.V. Oppenheim, R.W. Schafer, Discrete-Time Signal Processing, Prentice-Hall of India, New Delhi, 1998. [8] D.C. Reed, A.M. Zoubir, B. Boashash, Aircraft Light parameter estimation based on passive acoustic techniques using the polynomial Wigner–Ville distribution, J. Acoust. Soc. Am. 102 (1) (July 1997) 207–223.