Interpolation techniques for non-stationary signal analysis

PIh S0263-2241 (97)00030-4

ELSEVIER

Measurement Vol. 21, Nos. 1- 2, pp. 1 16, 1997 © 1997 Elsevier Science Limited. All rights reserved Printed in The Netherlands 0263-2241,/97 $17.00 +0.00

Interpolation techniques for non-stationary signal analysis Gregorio Andria, Filippo Attivissimo, Annie Lanzolla Dipartimento di Elettrotecnica ed Elettronica, Politecnico di Bari, Bari, Italy

Abstract

This paper deals with the investigation of the measurement accuracy of characteristic parameters of digitized non-sinusoidal signals. Errors caused by frequency variation and lack of synchronization are also investigated. The correct expression of these errors in terms both of typical parameters of the used window and characteristics of the signal are given to allow a suitable choice of the filter and of its parameters, which may increase the estimation accuracy. The relations between the stationary case errors and the non-stationary ones are carried out. Finally, a new interpolation technique which allows very accurate measurements of the instantaneous frequencies and amplitudes of the signal, even if heavily non-stationary components are present, is proposed. © 1997 Elsevier Science Ltd.

Kevwords: Time frequency distribution; Interpolation; Accuracy characterization; Non-stationary signal analysis

1. Introduction

Non-stationary signals are very common in various fields of interest, such as theory of systems, control of machines, signal processing, measurements on power systems, oceanography, seismology, acoustics and so on [1]. As a consequence, more and more scientific fields require efficient and robust on-line measurements of a number of physical quantities, which are generally non-stationary. Moreover, measuring of time-evolution of quantities characterizing the dynamic behaviour of complex systems gives rise to an efficient real-time control of some important process state-variables. Then, all these measurement and control processes require higher and higher performance in terms of flexibility, accuracy and speed in estimating the time-varying quantities of interest. Non-stationary signals are generally multicomponent; they have time-varying amplitudes and frequencies, called 'instantaneous characteristic parameters', whose real-time estimation is now considered very important by the measuring engineers [2]. A considerable variety of methodologies has

already been developed to measure the typical parameters of a multifrequency signal, both in time and in frequency domain [3]. Generally, the second approach is more attractive when high computational efficiency and real-time constraints are required even though it can have low accuracy due to leakage effects [4]. With reference to the last problem, several techniques have been proposed with good results, so that sensible improvements have been obtained in measurement accuracy and in the achieved frequency resolution [5,6]. Unfortunately, although these methods applied well to the stationary case, they gave rise to poor results when applied to the non-stationary case. In this paper, a new method is proposed to correct the estimates, coming from the well-known short-time Fourier transform (STFT). It is shown that both the undesired effects of long-range and short-range leakage caused by instantaneous frequency variations can be minimized, so that a considerable increase of the measurement accuracy is obtained, with a further low computational effort. After developing some windowing and interpo-

2

G. Andria et al.

lation theories relevant to non-stationary digitized signals, simple expressions of errors have been obtained for classical windows and fiat-top ones [5], in terms of type of frequency variation and of the consequential fractional spectral bin of frequency deviation. Finally, it has been shown that these errors are an extension of those relevant to the stationary case. Then, it can be proved that the non-stationary ones are obtained by the sum of the relevant stationary errors and by suitable interpolation polynomials.

following form: s ( t ) = ~ Xk(t)=~ Ak(t)e j°k(t) k=0,1 ..... K k

(1)

k

where K is the maximum harmonic order, and

Ak(t) and Ok(t) are called 'instantaneous amplitude' and 'instantaneous argument' of the kth harmonic component, respectively. According to several authors [2,9,10] it is possible to define the 'instantaneous frequency' of the kth harmonic component as: A-

1 dOk(t) 2re

(2)

dt

2. Preliminary considerations

2.1. Stationary and non-stationary signals The classical methods of noisy signal analysis are based on the fundamental assumption that the time series of samples are stationary or can be reduced to stationarity by some simple transformations. However, in real applications the stationary sample represents a mathematical idealization which may be only approximately valid in some cases; the analysis using traditional methods generally produces good results when the examined signals are almost stationary [4,5,7]. On the contrary, the application of these techniques to nonstationary cases will produce biased results; different approaches can be adopted with the aim of also giving, in this case, accurate analysis results and describing how the spectral content of a signal is changing in time. Due to both the simple and powerful basic idea and the relevant low computational effort, the STFT, natural extension of the Fourier analysis of the stationary case, is probably the more used method for analysis of signals whose fundamental parameters vary. Unfortunately, in the parameters estimation some undesirable effects are present which can affect the performance of this method [81. Now, let us consider a band-limited multifrequency analytic signal [9], whose time-varying amplitude and frequency can be represented as a linear combination of complex exponentials in the

The last expression represents the attempt to define a spectrum in the non-stationary case which may be considered as a 'local' power-frequency distribution. In this case, the whole time series cannot be expressed by means of sine and cosine functions with the associated conventional definitions of frequency and spectral component: then, a more general definition is necessary. Equation (2), which embodies the notion of instantaneous frequency as the derivative of the instantaneous argument, requires that the signal has to be replaced by its analytic version in the case of real values [11 ]. It is important to emphasize that this identification is consistent for a wide variety of signals but in some cases may produce erroneous results [9]. The quantity set ek(t)={fk(t); Ak(t)} is here considered as the set of the instantaneous 'characteristic parameters' of the multifrequency signal. In the cases in which Pk(t) do not change against time, i.e. fk(t)=fk;

Ak(t)=Ak

Vt>0

(3)

the multifrequency signal is considered here as stationary. If Eq. (3) is not verified, the signal is non-stationary. In a previous paper [7], the authors developed a new measurement methodology, characterized by a high level of accuracy, if applied to locally stationary signals. With the aim of extending the results to the general case of signals with frequency, however variable, the following expression of the

G. Andria et al.

harmonic instantaneous argument is considered: B

(4)

Ok(t)= ~ ebkt b b=0

by which for B = 0, 1, 2, 3 a continuous, stationary, linear FM and quadratic FM signal is given, respectively. Figure l(a) and (b) shows typical examples of variation laws of Ok(t) and fk(t) for different values of B. For the sake of analysis simplicity, in this work the errors relevant only to the linear variation of frequency will be examined; the extension to other kinds of variation requires only a small further computational effort. 2.2. Analysis of non-stationary signals via shorttime Fourier transform

The STFT of Eq. (1) is defined as [2]: S~(t,f)=

f

+oo

s(?)w*(7-t)e -j2~fr d7

(5)

-oo

where w(') is a suitable window function and the symbol '*' denotes the complex conjugate. As is well known, the STFT represents the instantaneous--or 'local'--spectrum of the signal centred around t. It is possible to consider its discrete version by sampling the signal at a frequency rate F and using a function window of finite length L according to: L-1

S'[(m,l)= ~ s ( n ) w * ( n - m ) e -jz~"uL

(6)

n=0

This time-frequency distribution is computationally efficient, but gives erroneous estimates of the characteristic parameters of the signal, owing to the finite length of the measurement interval, especially for heavily non-stationary signals, whose frequency components may change continuously [9]. Then, the periodic expansion of the digitized finite part of the signal exhibits discontinuities at the boundaries of the sampled time window of time length LAt. This phenomenon brings about negative effects in the spectra returned from the STFT-operation; they consist of spectral spreading in all the frequency domains, called 'long-range

3

leakage' which can heavily affect the analysis [4,5,12-14]. Figure 2(a) and (b) shows a typical digitized signal having two frequency components and its relevant leakage-affected power spectrum, respectively. On the basis of several previous works about the analysis of stationary and non-stationary signals, the authors can assert that an accurate evaluation of the set can be carried out if: (i) a suitable windowing operation is performed to minimize the spreading of each harmonic component and, consequentially, to reduce the longrange leakage effects; (ii) an efficient interpolation algorithm to remove the residual effects of the 'short-range leakage' is run. As is known [5,6, 12,13], the shortrange leakage consists of the residual spectral spreading of each frequency component fk(t) of the windowed signal. This phenomenon is not negligible in the range of frequency relevant to the spectral mainlobe of the used window, centred on fk(t), because of the generally fast asymptotic decay of the sidelobe of the window spectrum. Figure 3(a) and (b) shows the Hanning-windowed signal of Fig. 2(a) and the residual effect of 'shortrange leakage' on its relevant power spectrum, respectively.

3. Error analysis from STFT

Let Xk(t), the kth component of the multifrequency signal given in Eq. (1), having a constant amplitude Ak and an instantaneous argument given by: 0 k (t)

= C2kt z + Clk t + Cok

(7)

In this case, the relevant instantaneous frequency is, from Eq. (2): fk --

1 d0k(t) - =azkt+alk 2n dt

(8)

where a2k=C2k/n and alk=Clk/2rc represent the frequency slope and the initial frequency, respectively. The L-length discrete version of Xk(t),

4

G. Andria et al.

60( J J J /J J

L)

40C

/ J t f

B=3

/

f

/

B=2

O /

20C

1~

i

J

I

i

L

3

I

i

5 time (s)

(a)

..............

• • . . . , . . . . • • " "

r

i

_

7

9

xl0 4 2 1.6 / f

E

B=3

1.2

B=2

J

e

f

o.8 J J

j

~

... •.•

...

r

r

~:5...~. ....

0.4,

B=I

...

f f

_ ....

~

r ~

~ J

jr

r

0

2

4

(b)

time (s)

6

8

10

Fig. 1. Examples of variation laws for (a) the instantaneous argument, and (b) the instantaneous frequency, relevant to a nonstationary single-tone signal, for some values of B.

sampled at frequency rate F, according to the Shannon rule [4-6], is then obtained as:

xk(n)=Ake~"n/r(~-n+2atk)

(9)

by supposing null as the initial phase, without loss of generality. For the sake of simplicity, in Eq. (9) the time resolution At--1IF is omitted. As it is possible to show, the instantaneous hth spectral line of the relevant STFT is given by:

G. Andria et al.

5

150 t

~

50

'

I

O e.-

~e.-,

-50

-150 0

20

40

60

time (s)

(a)

100

80

> 6C e~

E

I

40

20

otTT

200

600

1000 frequency (Hz)

(b)

1400

1800

Fig. 2. Example relevant to an analytic non-stationary signal having two time-varying frequency components (its fundamental instantaneous frequency is linearly variable from 200 to 235 Hz in the observation interval): (a) time representation of its real part, and (b) relevant long-range leakage-affected power spectrum.

L-1

Xk,=(h)= ~ xk(n)e -j2(~/L)h" n=O

Ak L-1 l...a

L

n=o

t?jnn[(a2k/F2)n + (2alk/F) -- (2h/L)l

(10)

where 'ns' stands for n o n - s t a t i o n a r y case. As is known, Xkns(h) is relevant to the frequency hal, where A f = F/L is the frequency resolution. Let the hokth spectral line to be returned by S T F T operation be a line o f relative m a x i m u m in

6

G. Andria et al.

2 0 0

,

,

,

,

~

b

100

0

-10~

-201

' 10

' 30

(a)

' 50

70

time (s) 100

80 ;> 60

40

20

0;,,* l ,

200

(b)

-

v

.

.

.

.

.

600

.

.

.

1000 frequency (Hz)

1400

1800

Fig. 3. Example relevant to the analytic nonstationary signal of Fig. 2, windowed by the Hanning window (a) time representation of its real part, and (b) relevant short-range leakage-affectedpower spectrum.

the local spectrum relevant to the observation interval To = L A t . Assume now that the instantaneous frequency fk and consequently the m a x i m u m concentration of instantaneous energy of the kth harmonic component o f signal is located at the

relevant mean frequency ( f k } of the considered interval To. This assumption has already been discussed and proved in the literature (see particularly the 'tutorial' section in Ref. [9]), and is valid for every analytic signal. In this work we consider

G. Andria et al.

only linearly FM analytic signals, for the sake of simplicity. Then, the main frequency can be computed simply as the ratio between the frequency variation in the observation interval and its time duration. For other kinds of frequency changes, it is always possible to approximate as linear the frequency change in each sufficiently narrow observation interval, and the suggested method can also be successfully applied in this last case. Then, we can write: a2k

(11)

fk = (hok + 6k)Af = (fk) = alk + - -

2Af

where 6k represents the fractional spectral bin of frequency deviation relevant to instantaneous frequency. After a simple substitution, one obtains from Eq. (ll): hok = ~ L2 + al__kkL 2~ ~ F

-

7

It is easy to show that it depends on the window length, the slope a2k and the instantaneous frequency displacement 6k. The argument of exponentials in Eq. (15) is the sum of two terms: the first one takes into account the frequency variation of the signal, the latter one depends on the frequency displacement fig. In Fig. 4 the systematic error of amplitude is shown for different values of slope frequency; one can see that it increases as the slope frequency does. In the case of stationary signals, the error curve is symmetrical with respect to 6k; for signals with linearly variable frequency, instead, the error curves are asymmetrical, presenting a point of discontinuity for 6k=_+0.5, dependent on the sign of the frequency slope. Moreover, for non-stationary signals, an error occurs for 6k = 0 tOO, which represents the inaccuracy component due to the frequency variation.

(12)

t~ k

Finally, by substituting Eq. (12) into Eq. (10), it is possible to obtain an estimation "~k of the amplitude Ak from the following expression: Ak L-1 Xk,~(hok) = ~ ~ e j~"u"2k/r2'~"-L~+~2o~/L)I= ftk

n=O

(13) Consequently, the systematic error of amplitude, defined as:

4. Windowing and analysis of errors

A considerable reduction of the error expressed by Eq. (15) can be achieved by a suitable choice of the time-domain function window; its effect is to bring data smoothly to zero at the boundary, so that the discontinuity is reduced. The filtered version of the time-sequence (Eq. (9)) relevant to the kth component of the signal is given by: Yk (n) = X k w(n) = A k eJ(nn/F)[(a2kF)/n + 2a, k]

ak- Ak

ek,s- - Ak

(14)

i=0

(16)

can be written as: where w(n) is the discrete version of a cosine-type window, defined as [14]:

ek,, = ek,, (aZk,fk) = 1 -- Yns(a2k,fk) 2 eJnnt(a2k/F2)(n--L)+(2;~k/L)]l 11- - ~11 L-1 n=0

=

1 ~1

- L I

n=O

Ifikl<0.5

eJ"n[(a2k/FZ)(n--L)+(21Oklsign(a2k)/L)]l 16kl =0.5

(15) in the absence both of windowing and of appreciable interference between adjacent spectral responses.

[ 2rcni "~ w(n)= ~ ( - 1 ) i b i c o s ~ - - - £ - - )

i=0

(17)

As determined in Eq. (17) G and bi are the order and the/th coefficient of the used window, respectively. Now, the STFT of Eq. (16) is the result of the frequency-domain convolution product of the frequency spectrum of signal by the window

8

G. Andria et al.

i 0.4 ens

--

0.3

........ i ' : , ~ ~ ~ . a 2 k =

750i

. /////S

i

........

0.2

0.1

0 -0.4

-0.2

0

0.2

u~k

0.4

Fig. 4. Systematic error of amplitude for 32-length Boxcar window vs 6k for different positive slope values of linear frequency variation (-) and for stationary case (- • -).

spectrum W(h)= ~

(-1)i~[D(h-i)+D(h+i)l

(18)

i=1

where D(h) represents the Dirichelet kernel [14], according to:

ek.==ekn=(aZk,6k)= 1-- Tw.,(aZk,6k) bi = 1 - - ~ ( - - 1 ) i 2[Xkn,(6k+i)+Xk.s(~k--i)] i=0

1L-1

Yk,,(h)=Xk.=(h)s W(h)= ~ Y' yk(n)e E-j(2'~h"/L)] n=0

(19) It is possible to show that the instantaneous frequency estimation is still given by Eq. (1 I); after some algebra, one obtains the estimation of amplitude Akw for cosine-windowed signals as (see Appendix A):

~ (_l)ib,[~.,(hok_i) Yk"'(hok)= Ak 2 i=0 + Xk..(hok+,)] = .'~kw

values of the characteristic parameters of the analyzed signal:

(20)

where )(k..(h) is the normalized STFT of the signal filtered with the boxcar window. Finally, one can determine the following expression for the systematic error in amplitude measurement, in terms of the used cosine window and of the instantaneous

(21) In this case, the relevant error curves are practically symmetrical with respect to 6k, presenting an error even if the signal is synchronized. The windowing operation by flat-top windows, instead, gives rise to error curves characterized by a marked asymmetry with respect to 6k (Fig. 5). It is easy to show that for the boxcar window 7~==y.=, that Eq. (15) is a particular expression of Eq. (21). Moreover, for stationary or quasistationary signals, 7w~=approaches the scalloping loss [14] of the used windows and, consequently, ek,= grows into the error for stationary signals ek [14], relevant to a single tone, with a frequency equal to the mean value of the considered band. In Fig. 6 the amplitude errors for three different windows are shown, in the case of a non-stationary

G. Andria et al.

p

9

p

p

n

p

ens

-o.o4 t i -0.4

a2k = 1500

-0.2

!

!

0

0.2

~k

0.4

Fig. 5. Systematic error of amplitude for 32-length 3T-MS flat-top w i n d o w vs 6 k for different slope values of linear frequency variation ( ) and for stationary case (- • ).

0.4 ens x

•

x .

l" "x

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

'

.

.

.

.

.

.

t~

x

02 •

-

t

. . . . . .- . . . . . . . . . . . . . . . . . . . x

z

"x

-0.5

0

8k

Fig. 6. C o m p a r i s o n of systematic error of amplitude for Boxcar ( - - -), H a n n i n g (--) and 3T-MS flat-top ( values of 6, with a2k = 500 Hz s 1 and L = 32.

signal; as in the stationary case, the flat-top windows produce a negligible error for each value of 6k. Unfortunately, by using this kind of filter, it is

0.5 ) windows, for different

impossible to obtain an accurate estimation of the instantaneous frequency [ 15]. Finally, it is necessary to note that, even if the

10

G. Andria et al.

relationship Eq. (21) has been carried out for cosine windows, it is possible to extend this methodology to every kind of filter.

o f a2k~

di (a2k) ~--cl(a2k) = di2 a2k +dil a2k + dio i = 0,... M (25)

5. Stationary-time varying error relation In the previous section, the expression of the systematic error of amplitude in Eq. (21) is given in terms of the kind of frequency variation and of the frequency displacement 6k, when a cosine window is used. Now, it will be shown that it is possible to express ek,s as a sum of the amplitude error eke, relevant to a stationary signal, and of a suitable polynomial approximation pu(a2k,bk), which accounts for the negative effects of the timevarying parameters of the signal. By examining the curves of the amplitude error it can be noted that, for each used window, the decrease of slope a2k produces a reduction of the systematic error as far as necessary to be superimposed to the error relevant to the stationary case (a2k=0). Starting from this result, the possibility of separating the negative effects due to the frequency variation from the other ones caused by the asynchronous sampling of the signal, supposed stationary, has been investigated. Then, we can write: ek,, ( a2k,6k ) ~- eke(f k) + p M( a2k ,t~k)

(22)

where: eks = eks (6k) = 1 -- ?ws (6k) =l_6sin(rc6k) ~ (--1) i b____~/ bon i=o 62-i 2

(23)

according to the expression formulated in Refs [ 14,15]. The function ?w~(6k)represents the scalloping loss of the used window, and PM is a polynomial function of M-order, given by: M

PM(a2k,6k)= ~ di(a2k)6~ i=0

Both Eqs (23) and (24) are dependent on the used window. In Figs. 7 and 8 the coefficients di of the polynomial function PM obtained by Eq. (22) against the slope frequency are represented by solid lines for 32-length Hanning and Nuttall 3T-FD windows, respectively; the coefficients di obtained by Eq. (25) are represented by stars. It has been ascertained that the maximum absolute error relevant to the polynomial approximation Eq. (25) is less than 1 x 10 -4, for all the used windows and for every value of L. The authors determined that for class I-windows [16] the order of polynomial approximation is M = 2 , with dl(a2k)=O (Fig. 7); in this case the error of amplitude is simply given by: ek.,(a2k,6k)'~eks(6k)q-d2(a2k)6

2 + d0(aZk )

(26)

for flat-top windows [17], otherwise the order is M = 3 , with d~l¢ 1 (Fig. 8), so that one has: ek°. (a2k,6k) ~--eks(6k) + d3(a2k)63 + d2(a2k)62 +dl (a2k)6k + do(a2k)

(27)

Naturally, the polynomial approximation allows one to obtain for each window a set of polynomial coefficients d~j depending on the number of samples. The coefficients relevant to Hanning window are listed in Table 1 for some usual values of L. In Fig. 9, the amplitude errors obtained by using the polynomial approach are shown for three different windows. The very good agreement between simulation results and the approximation analysis clearly verifies the suitability of the proposed interpolation method; in fact the maximum absolute error is less than 8 x 10 -16. It is clear that the error of amplitude, expressed by Eq. (27), is univocally determined when an estimation of both 6k and a2k values is carried out.

(24)

where the terms d~ are obtained by finding the coefficients of a polynomial that fits the difference ekns(a2k,Ok)-eks(•k) according to the least-squares method; these ones can be expressed as a function

6. The interpolation technique In the previous section, it is shown that the use of class I-windowing filters produces better results

G. Andria et al.

11

0.15

0.05

-0.05 Z

500

i

i

1000

1500

a 2k

2000

Fig. 7. Representation of theoretical interpolation coefficients ( ) vs a2k for 32-length Hanning window. The stars represent the simulation results.

0.1

I

! d3

......

~

.........

I

x

~:

!-,.di . . . . . . . . . . . . .

-0.1

-0.2

500

1000

1500

a 2k

2000

Fig. 8. Representation of theoretical interpolation coefficients (-) vs for 32-length 3T-MS flat-top window. The stars represent the simulation results.

in the amplitude estimation, even though a residual component of error, depending on the frequency displacement 6k and the frequency slope a2k, is present. To obtain an accurate evaluation of characteristic parameters of the signal, a suitable inter-

polation technique is necessary, which, by starting from the frequency estimation, produces a good estimation of the instantaneous amplitude and phase of the considered component. It must be emphasized that this procedure is consistent,

G. Andria et al.

12

Table 1 Polynomial coefficients d o for the Hanning window for some values of L

do2 dm doo

dzz d21 d2o

L = 16

L = 32

L = 64

L = 128

L =256

- 1.2E- 9 - 1.0E-8 3.2E-6 2.2E-9 2.4E-8 - 7.4E- 6

- 1.5E- 8 -4.9E-6 1.9E- 3 2.8E-8 1.1E-5 --4.3E- 3

3.7E- 8 -3.7E-4 1.1E- 1 -2.0E-7 8.2E-4 -2.4E- 1

1.7E- 7 -5.7E-4 2.1E- l -4.8E-8 1.2E-4 5.3E- 1

9.0E- 8 -2.9E-4 -6.0E- 1 - 1.4E-9 3.5E-6 6.0E- 1

0.4 e

ns

i

0.3

0.2

0.1~

-0.4

-0.2

0

0.2

,_,Xk

0.4

Fig. 9. Simulation results (*) carried out to validate the theoretical values (-) of the polynomial approximation of systematic error by Boxcar, Hanning and 3T-FD Nuttall windows, for different values 6k, with a2k = 500 Hz s 1 and L = 32.

because the correct estimation of instantaneous frequency assures also correct amplitude and phase estimates [4]. If the kth component of the signal is located at the bin hog in the spectrum returned from STFT of the signal, it is necessary to determine the relevant fractional bin 6k. Then, one has: f ok = (hok + &k)Af

(28)

Now, by substituting Eqs(21) and (23) into Eq. (26) and after some algebra, it follows: ~ w , , ( a 2 k , 6 k ) ~ w s ( t ~ k ) - [ - d 2 ( a Z k ) 6 2 -[-do(azk )

(29)

where the first term relevant to the second member accounts for the effect of variation of ~k and the

others account for the effect of variation of frequency. In analogy with the interpolation algorithms proposed for stationary signals [5], it should be pointed out that, for each value of frequency slope, the correct evaluation of the fractional bin can be carried out by considering the spectral line of relative maximum in the spectrum Yk.s(azk,hok) and its adjacent spectral lines Yk.,(a2k,hok+l) and Yk.s(a2k,hok--1 ). Let us define now the new quantity ctk =

Y k , s ( a 2 k , h o k ) - Yk,s(a2k,hok d-

1)

Yk.s(a2k,hok) -- Yk.s(a2k,hok -- 1 )

(30)

G. Andria et al.

13

1.2 fi k =0.5

~k

~ k = .0.4i

0.8

. :

. . . .

8k._ 0. 3 . . . . . . .

.......

0.4

. . . . . . . . . . .

~k.=O, 2

.~k =O'li

1000

2000

a2k

3000

Fig. 10. Behaviourof quantity ctk vs the frequencyslope for differentvalues of @

that can be easily measured from a peak detector algorithm, applied to the local spectrum returned from the STFT. It has been determined that the largest spectral line generally presents an attenuation with respect to the relevant stationary case, but the ratio Eq. (30) is almost constant, for all the values of 6k in the range [ - 0 . 5 , +0.5]. Consequently, the value of ~k can be given by:

sured quantity ak and the unknown quantity ~k: ek =

Yk~ ( hok ) -- Yk~ ( hok + 1 ) Y~, ( hok ) - Yks ( hok -- 1 )

(31)

where Yk.~(hok), Yks(hok + 1) and Yk~(hok-- 1 ) represent the spectral lines relevant to the stationary signal with constant frequency equal to the mean value of the considered time-slice. Many simulation results verified this assertion; this is an important result in determining 6k independently on the variable characteristics of the signal. Figure 10 shows how ek do not depend on a2k, for different values of 3k. NOW, by substituting Eq. (23) into (31), after some algebra (see Appendix B), we can write the following important relationship between the mea-

(32)

(26k -- 1 )[6k + ( G + 1 )]

Then, the estimated value of 6k can be obtained by inverting Eq. (32): - ( 1 + 2G)(ek + 1 ) _+[( 1 + 2G) 2 ~k =

~k -~

(26k + 1 )[6k -- (G + 1 )]

x (ctk + 1 )2 + 8(~k __ 1)2(G+ 1 )]1/2 4(~k -- 1) ~k¢l 0

(33)

~k= 1

where only the plus sign has to be considered as physically significant. Finally, the frequency and amplitude estimation relevant to the kth component of signal can be determined by Eqs (28) and (20), respectively. For example, if the Hanning window is used for the analysis, 6k is given by: [9(1 + (Zk)2 -~- 16(1 --ek) 2] --3(1 + .~Zk) 6k = ~ 4(~k ~ ~ 0

~k V~ 1 ~k = 1 (34)

14

G. Andria et al.

0.4

5k 0.2

0

-0.2

i!!!!!ii'!

-0.4 i

0

(a)

10

i

i

20

30

I

40

m

50

I

140

130

120

110

(b)

0

10

20

30

40

m

50

Fig. 11. Behaviour of (a) &k and (b) fk, respectively, vs normalized time ( ) . Stars represent the estimate values by proposed interpolation.

The error of amplitude expressed by Eq. (26) for class I-windows and by Eq. (27) for flat-top windows is univocally determined when the value of

(~k does not depend on a2k (being 0~k independent from a2k), the frequency slope can be obtained on-line by the estimation of a2k is known. Because

G. Andria et al.

15

the two values of frequency relevant two consecutive time slots: fk(t +At) --fk(t) a2k = (35) At

permit the user to obtain a measurement instrument suitable for accurate real-time analysis of non-stationary signals.

where T is the sample period. To verify the proposed technique, a test signal has been used, with Ak = 1 V and a frequency linearly variable with slope a2k = 500 Hz s-1 from an initial value alk = 100 Hz, filtered by a 32-length Hanning window and sampled at a frequency rate F = 800 Hz. Figure 11 (a) and (b) shows the behaviour of the estimation of 6k andfk versus normalized time m = t/T, respectively. The maximum absolute error relevant 6k derived from this procedure is less than 6 × 10-4; for 3T-FD and 4T-FD Nuttall windows the maximum is reduced to 8 × 10 -6 and 3 × 10 -5. The result was less than 0.01% for the percentage error relevant to the frequency estimation.

References

7. Final remarks and conclusions

An analysis of the systematic components of error relevant to non-stationary signal analysis via STFT is carried out and the effects of different windowing functions are examined. It is shown how fiat-top windows produce a good estimation of amplitude but a very poor instantaneous frequency estimation, so that they do not have to be considered. On the contrary, using the class I-windows reduces the effects of amplitude inaccuracy caused by frequency variation and by the inevitable finite length of processed samples. Thus, an original and efficient STFT-based algorithm has been proposed; it gives an accurate determination of the fractional spectral bin of frequency deviation and, consequently, of the instantaneous frequency and amplitude. Numerical simulations confirmed the performance of the suggested method, also in the case of heavily nonstationary signals. It is necessary to emphasize that the noticeable computational efficiency and accuracy of the proposed algorithm, added to the considerable development that involves the fields of electronic technology and digital signal processing, may

1. Loughlin, P. J., Special Issue on time-frequency analysis. IEEE Proceedings, 1996, 84(1), 1195. 2. Hlawatsch, F. and Boudreaux-Bartels, G. F., Linear and quadratic time-frequency signal representations, IEEE S P Magazine, 1996, 21. 3. Kayk, S. M. and Marple, S. L., Spectrum analysis A modern perspective. IEEE Proceedings, 1981, 69 (11 ), 1380. 4. Schoukens, J., Pintelon, R. and Van Hamme, H., The interpolated fast Fourier transform: A comparative study. IEEE, 1992, IM-41(2), 226. 5. Andria, G., Savino, M. and Trotta, A., Windows and interpolation algorithms to improve electrical measurement accuracy. IEEE Transactions on Instrumentation and Measurement, 1989, IM-38(4), 856. 6. Offelli, C. and Petri, D., Interpolation techniques for realtime multifrequency waveform analysis. IEEE Transactions on Instrumentation and Measurement, 1990, IM-39( 1), 106. 7. Andria, G., Attivissimo, F. and Savino, M., Instantaneous power measurement in time-frequency domain, E T E P 6, 1996, 5. 8. Portnoff, M. R., Time frequency representations of digital signals and systems based on short-time Fourier trasform. IEEE Transactions on Acoustics, Speech and Signal Processing, 1980, 48(2), 55. 9. Cohen, L., Time-frequency distributions-.A review. IEEE Proceedings, 1989, 77(7), 941.

10. Boashash, B., Estimating and interpreting the instantaneous frequency of a signal~Part 1: Fundamentals. IEEE Proceedings, 1992, 80(4), 520. 11. Nuttall, A. H., On the quadrature approximation to the Hilbert trasform of modulated signals. IEEE Proceedings, 1966, 54(1), 1458. 12. Jain, V. K., Collins, W. L. and Davis, D. C., High-accuracy analog measurements via interpolated FFT. IEEE Transactions on Instrumentation and Measurement, 1979, IM-28, 113. 13. Grandke, T., Interpolation algorithms for discrete Fourier transforms of weighted signals. IEEE Transactions on Instrumentation and Measurement, 1983, IM-32, 350. 14. Harris, F. J., On the use of windows for harmonic analysis with the discrete Fourier transform. IEEE Proceedings, 1978, 66(1), 51. 15. Andria, G., Savino, M. and Trotta, A., FFT-based algorithms oriented to measurements on multifrequency signals. Measurement, 1993, 12(1 ), 25. 16. Rife, D. C. and Vincent, G. A., Use of the discrete Fourier transform in the measurement of frequencies and levels of tones. Bell Syst. Tech. J., 1970, 49(2), 197. 17. Salvatore, L. and Trotta, A., Flat-top windows for PWM

G. Andria et al.

16

waveform processing via DFT. IEEE Proceedings, 1988, 135(6), 346.

tionary signal and its right adjacent spectral line; these can be expressed in terms of the fractional spectral bin of frequency deviation 3k and the window coefficients [5] as follows:

Appendix A Let yk(n) be the kth component of windowed sequence Eq.(16); the relevant STFT may be computed, according to Eq. (19), as:

Yk~(hok)-- Yk,(hok + 1)=3k sin(r~6k) ~, (-- 1)i boT~ i=o bi sin[~z(6 k + 1)] ×- + --(6 k + 1) ~k2 -- i2 bo

L-1

Yk,8(h) = ~

yk(n)e -jt2'~hn/L)

6

n=O

x ~" (--1)'

= ~

i= 0

eJ(~./F)~(.2k/F +2.101

bi (t~k +

1)2

(B1) -- i 2

For the sake of simplicity, the Hanning window is used; by considering the relationship between its coefficients ( b o = - b 0 , Eq. (B1) can be rewritten

n=O

as;

(A1) By using Eulero's formula, after simple algebra, one can write:

_(3k + 1) sin[/r(3k + 1)] [

biAk L- 1

Ykns(h) =

(-- 1)~ ~ i=0

62k+ 1

Yk (hok)-- Yk,(hok + l ) = 3 k sin(rc~k) ~ [

1 ((~k+ 1)2

1 1 (3k + 1)2

n=O

(B2)

>( { eJTtn[(a2kn/F2)+2(alk/F)-(2/L)(h +i)] ..~ cjnn[(a2kn/FZ)+2(alk/F) - (2/L)(h- i)] }

(A2)

Now, by fixing:

Yk,(hok)-- Yks(hok + 1)=

1 L-1 ) ( k " ~ ( h ) = L n~--O eJrtn[(a2~:/F2)n+(2axk/F)-(2h/L)]

sinOz6k)(23 k + 1) m~k(3k + 1)(6k + 1)(3g +2)

(A3)

(B3) Analogously, the difference between the maximum and its left adjacent spectral line is given by:

it is possible to obtain:

Yk.,(h)=&~

After some algebra, it is possible to write Eq. (B2) in compact form as:

(--1)i~[Xk,,(h+i)+Xk,,(h--i)] Yk~(hok)-- Yk,(hok + 1 ) =

i=0

(A4) Finally, the amplitude estimation is obtained by substituting into Eq. (A4) the frequency in hok given in Eq. (12).

Let Yks (hok), Yks (hok+ 1) be the spectral line relevant to maximum in the spectrum of the sta-

1)

=3k(fk -- 1)(6k + 1)(6k --2) (B4)

Finally, by substituting Eqs (B3) and (B4) into Eq. (31) one has: -

Appendix B

sin(ruSk) ( 2 3 k --

(23k + 1)(3k --2)

(26k -- 1)(3k +2)

(B5)

In the case of a generic class I-windows of G order, it is easy to conclude that Eq. (B5) will be expressed by Eq. (32).