Optimal Settings for Measuring Frequency ... - Semantic Scholar
3276
IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 58, NO. 9, SEPTEMBER 2009
Optimal Settings for Measuring Frequency Response Functions With Weighted Overlapped Segment Averaging Jérôme Antoni and Johan Schoukens, Fellow, IEEE
Abstract—This paper investigates the measurement errors involved in estimating frequency response functions (FRFs)—and related quantities such as the coherence function—from weighted overlapped segment averaging, a technique that has become a standard in modern data analyzers due to its computational advantages. Particular attention is paid to leakage errors, for which this technique has frequently been criticized. Our main result is that a half-sine or diff window with about 2/3 overlap achieves the best compromise to reduce leakage errors in the case of stationary random excitations, and this is independently of the system FRF. This conclusion is to be contrasted with the customary habit of using a Hanning window with 1/2 overlap. The same reasoning confirms that a rectangular window without overlap minimizes measurement errors in the special case of multisine excitations. Moreover, practical formulas are provided to the reader for computing the bias and variance of the frequency response and coherence functions in the general case. Index Terms—Frequency response function (FRF), optimal window, spectral leakage, transient effects, weighted overlapped segment averaging (WOSA).
I. I NTRODUCTION
T
HE MEASUREMENT of frequency response functions (FRFs) from experimental data is of continuous interest in numerous engineering fields and, more specifically, in any application where one is concerned with inferring the structure of an unknown linear system from its measured inputs and outputs. Although the parametric identification of the system is usually the ultimate objective—on which there exists a vast literature—the nonparametric estimation of FRFs offers a simple and yet very useful preliminary description of the system in the frequency domain. For instance, measured FRFs—or related quantities such as the coherence function—can be used to better understand a system order or to check its (non)linearity. Modern FRF measurement techniques are extensively based on the discrete Fourier transform (DFT), owing to its very efficient computation by means of the fast Fourier transform (FFT) algorithm. This flexibility, plus the constantly increasing capacity of computational devices, has led to the proposal of several estimators whose relative optimality strongly depends on the statistical nature of the excitation signals (i.e., inputs; see Manuscript received May 6, 2008; revised April 9, 2009. Current version published August 12, 2009. J. Antoni is with the Laboratory Roberval of Mechanics, University of Technology of Compiègne, 60205 Compiègne, France (e-mail: [email protected]). J. Schoukens is with the Department of Electrical Engineering, Vrije Universiteit Brussel, 1050 Brussels, Belgium. Digital Object Identifier 10.1109/TIM.2009.2022376
[1] for a comprehensive overview on the subject). Surely, the best possible choice is to use the maximum-likelihood estimator together with periodic excitation signals, since only in that case can the DFT be exactly used without finite-length effects. There may be, however, some technical or psychological reasons that lend the experimenter to prefer other types of excitations. Among those, a popular subclass is encompassed by stationary random signals (e.g., binary signal, random uniform, and random Gaussian). Different FRF estimators have been proposed over the years in this context, which essentially find their roots in experimental spectral analysis. They all rely on various frequency smoothing or averaging strategies, with the fundamental assumption (hope) that the estimator variability can so be stabilized without distorting too much the smooth structure of the FRF. In any case, this obviously leads to a difficult bias/variance tradeoff, which nonparametric FRF measurement techniques have always been criticized for. In this paper, we focus on the so-called weighted overlapped segment averaging (WOSA) or Welch’s procedure. The WOSA procedure consists of segmenting the available records into several short-length segments, computing their DFTs, and then averaging their products to get a stable estimator of the system FRF. Although the idea of the WOSA procedure is quite old [2], it has surprisingly received very little attention in the scientific literature. When it comes to the FRF measurement, its exact statistical performance and optimal settings remain virtually unknown. It is the objective of this communication to partially fill in these gaps. In particular, we address important and practical questions. 1) Which window shape should we use to taper the data? 2) How much can we gain from overlapping adjacent segments? 3) Which optimal amount of overlap should we set between adjacent segments? The answers to these questions will obviously provide the experimenter with some useful guidelines to make optimal usage of the WOSA technique in FRF measurements. This paper is organized as follows. In Section II the principle and limitations of the WOSA procedure are reminded. The statistical errors involved by the WOSA procedure are quantified in Section III, in terms of systematic and stochastic errors due to leakage and measurement noise. From these results, the optimal settings of the WOSA parameters are then deduced in Section IV. Section V shows that the conclusions remain unchanged for the WOSA estimation of other spectral
0018-9456/$26.00 © 2009 IEEE
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
ANTONI AND SCHOUKENS: OPTIMAL SETTINGS FOR MEASURING FRFs WITH WOSA
Fig. 1.
3277
Plant model.
quantities such as the coherence function. Section VI proves that the same reasoning can be applied to the special case of multisine excitations, although the optimal settings are then found to be very different. Finally, Section VII presents some numerical experiments. The proofs of our main results are sketched in the Appendix. II. FRF M EASUREMENT W ITH THE WOSA P ROCEDURE A. Classical Approach Let us consider a pair of input/output stationary random signals {u(n); y(n)}, n ∈ Z, and an assumed linear time-invariant relationship, i.e., y(n) =
∞
g(k)u(n − k) + νy (n)
Ywi (ωk ) =
where g(k) stands for the (causal) system impulse response, and νy (n) for the additive measurement noise at the output (see Fig. 1). The measurement of the system con FRF then −jωk g(k)e = sists of estimating the values of G(ω) = ∞ k=0 F{g(k)} at some specific frequencies ωk in a finite set Ω from L−1 the finite-length records {u(n)}L−1 n=0 and {y(n)}n=0 . The classical approach to this problem is to estimate G(ω) as the ratio of the smoothed cross-periodogram to the smoothed periodogram, i.e., ∗ M 2πi 2πi ω U ω g Y − − i L k k L i=−M L L ˆ k) = (2) G(ω M 2πi 2 | i=−M gi |UL ωk − L where L−1
y(n)e−jωk n = DFT{y(n)},
n=0
idea of the WOSA estimator (5) is to arrive at the same result by averaging over segments of data. Namely, let us denote
(1)
k=0
YL (ωk ) =
Fig. 2. Principle of the segmenting and windowing process. The signal of length L is divided into M (weighted) segments of length N that are overlapping by N − R samples.
ωk =
2πk L (3)
is the DFT of {y(n)}L−1 n=0 —and similarly for UL (ωk )—and where the gi ’s, i = −M, . . . , M , are the taps of a suitably designed low-pass filter. The FRF estimator (2) is directly inherited from the early works of Blackman and Tukey in power spectrum estimation [3] (it is also known as the Blackman ˆ k ) is the and Tukey procedure). Note that the so-defined G(ω weighted mean-square estimator of G(ωk ) given the set of (noisy) observations {YL (ωk − 2πi/L)}M i=−M and the set of re. The statistical performance gressors {UL (ωk − 2πi/L)}M i=−M of this approach has been thoroughly studied in references such as [4]–[6]. B. WOSA Procedure Whereas the Blackman and Tukey estimator (2) controls stability by smoothing over adjacent frequency lines, the basic
iR+N −1
w(n − iR)y(n)e−jωk n ,
n=iR
ωk =
2πk N (4)
−1 the short-time DFT of the segment of data {y(n)}iR+N ann=iR chored at time iR and tapered by the smoothed N -long window −1 {w(n)}N n=0 ; the WOSA FRF estimator is then obtained as the ratio of the averaged cross-periodogram to the averaged periodogram, i.e., M −1 Ywi (ωk )Uw∗ i (ωk ) ˆ G(ωk ) = i=0 (5) M −1 2 i=0 |Uwi (ωk )|
with M = (L − N )/R + 1 as the total number of data segments (x standing for the greatest integer less than or equal to x) and 1 − R/N as the fraction of overlap between adjacent segments (see Fig. 2). The WOSA procedure seems to have become a standard in commercial data analyzers in virtue of its simplicity and many computational advantages. 1) It can be very efficiently computed by means of several short FFTs (which is always faster than computing one very long FFT). 2) It requires few memory allocation, and it can actually be recursively computed in almost real time (this is the solution implemented in current commercial data analyzers). 3) It is robust against nonstationary, outlying, or missing data since the process of segmenting allows for easy detrending and removal of corrupted segments. The performance of WOSA has been investigated in some early papers within the context of spectrum estimation. In a series of papers [7]–[9], Nutall showed that the WOSA procedure can be made nearly as efficient (in the sense of achieving the same statistical stability given a similar frequency resolution) as the Blackman and Tukey procedure provided that sufficient
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
3278
IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 58, NO. 9, SEPTEMBER 2009
Fig. 3. Transient effects with the (left) rectangular window and (right) Hanning window.
overlap between adjacent segments is applied. From numerical simulations, Nutall reported that ideally 62.5% overlap should be imposed, but that 50% is more practical and almost just as efficient. This recommendation seems to have become the standard nowadays. In another series of papers, the key role of tapering the segments with a smooth data window (as opposed to the rectangular window) was clearly demonstrated (see [10] and the references therein). Nevertheless, there has been no unique answer so far to the question of an optimal window shape; instead, many “good” windows have been proposed over the years, each with its own peculiarities [11]–[15]. Surprisingly, there are very few related results when it comes to applying WOSA for the measurement of FRFs [16], [17]. Indeed, this issue happens to be more involved than that of spectral analysis. One major peculiarity of the WOSA estimator when used for FRF measurements is that the use of several segments of reduced length makes it prone to severe leakage errors. If the mechanism of leakage due to Fourier transforming a truncated time series is well known in spectral analysis, its implication in the WOSA technique is more subtle and requires some explanations. It stems from the inherent but erroneous assumption that the tapered output segments are identical to the system responses of the tapered input segments, despite of the nonlinearity of the windowing process. This produces transient artifacts in the signals at the onset and end of each segment, which will perturb the FRF estimation process. As illustrated in Fig. 3, a typical solution to reduce the transient magnitude is to use a smooth data window—e.g., the Hanning window as opposed to the rectangular window—but this is generally at the price of an increase in the transient duration so that some compromise has to be found. Another intuitive solution is to increase overlap between adjacent segments so that transient artifacts got averaged out to a further extent. We now address such questions as whether there is an optimal window shape to taper the data with, as well as how much can be gained from overlapping adjacent segments. The answers to these questions are nontrivial, and when it comes to FRF measurements, the approach to obtain them is quite different from that used in classical spectral analysis. The proposed approach proceeds by first analyzing the measurement errors involved in the WOSA procedure. These are shown to stem from three sources, such as the following: 1) systematic errors; 2) stochastic errors due to spectral leakage; and 3) stochastic errors due to additive measurement noise.
This can be summarized by the following remarkably simple formula (see the proof in the Appendix) ˆ k ) = G(ωk ) + βL (ωk ) + βN (ωk ) G(ω
(6)
where βL (ωk ) is a nonzero-mean stochastic error accounting for spectral leakage, and βN (ωk ) is a zero-mean stochastic error accounting for measurement noise. We then pay particular attention to the leakage error βL (ωk ), the statistical analysis of which provides closed-form (asymptotic) expressions for the bias and the variance of the WOSA FRF as a function of the data-window shape and length, the number of segments, and the amount of overlap between adjacent segments. Incidentally, (6) proves that leakage affects both the bias and ˆ k ) since the variability of G(ω ˆ k ) = G(ωk ) + E {βL (ωk )} (7) E G(ω ˆ k ) = Var {βL (ωk )} + Var {βN (ωk )} . (8) Var G(ω It seems that only the deterministic nature of leakage (i.e., the first term on the right-hand side of (7) has been recognized in classical references on the subject) and that its random contribution has been ignored from variance analyses [2], [5]. However, there are many instances where the latter dominates the former [17]. III. S TATISTICAL E RRORS A. Assumptions and Notations In this section we assume, for simplicity, that A1) the input signal u(n) is real, zero mean, and stationary; A2) its correlation length is much smaller than the system time constant1 τc =
∞ 1 2 2 k=0 |g(k)|k ∞ k=0 |g(k)|
(9)
as expected from a well-designed input signal [1]; A3) the measurement noise is independent of u(n); 1 Note that this definition of the time-constant encompasses both notions of pure delay and of group delay of the impulse response of the system.
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
ANTONI AND SCHOUKENS: OPTIMAL SETTINGS FOR MEASURING FRFs WITH WOSA
3279
TABLE I CORRECTION FACTORS FOR THE HANNING WINDOW
TABLE II CORRECTION FACTORS FOR THE HALF-SINE WINDOW Fig. 4. (Left) Hanning and (right) half-sine windows v(τ ) and their corresponding (normalized) autocorrelation functions κ(τ ).
A4) the system time constant τc is much smaller than the data-window length (as classically recommended with WOSA and more generally with any nonparametric measurement of FRFs). Furthermore, to propose concise expressions, let us introduce the continuous window shape v(τ ) defined over the support set [0; 1[. The discrete N -long data window w(n) used in WOSA [e.g., in (5)] is obtained by sampling v(τ ) with N samples per unit length, i.e., n = 0, . . . , N − 1.
w(n) = v(n/N ),
of segments, and the window shape v(τ ). Note that Cv2 (θ) exists only for windows having a bounded derivative, a requirement that precludes, in particular, the rectangular window.
(10) B. Systematic Error
Let us also denote 1−|τ
|
v (|τ | + x) v(x)dx
R2v (τ ) =
(11)
0
Under Assumptions A1–A4 and provided that v(τ ) has a ˆ k )} of (7) bounded derivative, the systematic error Bias{G(ω can be shown to have the asymptotic expression2 (see the Appendix)
the autocorrelation function of v(τ ), and
E {βL (ωk )} −−−−→ N τc
R2v (τ ) κ(τ ) = R2v (0)
(12)
its normalized version (see Fig. 4). Finally, let 100(1 − θ) be the percentage of overlap between two adjacent segments, where θ = R/N (see Fig. 2). With these notations, we are now in position to define the two correction factors Cv2 (θ)
M (θ)−1
= − κ (0) − 2
i=1
i 1− M(θ)
κ (iθ)κ(iθ) (13)
Dv2 (θ) = 1 + 2
M (θ)−1
i=1
1−
i M (θ)
κ(iθ)2
(14)
which are explicit functions of the percentage of overlap, the number M (θ) = (L − N )/N θ + 1
(15)
Cv2 (1) G (ωk ) ∼ O(N −2 ) 2N 2
(16)
where the correction factor Cv2 (1) = −κ (0) is as given in (13) with θ = 1 and may be interpreted as the curvature of the autocorrelation function of the window at lag 0—see the first row in Tables I and II for some particular values—and G (ωk ) is the second derivative of G(ωk ). Expression (16) is instructive in two instances: First, since the bias is an O(N −2 ), long data windows should be used. Second, since it is proportional to −κ (0), data windows implying minimum curvatures should be used. This issue will be addressed in the next section. C. Stochastic Errors From (8), the total variance of the WOSA FRF estimator decomposes into one contribution due to leakage and one contribution due to measurement noise. We now provide the general expressions for these two contributions given an arbitrary 2 The
notation larger than τc .”
−−−−−−−−−→(·) N τc
means “converges to (·) as N becomes much
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
3280
IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 58, NO. 9, SEPTEMBER 2009
percentage of overlap and an arbitrary window shape. The sketch of the proofs is provided in the Appendix. 1) Stochastic Errors Due to Leakage: Under Assumptions A1–A4 and provided that v(τ ) has a bounded ˆ k ) due to leakage is (see the derivative, the variance of G(ω Appendix) Var {βL (ωk )} −−−−→ N τc
Cv2 (θ) 2 |G (ωk )| , ωk = 0 mod (π) M(θ)N 2 ∼ θCv2 (θ) · O(L−1 N −1 )
(17)
with the correction factor Cv2 (θ) as given in (13). Equation (17) can be somehow simplified when particular cases are considered: we deemed it useful to provide the reader with the equivalent formula reported in Tables I and II for the typical cases of a Hanning and a Half-sine window, respectively, with 0, 1/2, and 2/3 overlap. Since the variance is shown to be an O(L−1 N −1 ), one solution to reduce stochastic errors due to leakage (given a fixed signal length L) is to use long data windows. Another solution is to set the fraction of overlap θ = R/N that minimizes Cv2 (θ). This issue will be addressed in the next section. 2) Stochastic Errors Due to Measurement Noise: Under ˆ k ) due to measureAssumptions A1–A3, the variance of G(ω ment noise is (see the Appendix) Var {βN (ωk )} ≈
Dv2 (θ) · NSRout (ωk ) M (θ) ∼ θDv2 (θ) · O(L−1 N )
(18)
with the correction factor Dv2 (θ) as given in (14), and where NSRout (ωk ) = S2νy (ωk )/S2u (ωk ) stands for the output noiseto-signal ratio (NSR) at frequency ωk . Again, the reader may find it convenient to use the equivalent formulas reported in Tables I and II for the typical cases of a Hanning and a halfsine window, respectively, with 0, 1/2, and 2/3 overlap. Expression (18) is similar to that reported in [2]. However, we emphasize again the fact that it cannot alone explain the variability of the WOSA FRF estimator. IV. O PTIMAL S ETTINGS OF THE WOSA P ARAMETERS We now make use of the preceding results (16)–(18) to find the optimal parameters to be used with the WOSA procedure. Our first result concerns the optimal window shape that is found to minimize the errors due to leakage (both systematic and stochastic). Our second result specifies how much overlap should be used to minimize the overall stochastic errors (both leakage and measurement noise). A. Optimal Window Shape The optimal window shape should minimize the following mean-square error (MSE): MSE(ωk ) = |E {βL (ωk )}|2 +Var {βL (ωk )}+Var {βN (ωk )} . For the sake of simplicity and to avoid having the optimum depending on the unknown signal-to-noise ratio (SNR) and
Fig. 5. Cv2 (θ)/M(θ) as a function of the amount of overlap (1 − θ) for the (dotted line) half-sine and (continuous line) Hanning windows, with L = 10 000 and N = 100. Note that M (θ) = (L − N )/N θ + 1.
actual FRF, we consider only the minimization of the errors due to leakage (i.e., the first two terms on the right-hand side of the preceding equation). This amounts to finding the optimal window in the noise-free case. Turning back to (16) and (17), it is shown that, in the particular case of no overlap (θ = 1), such a window should minimize the curvature −κ (0) of its autocorrelation function. In the family of real symmetric windows, it can be shown (e.g., by standard variational analysis) that such a minimum is achieved by the half-sine window shape (see Fig. 4), i.e., v(τ )−−−−→ sin(πτ ), N τc
0 ≤ τ ≤ 1.
(19)
The half-sine window has several appealing properties. 1) It has the same minimum-bias performance as the recently proposed (complex and asymmetric) “diff” window [17]. 2) Its optimality is asymptotically (when N τc ) independent of the system FRF. 3) It is incidentally very similar to the parabolic window historically used in Welch’s original paper [2]. Moreover, it can be shown that the optimality of the half-sine window still approximatively holds when overlap is increased (θ < 1). This is because the value of Var{βL (ωk )} becomes almost independent of the window shape as θ → 0, as shown in Fig. 5. In fact, the superiority of the half-sine window over the Hanning window can be checked in Fig. 4, wherein the larger curvature of κ(τ ) at τ = 0 implies a smaller bias according to (16), and in Fig. 5, wherein the smaller correction factor Cv2 (θ) implies a smaller variance according to (17). Precisely, Tables I and II) indicate a reduction of 1.33 on the bias and, more modestly, of 1.33, 1.39, and 1.09 on the variance with 0, 1/2, and 2/3 overlap respectively. This superiority may be
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
ANTONI AND SCHOUKENS: OPTIMAL SETTINGS FOR MEASURING FRFs WITH WOSA
3281
Moreover, Figs. 5 and 6 show that not much more can be gained by forcing adjacent segments to overlap more than about 2/3 and 1/2, respectively, because the added information then becomes extremely redundant. This means that, depending on the window shape and whether leakage noise or measurement noise predominates, the optimal fraction of overlap is somewhere between 1/2 and 2/3. Since the actual SNR will rarely be known, a conservative choice is about 2/3 overlap, a limit above which the extra computational cost involved by forcing more overlap is surely not justified. More accurate settings may be formulated on a case-by-case basis, depending on both the used window and the leakage- or noise-free case, but we believe that the simple and conservative “2/3 rule” is just good enough in practice—it can actually be shown to conservatively hold for all smooth (i.e., with a bounded derivative) data windows, just as observed here above on the Hanning and half-sine windows.3 V. E XTENSION TO O THER S PECTRAL F UNCTIONS Fig. 6. Dv2 (θ)/M(θ) as a function of the amount of overlap (1 − θ) for the (dotted line) half-sine, (continuous line) Hanning, and (gray line) rectangular windows, with L = 10 000 and N = 100. Note that M (θ) = (L − N )/N θ + 1.
explained by the fact that, in the frequency domain, the half-sine window achieves a better balance between local leakage (width of the main lobe) and global leakage (roll-off of the secondary lobes). Finally, it should be emphasized that the half-sine window is no longer optimal—in the sense of minimizing the overall MSE—in the presence of (significant) measurement noise. It can be seen from (14) and (18) that the reduction of Var{βN (ωk )} requires a window shape with small inertia (a Dirac impulse in the limit), which is conflicting with the reduction of the leakage errors. Therefore, the optimum window will be a compromise in the general case, depending on both the unknown SNR and FRF. Addressing this difficult issue is outside the scope of this paper and is also unlikely to provide the end user with any simple recommendation. B. Optimal Amount of Overlap Having set the window shape that minimizes the effect of leakage errors, it is still possible to significantly reduce the MSE by acting on the overlap. As mentioned in the preceding section, the intuition is that overlapping will average out stochastic errors to a further extent by providing a larger number M (θ) of segments. Indeed, turning back to (13) and (14) and their illustration in Figs. 5 and 6, it is shown that Cv2 (θ)/M(θ) and Dv2 (θ)/M(θ) are globally decreasing functions of the amount of overlap (1 − θ). More precisely, Tables I and II show that from 0 to 2/3 overlap, the number of windows is roughly tripled, whereas the correction factors Cv2 (θ) and Dv2 (θ) are only slightly changed. For instance, with a Hanning window, this yields a reduction by 3.8 for the variance due to leakage and by 2 for the variance due to measurement noise. With a half-sine window, the reduction is by 3 for the variance due to leakage and by 1.7 for the variance due to measurement noise.
It is a natural question to investigate whether the optimal settings found hitherto for measuring FRFs also apply to the measurement of other spectral functions. Of particular interest is the WOSA estimator of the (squared-magnitude) coherence function and the alternative WOSA measurement of FRFs—the so-called “H2 ” function. A. Coherence Function The coherence function 2 M −1 i=0 Ywi (ωk )Uw∗ i (ωk ) 2 γˆ (ωk ) = M −1 2 M −1 2 i=0 |Ywi (ωk )| i=0 |Uwi (ωk )|
(20)
is a useful quantity for checking the linearity of the system under study. For the sake of simplicity, we only address here the impact of leakage noise on its measurement (i.e., the noisefree case), which, to our knowledge, has never been investigated before in the literature. Investigation of the combined effect of leakage and measurement noise is left for future research. Under Assumptions A1–A4 and provided that v(τ ) has a bounded derivative, the bias and variance of γˆ 2 (ωk ) due to leakage are (see the Appendix)
C 2 (1) |G (ωk )|2 Bias γˆ 2 (ωk ) −−−−→ v 2 N τc 2N |G(ωk )|2
E 2 (θ) |G (ωk )|4 Var γˆ 2 (ωk ) −−−−→ v , N τc M (θ)N 4 |G(ωk )|4
(21) ωk = 0 mod (π)
∼ θEv2 (θ) · O(L−1 N −3 )
(22)
where Ev2 (θ) = κ (0)2 + 2
M (θ)−1
i=1
1−
i M (θ)
κ (iθ)2 .
(23)
3 The formal proof of this assertion is only of technical interest and will be published elsewhere; it can actually easily be verified by computing (13) and (14) for different window candidates.
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
3282
IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 58, NO. 9, SEPTEMBER 2009
where Ak are real deterministic amplitudes, φk are independent random phases uniformly distributed over [0; 2π[ (whose outcomes are fixed for a given experiment), and Ω = {2πk/N, k = 0, . . . , N − 1} is a set of frequency lines. Assuming the model plant of Fig. 1, the systematic error is then found to have the expression
Fig. 7. Noise assumption for the H2 estimator.
Interestingly, the variance due to leakage is now an O(N −3 ) for the coherence function instead of an O(N −1 ), as previously found for the FRF.
ˆ k) ≈ Bias G(ω
N −1 l=0
|Al |2 |W (ωl − ωk )|2 [G(ωl ) − G(ωk )] N −1 l=0
|Al |2 |W (ωl − ωk )|2 (28)
B. H2 Function The H2 function is the mean-square estimator of a FRF in the presence of input noise only, as illustrated in Fig. 7. Denoting by x(n) = u(n) + νx (n) the noisy measurement of the input signal, it is given by M −1 |Ywi (ωk )|2 ˆ G2 (ωk ) = M −1i=0 . (24) ∗ i=0 Xwi (ωk )Ywi (ωk ) Under Assumptions A1–A4 and provided that v(τ ) has a ˆ 2 (ωk ) due to bounded derivative, the bias and variance of G leakage are (see the Appendix) 2 2 C (1) (ω )| |G k ˆ 2 (ωk ) −−−−→ v ) G (ωk )+2 Bias G N τc 2N 2 G(ωk )∗ (25) 2 2 ˆ 2 (ωk ) = Cv (θ) |G (ωk )|2 + Dv (θ) · NSRin (ωk ), Var G M (θ)N 2 M (θ) ωk = mod(π)
(26)
where NSRin (ωk ) = S2νx (ωk )/S2u (ωk ) stands for the input NSR at frequency ωk , and where the correction factors Cv2 (θ) ˆ k ) (output-noiseand Dv2 (θ) are identical to those found for G(ω only case). C. Conclusion It is noteworthy that (21), (22), (25), and (26) all involve the same structure as (16) and (17). In conclusion, those settings ˆ k )—i.e., a half-sine or that were found to be optimal for G(ω a diff window with about 2/3 overlap—are also optimal for measuring the coherence function and the H2 function. VI. M ULTISINE E XCITATIONS So far, our analysis has focused on the case where the excitation signal is random stationary. A related question of interest is whether the same reasoning can be followed to find the optimal settings in the case of a multisine excitation. For that purpose, let us write the multisine excitation signal of period N as Ak ejωk n+φk (27) u(n) = u(n + N ) = ωk ∈Ω
with W (ω) as the DTF of the data window w(n). Clearly, the bias is exactly zero provided that W (ωl − ωk ) = 0 for any ωl = ωk , that is, provided that W (ω) has zeros at all frequencies ω = 2πk/N , k = 1, . . . , N − 1. For an N -long window, this can only be achieved by the rectangular shape. Similarly, the expression of the variance due to leakage noise can be found to depend on factors involving W (ωl − ωk ), so that those windows that cancel the systematic error also cancel the stochastic leakage errors. Finally, the variance due to measurement noise is found to be 2 ˆ k ) ≈ Fv (θ) · Bw · NSRout (ωk ) Var G(ω M (θ) ∼ θFv2 (θ) · O(L−1 )
(29)
where NSRout (ωk ) = S2νy (ωk )/|Ak |2 stands for the output NSR, Bw = N R2v (0)/|W (0)|2 stands for the bandwidth of the data window (e.g., Bw = 1/N for the rectangular window), and Fv2 (θ) = 1 + 2
M (θ)−1
i=1
1−
i M (θ)
κ(iθ)
(30)
is a positive correction factor very similar to (14) but without the power 2 on κ(iθ). Contrary to the previous cases, it can be checked that Fv2 (θ) is now minimum at θ = 1, that is, when there is no overlap. Moreover, the window shape that minimizes (29) is that with the minimum bandwidth Bw , which is again the rectangular window. This is illustrated in Fig. 8. It is also shown in this figure that overlapping the rectangular window has some unexpected consequences due to the “bouncing” behavior of Fv2 (θ): the worst choice is 1/4 overlap; suboptimal choices are 1/2, 2/3, 3/4, . . . overlap, but with an increased computational cost. This is summarized in Table III. In conclusion, the WOSA optimal settings for measuring FRFs in the case of multisine excitations with period N is a rectangular window of length N without overlap. Such a choice completely cancels the impact of leakage noise and minˆ k )} ≈ imizes that of measurement noise (in that case, Var{G(ω NSRout /L). As expected, this conclusion is fully coherent with that of the maximum-likelihood approach presented in [1].
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
ANTONI AND SCHOUKENS: OPTIMAL SETTINGS FOR MEASURING FRFs WITH WOSA
3283
Fig. 9. Simulated system. TABLE IV MEASUREMENT ERRORS AT THE RESONANCE FREQUENCY OF A ONE-DEGREE-OF-FREEDOM OSCILLATOR Fig. 8. Bw Fv2 (θ)/M(θ) as a function of the amount of overlap (1 − θ) for the (continuous line) rectangular and (dotted line) half-sine windows, with L = 10 000 and N = 100. Note that M (θ) = (L − N )/N θ + 1. TABLE III CORRECTION FACTORS FOR THE RECTANGULAR WINDOW
with i as the index of the experiment. In each case, the window length was set to N = 210 , which is about four times as great as the system time constant τc . A. Validation of Bias and Variance Formulas VII. N UMERICAL V ERIFICATIONS The results established in the preceding sections are now verified on simulations. The simulated system is a one-degreeof-freedom oscillator with the transfer function G(z) =
1 − 2z −1 + z −2 1 − 1.6079z −1 + 0.9875z −2
(31)
which shows a sharp resonance at ω0 = 0.2π. According to (9), its corresponding time constant τc is about 223 samples. The system was excited by a random Gaussian input of unit variance and of length L = 214 , as illustrated in Fig. 9. Note that significant bias and variance are expected at the resonance frequency ω0 according to (16) and (17), which involve the first and second derivatives of G(ω) (Table IV). To make the computation of experimental bias and variance possible, the same experiment was repeated 1000 times but with different (independent) inputs, so that, concerning the measurement of the FRF, for instance, its experimental bias and variance could be computed as follows: ˆ k) = G(ω Bias
ˆ k) = G(ω Var
1000 1 ˆ G(ωk ; i) − G(ωk ) 1000 i=1
2 1000 1000 1 ˆ 1 ˆ G(ωk ; i) G(ωk ; i) − 1000 i=1 1000 i=1
Our primary concern is to check the validity of the main results (16), (17), (21), (22), (25), and (26). No measurement noise was added to better emphasize the impact of leakage noise only on the systematic and stochastic errors. The parameters for the WOSA measurements were set with a half-sine window and 2/3 overlap. Fig. 10 compares the experimental bias and standard deviation with the corresponding formulas provided in this paper, for the coherence function, the classical FRF estimate, and the “H2 ” function. A very good match is generally noticeable. As expected from our formulas, the maximum errors occur in the vicinity of the resonance frequency where the leakage effect is the strongest. Concerning the FRFs, the bias error becomes so considerable in that frequency region that it overruns the stochastic errors (i.e., standard deviation; see Fig. 10(b) and (c)]. As for the coherence function [see Fig. 10(a)], its bias error largely dominates over the whole frequency range; this is because its variance was found to be very small: an O(N −3 ) instead of an O(N −1 ) for the FRFs [see (17), (22), and (26)]. Finally, it is noteworthy that the bias of the “H2 ” function suddenly drops at exactly ω0 . Again, this behavior can be well predicted by plugging (31) into (25), and it suggests that “H2 ” may be an excellent estimator of the resonance amplitude (this is true, in general, for any resonant system). B. Validation of the Optimal Settings Our second concern is to experimentally check the impact of overlap and of the window shape on the WOSA measurement.
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
3284
IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 58, NO. 9, SEPTEMBER 2009
Fig. 10. Comparison of (continuous lines) experimental and (dotted lines) theoretical biases and standard deviations in the case of a half-sine window ˆ 0 ). (c) G ˆ 2 (ω0 ). Bold lines depict the with 2/3 overlap. (a) γ ˆ 2 (ωk ). (b) G(ω actual quantities to be measured.
A first experiment was conducted to compare the estimation errors in the cases without overlap and with 2/3 overlap. No measurement noise was added, and a half-sine window was used. The experimental biases and standard deviations are shown in Fig. 11. It is obvious from these results that 2/3 overlap can reduce the standard deviation of stochastic errors by almost a factor of 2. Precisely, the case without overlap yields M (θ = 1) = 16 segments and, according to Table I, a correction factor Cv2 (1) = π 2 , whereas the case with 2/3 overlap yields M (θ = 1/3) = 46 and Cv2 (1/3) = 9.61(1 + 0.13/46); thus, the reduction in standard deviation is like the square root of Cv2 (1)/M (1) over Cv2 (1/3)/M (1/3), that is, 1.7. In addition, it is checked that overlap has no effect on the bias error, as predicted from (16), (21), and (25). A second experiment was conducted to obtain more quantitative results with respect to both overlap and the window shape. To summarize each simulation by a single couple of points, the experimental bias and variance results were integrated over the whole positive frequency axis, i.e.,
Bias2 =
L/2−1
2 ˆ k ) G(ω Bias
k=1 L/2−1
Var =
k=1
ˆ k) . G(ω Var
Fig. 11. Comparison of experimental biases and standard deviations (dotted lines) without overlap and (continuous lines) with 2/3 overlap. (a) γ ˆ 2 (ωk ). ˆ ˆ (b) G(ω0 ). (c) G2 (ω0 ). Bold lines depict the actual quantities to be measured.
Note that, according to these definitions, the total MSE is simply obtained as Bias2 + Var. The corresponding results are displayed in Fig. 12. As expected, it is shown that the squared bias is (statistically) constant against overlap, whereas the variance is a monotonically decreasing curve, with its minimum nearly reached at 2/3 overlap. This is independently of the window shape. These experimental curves fairly well reproduce the trends of Fig. 5, notwithstanding a different scaling factor relating to the system characteristics. Indeed, from (17), it theoretically holds that Var
L/2−1 Cv2 (θ) 2 |G (ωk )| M (θ)N 2 k=1
where Cv2 (θ)/M (θ) is the quantity displayed in Fig. 5. This again validates (17) and its implication about the optimality of the 2/3 rule. Finally, Fig. 12 demonstrates that the half-sine window is much better than the Hanning window as long as leakage errors (bias and variance) are of concern. A reduction of about 1.8 is observed on the squared bias, which is in perfect accordance with what can be predicted from (16) and Tables I and II, i.e., 2 2 C (1)Hanning 2 = 4 1.778. = v2 Cv (1)half -sine 3 Bias2 (half-sine) Bias2 (Hanning)
The reduction on the variance term is less significant, particularly at high levels of overlap. Equation (17) and Tables I and II,
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
ANTONI AND SCHOUKENS: OPTIMAL SETTINGS FOR MEASURING FRFs WITH WOSA
Fig. 12. Integrated squared bias and variance due to leakage with respect to the amount of overlap for the Hanning and half-sine windows.
in the case of no overlap, give Var(Hanning) C 2 (1)Hanning 4 = v2 = 1.33 Cv (1)half -sine 3 Var(half-sine) which fairly well compares with the reduction of about 1.4, as shown in Fig. 12, and in the case of 2/3 overlap, i.e., 0.27 10.50 1 + 2 M (1/3) Var(Hanning) Cv (1/3)Hanning = = Var(half-sine) Cv2 (1/3)half -sine 9.61 1 + 0.13 M (1/3)
1.10, which again is coherent with Fig. 12. Very similar results were obtained for the measurement of the coherence and the H2 functions, which are not displayed here due to the lack of space. C. Investigation of the Effect of Measurement Noise Our last concern is to check the scope of validity of our conclusions when the SNR becomes low. To do so, the total MSE Bias2 + Var was computed for different output SNRs (see Fig. 1), and the values were reported for both the Hanning and half-sine windows at various amounts of overlap. The experimental results are displayed in Fig. 13. It is shown that the half-sine window still outperforms the Hanning window, until the SNR reaches the unrealistically low value of −40 dB. At this stage, measurement noise completely overwhelms leakage noise, and the curves of Fig. 13(d) actually reproduce the trends of Fig. 6. This demonstrates that, at least for the very resonant system simulated in this section, leakage errors are very predominant (actually mainly the bias) that the half-sine window remains a very good choice, even in the presence of (realistic levels of) additive noise. VIII. C ONCLUSION The measurement of FRFs by means of the WOSA procedure suffers from three types of errors, namely, a systematic error
3285
Fig. 13. Overall MSE as a function of the SNR for the (· · · • · · ·) Hanning and (· · · ◦ · · ·) half-sine windows. (a) SNR = ∞ dB. (b) SNR = 0 dB. (c) SNR = −20 dB. (d) SNR = −40 dB.
due to leakage, a stochastic error due to leakage noise, and a stochastic error due to measurement noise. We have proven that the optimal real symmetric data-window shape that globally minimizes the leakage errors is the half-sine window, which incidentally achieves all the same properties as the (complexvalued and nonsymmetric) diff window recently proposed in [17]. Moreover, the stochastic errors can significantly be reduced by increasing the percentage of overlap between adjacent segments. The ultimate variance reduction that can achieved this way is about 3 on the leakage errors and 1.7 on the measurement noise errors with the half-sine window. An important practical result is that, in all cases, this is nearly achieved by setting about 2/3 overlap independently of the system FRF and of the data window. The “2/3 law” should replace the customary 1/2 law used in commercial data analyzers without increasing too much the computational demand. Moreover, it has been shown that the same recommendations apply as well to the measurement of other spectral quantities and, in particular, to the coherence function, which is often used in conjunction with FRF measurements. As a final note, we emphasize again that the optimality of these settings applies provided that the data-window length is set larger than the system time constant—which is in any case the usual recommendation with WOSA—and only with stationary random excitations. The special case of multisine excitations was verified to yield completely different conclusions. A PPENDIX The sketch of our proofs starts from the demonstration of relation (6). From Fig. 3, the leakage error due g(k)u(n − to segment i is twi (n) = y(n)w(n − iR) − ∞ k=0 k)w(n − k − iR), i.e.,
Twi (ωk ) =
Wi (ωk − ν) (G(ν) − G(ωk )) dU (ν)
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
(32)
3286
IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 58, NO. 9, SEPTEMBER 2009
where Wi (ωk − ν) = W (ωk − ν)e−jiR(ωk −ν) , and dU (ν) is the spectral increment of process u(n) at frequency ν [6]. Therefore, Ywi (ωk ) = G(ωk )Uwi (ωk ) + Twi (ωk ) + Nwi (ωk ), from which (6) immediately follows with M −1
i=0 βL (ωk ) = M −1 i=0
M −1
βN (ωk ) = i=0 M −1 i=0
Twi (ωk )Uwi (ωk )∗
(33)
Uwi (ωk )Uwi (ωk )∗
D. Equations (21) and (22)
∗
Nwi (ωk )Uwi (ωk )
Uwi (ωk )Uwi (ωk )∗
expected value of the denominator. Here, again, only the former quantity needs to be evaluated. From P2, P3 with p = 0, P4, P1, −1 ∗ and Assumptions A1–A3, Var{ M i=0 Nwi (ωk )Uwi (ωk ) } = 2 2 2 2 S2u (ωk )S2ν (ωk )(2π) M Dv (θ)R2v (0)N . Taking the ratio and using definition (12), (18) follows.
.
(34)
Now, from Assumption A4, Twi (ωk ) may be approximated as
1 2 Wi (ωk − ν) (ν − ωk )G (ν) − (ωk − ν) G (ν) dU (ν). 2 (35)
The evaluation of (21) and (22) starts with replacing Ywi (ωk ) in (20) by its expansion G(ωk )Uwi (ωk ) + Twi (ωk ). From Assumption A4, the magnitude of the transient term Twi (ωk ) is much smaller than Uwi (ωk ) on the average. Expanding ratio (20) to the second order with respect to Twi (ωk ) then yields the approximation M −1
2
γˆ (ωk ) 1 − i=0 M −1 i=0
All proofs proceed by making use of (33) and (34) and a set of properties: (P1): (P2): (P3): (P4):
E{dU (ν1 )dU (ν2 )} = S2u (ν1 )δ(ν1 + ν2 )dν1 dν2 . E{dU (ν1 )dN (ν2 )} = 0. (p) |W (ν)|2 ν p ≈ (−j)p 2πN (1−p) R2v (0), p = 0, 1, 2. The input signal u(n) is assumed to be Gaussian.
Property P1 follows from Assumption A1, Property P2 from Assumption A3, and Property P3 from the assumed differentiability of the window shape v(τ ). Property P4 is invoked here for simplicity, although not fundamentally necessary since the DFT makes (most) signal tend to Gaussianity in virtue of the central limit theorem [6].
From Assumption A4, the expected value of the ratio (33) tends to the ratio of the expected values. Namely, from P1, with p = 0, and Assumptions A1 and P3 −1 ∗ } = S2u (ωk )2πN R2v (0)M , A2, E{ M i=0 Uwi (ωk )Uwi (ωk ) −1 ∗ and from P3 with p = 2, E{ M i=0 Twi (ωk )Uwi (ωk ) } = −1 −S2u (ωk )πN R2v (0)M G (ω). Taking the ratio and using definition (12), (16) immediately follows. B. Equation (17) From Assumption A4, the variance of the ratio (33) tends to the ratio of the variance of the numerator to the square of the expected value of the denominator. Since the latter has just been evaluated in the preceding discussion, we only need to work out the former. Namely, from P1, P3−1with p = 0, 1, P4,∗ and Assumptions A1 and A2, Var{ M i=0 Twi (ωk )Uwi (ωk ) } = 2 2 (ωk )(2π)2 M Cv2 (θ)R2v (0)|G (ω)|2 . Taking the ratio and S2u using definition (12), (17) immediately follows.
C. Equation (18) From Assumption A4, the variance of the ratio (34) tends to the ratio of the variance of the numerator to the square of the
|Ywi (ωk )|2
.
(36)
The rest of the proof follows as for the bias and variance of ˆ k ) by taking the expected value and the variance of (36), G(ω simplifying the statistics of a ratio as the ratio of the statistics due to Assumption A4, and using Assumptions A1–A2 and Properties P1–P4.
E. Equations (25) and (26) The evaluation of (25) and (26) proceeds as for the H2 function; expansion of the ratio (24) to the second order with respect to Twi (ωk ) then yields the approximation
A. Equation (16)
|Twi (ωk )|2
G2 (ωk ) G(ωk ) ·
M −1
1+
∗ i=0 Twi (ωk )Ywi (ωk ) M −1 2 i=0 |Ywi (ωk )|
. (37)
The rest of the proof follows as before, by taking the expected value and the variance of (37), simplifying according to Assumption A1–A4, and using Properties P1–P4. R EFERENCES [1] R. Pintelon and J. Schoukens, System Identification—A Frequency Domain Approach. Piscataway, NJ: IEEE Press, 2001. [2] P. Welch, “The use of the fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms,” IEEE Trans. Audio Electroacoust., vol. AU-15, no. 2, pp. 70–73, Jun. 1967. [3] R. B. Blackman and J. W. Tukey, The Measurement of Power Spectra From the Point of View of Communications Engineering. New York: Dover, 1959. [4] G. M. Jenkins and D. G. Watts, Spectral Analysis and Its Applications. Boca Raton, FL: Emerson-Adams Press. [5] J. Bendat and A. Piersol, Random Data: Analysis and Measurement Procedures, 2nd ed. New York: Wiley-Interscience, 1986. [6] D. R. Brillinger, Time Series: Data Analysis and Theory. New York: SIAM, 1975. [7] A. H. Nuttall and G. C. Carter, “A generalized framework for power spectral estimation,” IEEE Trans. Acoust., Speech, Signal Process., vol. ASSP-28, no. 3, pp. 334–335, Jun. 1980. [8] A. H. Nuttall, “Spectral estimation of means of overlapped FFT processing of windowed data,” Naval Underwater Syst. Center, New London, CT, Rep. 4169. suppl. TR4169S.
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
ANTONI AND SCHOUKENS: OPTIMAL SETTINGS FOR MEASURING FRFs WITH WOSA
[9] A. H. Nuttall, “On the weighted overlapped segment averaging method for power spectral estimation,” Proc. IEEE, vol. 68, no. 10, pp. 1352–1354, Oct. 1980. [10] D. R. Brillinger, “The key role of tapering in spectrum estimation,” IEEE Trans. Acoust., Speech, Signal Process., vol. ASSP-29, no. 5, pp. 1075– 1076, Oct. 1981. [11] F. J. Harris, “On the use of windows for harmonic analysis with the discrete Fourier transform,” Proc. IEEE, vol. 66, no. 1, pp. 51–83, Jan. 1978. [12] A. Eberhard, “An optimal discrete window for the calculation of power spectra,” IEEE Trans. Audio Electroacoust., vol. AU-21, no. 1, pp. 37–43, Feb. 1973. [13] N. Geçkinli and D. Yavuz, “Some novel windows and a concise tutorial comparison of window families,” IEEE Trans. Acoust., Speech, Signal Process., vol. ASSP-26, no. 6, pp. 501–507, Dec. 1978. [14] A. Nuttall, “Some windows with very good sidelobe behavior,” IEEE Trans. Acosut., Speech, Signal Process., vol. ASSP-29, no. 6, pp. 84–91, Feb. 1981. [15] J. Le Roux and J. Ménez, “A cost minimization approach for optimal window design in spectral analysis of sampled signals,” IEEE Trans. Signal Process., vol. 40, no. 4, pp. 996–999, Apr. 1992. [16] J. L. Douce and L. Balmer, “Statistics of frequency-response estimates,” Proc. Inst. Elect. Eng., vol. 137, no. 5, pp. 290–296, Sep. 1990. [17] J. Schouken, Y. Rolain, and R. Pintelon, “Analysis of windowing/leakage effects in frequency response function measurements,” in Proc. 16th IFAC World Congr., Prague, Czech Republic, Jul. 4–8, 2005.
3287
Jérôme Antoni was born in Strasbourg, France, on November 11, 1972. He received the Ph.D. degree (with highest honors) in signal processing from the Polytechnic Institute of Grenoble, Grenoble, France, in 2000. He is currently a Lecturer with the University of Technology of Compiègne, Compiègne, France. His research interests include nonstationary signal analysis, system identification, and signal separation. His applications are concerned with noise and vibration analysis, modal analysis, and system diagnosis.
Johan Schoukens (M’90–SM’92–F’97) received the degree in engineering and the Ph.D. degree in applied sciences from the Vrije Universiteit Brussel (VUB), Brussels, Belgium, in 1980 and 1985, respectively. He is currently a Professor with the VUB. The prime factors of his research are in the field of system identification for linear and nonlinear systems. Dr. Schoukens was the recipient of the Best Paper Award and the Society Distinguished Service Award from the IEEE Instrumentation and Measurement Society in 2002 and 2003, respectively.
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 58, NO. 9, SEPTEMBER 2009
Optimal Settings for Measuring Frequency Response Functions With Weighted Overlapped Segment Averaging Jérôme Antoni and Johan Schoukens, Fellow, IEEE
Abstract—This paper investigates the measurement errors involved in estimating frequency response functions (FRFs)—and related quantities such as the coherence function—from weighted overlapped segment averaging, a technique that has become a standard in modern data analyzers due to its computational advantages. Particular attention is paid to leakage errors, for which this technique has frequently been criticized. Our main result is that a half-sine or diff window with about 2/3 overlap achieves the best compromise to reduce leakage errors in the case of stationary random excitations, and this is independently of the system FRF. This conclusion is to be contrasted with the customary habit of using a Hanning window with 1/2 overlap. The same reasoning confirms that a rectangular window without overlap minimizes measurement errors in the special case of multisine excitations. Moreover, practical formulas are provided to the reader for computing the bias and variance of the frequency response and coherence functions in the general case. Index Terms—Frequency response function (FRF), optimal window, spectral leakage, transient effects, weighted overlapped segment averaging (WOSA).
I. I NTRODUCTION
T
HE MEASUREMENT of frequency response functions (FRFs) from experimental data is of continuous interest in numerous engineering fields and, more specifically, in any application where one is concerned with inferring the structure of an unknown linear system from its measured inputs and outputs. Although the parametric identification of the system is usually the ultimate objective—on which there exists a vast literature—the nonparametric estimation of FRFs offers a simple and yet very useful preliminary description of the system in the frequency domain. For instance, measured FRFs—or related quantities such as the coherence function—can be used to better understand a system order or to check its (non)linearity. Modern FRF measurement techniques are extensively based on the discrete Fourier transform (DFT), owing to its very efficient computation by means of the fast Fourier transform (FFT) algorithm. This flexibility, plus the constantly increasing capacity of computational devices, has led to the proposal of several estimators whose relative optimality strongly depends on the statistical nature of the excitation signals (i.e., inputs; see Manuscript received May 6, 2008; revised April 9, 2009. Current version published August 12, 2009. J. Antoni is with the Laboratory Roberval of Mechanics, University of Technology of Compiègne, 60205 Compiègne, France (e-mail: [email protected]). J. Schoukens is with the Department of Electrical Engineering, Vrije Universiteit Brussel, 1050 Brussels, Belgium. Digital Object Identifier 10.1109/TIM.2009.2022376
[1] for a comprehensive overview on the subject). Surely, the best possible choice is to use the maximum-likelihood estimator together with periodic excitation signals, since only in that case can the DFT be exactly used without finite-length effects. There may be, however, some technical or psychological reasons that lend the experimenter to prefer other types of excitations. Among those, a popular subclass is encompassed by stationary random signals (e.g., binary signal, random uniform, and random Gaussian). Different FRF estimators have been proposed over the years in this context, which essentially find their roots in experimental spectral analysis. They all rely on various frequency smoothing or averaging strategies, with the fundamental assumption (hope) that the estimator variability can so be stabilized without distorting too much the smooth structure of the FRF. In any case, this obviously leads to a difficult bias/variance tradeoff, which nonparametric FRF measurement techniques have always been criticized for. In this paper, we focus on the so-called weighted overlapped segment averaging (WOSA) or Welch’s procedure. The WOSA procedure consists of segmenting the available records into several short-length segments, computing their DFTs, and then averaging their products to get a stable estimator of the system FRF. Although the idea of the WOSA procedure is quite old [2], it has surprisingly received very little attention in the scientific literature. When it comes to the FRF measurement, its exact statistical performance and optimal settings remain virtually unknown. It is the objective of this communication to partially fill in these gaps. In particular, we address important and practical questions. 1) Which window shape should we use to taper the data? 2) How much can we gain from overlapping adjacent segments? 3) Which optimal amount of overlap should we set between adjacent segments? The answers to these questions will obviously provide the experimenter with some useful guidelines to make optimal usage of the WOSA technique in FRF measurements. This paper is organized as follows. In Section II the principle and limitations of the WOSA procedure are reminded. The statistical errors involved by the WOSA procedure are quantified in Section III, in terms of systematic and stochastic errors due to leakage and measurement noise. From these results, the optimal settings of the WOSA parameters are then deduced in Section IV. Section V shows that the conclusions remain unchanged for the WOSA estimation of other spectral
0018-9456/$26.00 © 2009 IEEE
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
ANTONI AND SCHOUKENS: OPTIMAL SETTINGS FOR MEASURING FRFs WITH WOSA
Fig. 1.
3277
Plant model.
quantities such as the coherence function. Section VI proves that the same reasoning can be applied to the special case of multisine excitations, although the optimal settings are then found to be very different. Finally, Section VII presents some numerical experiments. The proofs of our main results are sketched in the Appendix. II. FRF M EASUREMENT W ITH THE WOSA P ROCEDURE A. Classical Approach Let us consider a pair of input/output stationary random signals {u(n); y(n)}, n ∈ Z, and an assumed linear time-invariant relationship, i.e., y(n) =
∞
g(k)u(n − k) + νy (n)
Ywi (ωk ) =
where g(k) stands for the (causal) system impulse response, and νy (n) for the additive measurement noise at the output (see Fig. 1). The measurement of the system con FRF then −jωk g(k)e = sists of estimating the values of G(ω) = ∞ k=0 F{g(k)} at some specific frequencies ωk in a finite set Ω from L−1 the finite-length records {u(n)}L−1 n=0 and {y(n)}n=0 . The classical approach to this problem is to estimate G(ω) as the ratio of the smoothed cross-periodogram to the smoothed periodogram, i.e., ∗ M 2πi 2πi ω U ω g Y − − i L k k L i=−M L L ˆ k) = (2) G(ω M 2πi 2 | i=−M gi |UL ωk − L where L−1
y(n)e−jωk n = DFT{y(n)},
n=0
idea of the WOSA estimator (5) is to arrive at the same result by averaging over segments of data. Namely, let us denote
(1)
k=0
YL (ωk ) =
Fig. 2. Principle of the segmenting and windowing process. The signal of length L is divided into M (weighted) segments of length N that are overlapping by N − R samples.
ωk =
2πk L (3)
is the DFT of {y(n)}L−1 n=0 —and similarly for UL (ωk )—and where the gi ’s, i = −M, . . . , M , are the taps of a suitably designed low-pass filter. The FRF estimator (2) is directly inherited from the early works of Blackman and Tukey in power spectrum estimation [3] (it is also known as the Blackman ˆ k ) is the and Tukey procedure). Note that the so-defined G(ω weighted mean-square estimator of G(ωk ) given the set of (noisy) observations {YL (ωk − 2πi/L)}M i=−M and the set of re. The statistical performance gressors {UL (ωk − 2πi/L)}M i=−M of this approach has been thoroughly studied in references such as [4]–[6]. B. WOSA Procedure Whereas the Blackman and Tukey estimator (2) controls stability by smoothing over adjacent frequency lines, the basic
iR+N −1
w(n − iR)y(n)e−jωk n ,
n=iR
ωk =
2πk N (4)
−1 the short-time DFT of the segment of data {y(n)}iR+N ann=iR chored at time iR and tapered by the smoothed N -long window −1 {w(n)}N n=0 ; the WOSA FRF estimator is then obtained as the ratio of the averaged cross-periodogram to the averaged periodogram, i.e., M −1 Ywi (ωk )Uw∗ i (ωk ) ˆ G(ωk ) = i=0 (5) M −1 2 i=0 |Uwi (ωk )|
with M = (L − N )/R + 1 as the total number of data segments (x standing for the greatest integer less than or equal to x) and 1 − R/N as the fraction of overlap between adjacent segments (see Fig. 2). The WOSA procedure seems to have become a standard in commercial data analyzers in virtue of its simplicity and many computational advantages. 1) It can be very efficiently computed by means of several short FFTs (which is always faster than computing one very long FFT). 2) It requires few memory allocation, and it can actually be recursively computed in almost real time (this is the solution implemented in current commercial data analyzers). 3) It is robust against nonstationary, outlying, or missing data since the process of segmenting allows for easy detrending and removal of corrupted segments. The performance of WOSA has been investigated in some early papers within the context of spectrum estimation. In a series of papers [7]–[9], Nutall showed that the WOSA procedure can be made nearly as efficient (in the sense of achieving the same statistical stability given a similar frequency resolution) as the Blackman and Tukey procedure provided that sufficient
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
3278
IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 58, NO. 9, SEPTEMBER 2009
Fig. 3. Transient effects with the (left) rectangular window and (right) Hanning window.
overlap between adjacent segments is applied. From numerical simulations, Nutall reported that ideally 62.5% overlap should be imposed, but that 50% is more practical and almost just as efficient. This recommendation seems to have become the standard nowadays. In another series of papers, the key role of tapering the segments with a smooth data window (as opposed to the rectangular window) was clearly demonstrated (see [10] and the references therein). Nevertheless, there has been no unique answer so far to the question of an optimal window shape; instead, many “good” windows have been proposed over the years, each with its own peculiarities [11]–[15]. Surprisingly, there are very few related results when it comes to applying WOSA for the measurement of FRFs [16], [17]. Indeed, this issue happens to be more involved than that of spectral analysis. One major peculiarity of the WOSA estimator when used for FRF measurements is that the use of several segments of reduced length makes it prone to severe leakage errors. If the mechanism of leakage due to Fourier transforming a truncated time series is well known in spectral analysis, its implication in the WOSA technique is more subtle and requires some explanations. It stems from the inherent but erroneous assumption that the tapered output segments are identical to the system responses of the tapered input segments, despite of the nonlinearity of the windowing process. This produces transient artifacts in the signals at the onset and end of each segment, which will perturb the FRF estimation process. As illustrated in Fig. 3, a typical solution to reduce the transient magnitude is to use a smooth data window—e.g., the Hanning window as opposed to the rectangular window—but this is generally at the price of an increase in the transient duration so that some compromise has to be found. Another intuitive solution is to increase overlap between adjacent segments so that transient artifacts got averaged out to a further extent. We now address such questions as whether there is an optimal window shape to taper the data with, as well as how much can be gained from overlapping adjacent segments. The answers to these questions are nontrivial, and when it comes to FRF measurements, the approach to obtain them is quite different from that used in classical spectral analysis. The proposed approach proceeds by first analyzing the measurement errors involved in the WOSA procedure. These are shown to stem from three sources, such as the following: 1) systematic errors; 2) stochastic errors due to spectral leakage; and 3) stochastic errors due to additive measurement noise.
This can be summarized by the following remarkably simple formula (see the proof in the Appendix) ˆ k ) = G(ωk ) + βL (ωk ) + βN (ωk ) G(ω
(6)
where βL (ωk ) is a nonzero-mean stochastic error accounting for spectral leakage, and βN (ωk ) is a zero-mean stochastic error accounting for measurement noise. We then pay particular attention to the leakage error βL (ωk ), the statistical analysis of which provides closed-form (asymptotic) expressions for the bias and the variance of the WOSA FRF as a function of the data-window shape and length, the number of segments, and the amount of overlap between adjacent segments. Incidentally, (6) proves that leakage affects both the bias and ˆ k ) since the variability of G(ω ˆ k ) = G(ωk ) + E {βL (ωk )} (7) E G(ω ˆ k ) = Var {βL (ωk )} + Var {βN (ωk )} . (8) Var G(ω It seems that only the deterministic nature of leakage (i.e., the first term on the right-hand side of (7) has been recognized in classical references on the subject) and that its random contribution has been ignored from variance analyses [2], [5]. However, there are many instances where the latter dominates the former [17]. III. S TATISTICAL E RRORS A. Assumptions and Notations In this section we assume, for simplicity, that A1) the input signal u(n) is real, zero mean, and stationary; A2) its correlation length is much smaller than the system time constant1 τc =
∞ 1 2 2 k=0 |g(k)|k ∞ k=0 |g(k)|
(9)
as expected from a well-designed input signal [1]; A3) the measurement noise is independent of u(n); 1 Note that this definition of the time-constant encompasses both notions of pure delay and of group delay of the impulse response of the system.
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
ANTONI AND SCHOUKENS: OPTIMAL SETTINGS FOR MEASURING FRFs WITH WOSA
3279
TABLE I CORRECTION FACTORS FOR THE HANNING WINDOW
TABLE II CORRECTION FACTORS FOR THE HALF-SINE WINDOW Fig. 4. (Left) Hanning and (right) half-sine windows v(τ ) and their corresponding (normalized) autocorrelation functions κ(τ ).
A4) the system time constant τc is much smaller than the data-window length (as classically recommended with WOSA and more generally with any nonparametric measurement of FRFs). Furthermore, to propose concise expressions, let us introduce the continuous window shape v(τ ) defined over the support set [0; 1[. The discrete N -long data window w(n) used in WOSA [e.g., in (5)] is obtained by sampling v(τ ) with N samples per unit length, i.e., n = 0, . . . , N − 1.
w(n) = v(n/N ),
of segments, and the window shape v(τ ). Note that Cv2 (θ) exists only for windows having a bounded derivative, a requirement that precludes, in particular, the rectangular window.
(10) B. Systematic Error
Let us also denote 1−|τ
|
v (|τ | + x) v(x)dx
R2v (τ ) =
(11)
0
Under Assumptions A1–A4 and provided that v(τ ) has a ˆ k )} of (7) bounded derivative, the systematic error Bias{G(ω can be shown to have the asymptotic expression2 (see the Appendix)
the autocorrelation function of v(τ ), and
E {βL (ωk )} −−−−→ N τc
R2v (τ ) κ(τ ) = R2v (0)
(12)
its normalized version (see Fig. 4). Finally, let 100(1 − θ) be the percentage of overlap between two adjacent segments, where θ = R/N (see Fig. 2). With these notations, we are now in position to define the two correction factors Cv2 (θ)
M (θ)−1
= − κ (0) − 2
i=1
i 1− M(θ)
κ (iθ)κ(iθ) (13)
Dv2 (θ) = 1 + 2
M (θ)−1
i=1
1−
i M (θ)
κ(iθ)2
(14)
which are explicit functions of the percentage of overlap, the number M (θ) = (L − N )/N θ + 1
(15)
Cv2 (1) G (ωk ) ∼ O(N −2 ) 2N 2
(16)
where the correction factor Cv2 (1) = −κ (0) is as given in (13) with θ = 1 and may be interpreted as the curvature of the autocorrelation function of the window at lag 0—see the first row in Tables I and II for some particular values—and G (ωk ) is the second derivative of G(ωk ). Expression (16) is instructive in two instances: First, since the bias is an O(N −2 ), long data windows should be used. Second, since it is proportional to −κ (0), data windows implying minimum curvatures should be used. This issue will be addressed in the next section. C. Stochastic Errors From (8), the total variance of the WOSA FRF estimator decomposes into one contribution due to leakage and one contribution due to measurement noise. We now provide the general expressions for these two contributions given an arbitrary 2 The
notation larger than τc .”
−−−−−−−−−→(·) N τc
means “converges to (·) as N becomes much
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
3280
IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 58, NO. 9, SEPTEMBER 2009
percentage of overlap and an arbitrary window shape. The sketch of the proofs is provided in the Appendix. 1) Stochastic Errors Due to Leakage: Under Assumptions A1–A4 and provided that v(τ ) has a bounded ˆ k ) due to leakage is (see the derivative, the variance of G(ω Appendix) Var {βL (ωk )} −−−−→ N τc
Cv2 (θ) 2 |G (ωk )| , ωk = 0 mod (π) M(θ)N 2 ∼ θCv2 (θ) · O(L−1 N −1 )
(17)
with the correction factor Cv2 (θ) as given in (13). Equation (17) can be somehow simplified when particular cases are considered: we deemed it useful to provide the reader with the equivalent formula reported in Tables I and II for the typical cases of a Hanning and a Half-sine window, respectively, with 0, 1/2, and 2/3 overlap. Since the variance is shown to be an O(L−1 N −1 ), one solution to reduce stochastic errors due to leakage (given a fixed signal length L) is to use long data windows. Another solution is to set the fraction of overlap θ = R/N that minimizes Cv2 (θ). This issue will be addressed in the next section. 2) Stochastic Errors Due to Measurement Noise: Under ˆ k ) due to measureAssumptions A1–A3, the variance of G(ω ment noise is (see the Appendix) Var {βN (ωk )} ≈
Dv2 (θ) · NSRout (ωk ) M (θ) ∼ θDv2 (θ) · O(L−1 N )
(18)
with the correction factor Dv2 (θ) as given in (14), and where NSRout (ωk ) = S2νy (ωk )/S2u (ωk ) stands for the output noiseto-signal ratio (NSR) at frequency ωk . Again, the reader may find it convenient to use the equivalent formulas reported in Tables I and II for the typical cases of a Hanning and a halfsine window, respectively, with 0, 1/2, and 2/3 overlap. Expression (18) is similar to that reported in [2]. However, we emphasize again the fact that it cannot alone explain the variability of the WOSA FRF estimator. IV. O PTIMAL S ETTINGS OF THE WOSA P ARAMETERS We now make use of the preceding results (16)–(18) to find the optimal parameters to be used with the WOSA procedure. Our first result concerns the optimal window shape that is found to minimize the errors due to leakage (both systematic and stochastic). Our second result specifies how much overlap should be used to minimize the overall stochastic errors (both leakage and measurement noise). A. Optimal Window Shape The optimal window shape should minimize the following mean-square error (MSE): MSE(ωk ) = |E {βL (ωk )}|2 +Var {βL (ωk )}+Var {βN (ωk )} . For the sake of simplicity and to avoid having the optimum depending on the unknown signal-to-noise ratio (SNR) and
Fig. 5. Cv2 (θ)/M(θ) as a function of the amount of overlap (1 − θ) for the (dotted line) half-sine and (continuous line) Hanning windows, with L = 10 000 and N = 100. Note that M (θ) = (L − N )/N θ + 1.
actual FRF, we consider only the minimization of the errors due to leakage (i.e., the first two terms on the right-hand side of the preceding equation). This amounts to finding the optimal window in the noise-free case. Turning back to (16) and (17), it is shown that, in the particular case of no overlap (θ = 1), such a window should minimize the curvature −κ (0) of its autocorrelation function. In the family of real symmetric windows, it can be shown (e.g., by standard variational analysis) that such a minimum is achieved by the half-sine window shape (see Fig. 4), i.e., v(τ )−−−−→ sin(πτ ), N τc
0 ≤ τ ≤ 1.
(19)
The half-sine window has several appealing properties. 1) It has the same minimum-bias performance as the recently proposed (complex and asymmetric) “diff” window [17]. 2) Its optimality is asymptotically (when N τc ) independent of the system FRF. 3) It is incidentally very similar to the parabolic window historically used in Welch’s original paper [2]. Moreover, it can be shown that the optimality of the half-sine window still approximatively holds when overlap is increased (θ < 1). This is because the value of Var{βL (ωk )} becomes almost independent of the window shape as θ → 0, as shown in Fig. 5. In fact, the superiority of the half-sine window over the Hanning window can be checked in Fig. 4, wherein the larger curvature of κ(τ ) at τ = 0 implies a smaller bias according to (16), and in Fig. 5, wherein the smaller correction factor Cv2 (θ) implies a smaller variance according to (17). Precisely, Tables I and II) indicate a reduction of 1.33 on the bias and, more modestly, of 1.33, 1.39, and 1.09 on the variance with 0, 1/2, and 2/3 overlap respectively. This superiority may be
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
ANTONI AND SCHOUKENS: OPTIMAL SETTINGS FOR MEASURING FRFs WITH WOSA
3281
Moreover, Figs. 5 and 6 show that not much more can be gained by forcing adjacent segments to overlap more than about 2/3 and 1/2, respectively, because the added information then becomes extremely redundant. This means that, depending on the window shape and whether leakage noise or measurement noise predominates, the optimal fraction of overlap is somewhere between 1/2 and 2/3. Since the actual SNR will rarely be known, a conservative choice is about 2/3 overlap, a limit above which the extra computational cost involved by forcing more overlap is surely not justified. More accurate settings may be formulated on a case-by-case basis, depending on both the used window and the leakage- or noise-free case, but we believe that the simple and conservative “2/3 rule” is just good enough in practice—it can actually be shown to conservatively hold for all smooth (i.e., with a bounded derivative) data windows, just as observed here above on the Hanning and half-sine windows.3 V. E XTENSION TO O THER S PECTRAL F UNCTIONS Fig. 6. Dv2 (θ)/M(θ) as a function of the amount of overlap (1 − θ) for the (dotted line) half-sine, (continuous line) Hanning, and (gray line) rectangular windows, with L = 10 000 and N = 100. Note that M (θ) = (L − N )/N θ + 1.
explained by the fact that, in the frequency domain, the half-sine window achieves a better balance between local leakage (width of the main lobe) and global leakage (roll-off of the secondary lobes). Finally, it should be emphasized that the half-sine window is no longer optimal—in the sense of minimizing the overall MSE—in the presence of (significant) measurement noise. It can be seen from (14) and (18) that the reduction of Var{βN (ωk )} requires a window shape with small inertia (a Dirac impulse in the limit), which is conflicting with the reduction of the leakage errors. Therefore, the optimum window will be a compromise in the general case, depending on both the unknown SNR and FRF. Addressing this difficult issue is outside the scope of this paper and is also unlikely to provide the end user with any simple recommendation. B. Optimal Amount of Overlap Having set the window shape that minimizes the effect of leakage errors, it is still possible to significantly reduce the MSE by acting on the overlap. As mentioned in the preceding section, the intuition is that overlapping will average out stochastic errors to a further extent by providing a larger number M (θ) of segments. Indeed, turning back to (13) and (14) and their illustration in Figs. 5 and 6, it is shown that Cv2 (θ)/M(θ) and Dv2 (θ)/M(θ) are globally decreasing functions of the amount of overlap (1 − θ). More precisely, Tables I and II show that from 0 to 2/3 overlap, the number of windows is roughly tripled, whereas the correction factors Cv2 (θ) and Dv2 (θ) are only slightly changed. For instance, with a Hanning window, this yields a reduction by 3.8 for the variance due to leakage and by 2 for the variance due to measurement noise. With a half-sine window, the reduction is by 3 for the variance due to leakage and by 1.7 for the variance due to measurement noise.
It is a natural question to investigate whether the optimal settings found hitherto for measuring FRFs also apply to the measurement of other spectral functions. Of particular interest is the WOSA estimator of the (squared-magnitude) coherence function and the alternative WOSA measurement of FRFs—the so-called “H2 ” function. A. Coherence Function The coherence function 2 M −1 i=0 Ywi (ωk )Uw∗ i (ωk ) 2 γˆ (ωk ) = M −1 2 M −1 2 i=0 |Ywi (ωk )| i=0 |Uwi (ωk )|
(20)
is a useful quantity for checking the linearity of the system under study. For the sake of simplicity, we only address here the impact of leakage noise on its measurement (i.e., the noisefree case), which, to our knowledge, has never been investigated before in the literature. Investigation of the combined effect of leakage and measurement noise is left for future research. Under Assumptions A1–A4 and provided that v(τ ) has a bounded derivative, the bias and variance of γˆ 2 (ωk ) due to leakage are (see the Appendix)
C 2 (1) |G (ωk )|2 Bias γˆ 2 (ωk ) −−−−→ v 2 N τc 2N |G(ωk )|2
E 2 (θ) |G (ωk )|4 Var γˆ 2 (ωk ) −−−−→ v , N τc M (θ)N 4 |G(ωk )|4
(21) ωk = 0 mod (π)
∼ θEv2 (θ) · O(L−1 N −3 )
(22)
where Ev2 (θ) = κ (0)2 + 2
M (θ)−1
i=1
1−
i M (θ)
κ (iθ)2 .
(23)
3 The formal proof of this assertion is only of technical interest and will be published elsewhere; it can actually easily be verified by computing (13) and (14) for different window candidates.
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
3282
IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 58, NO. 9, SEPTEMBER 2009
where Ak are real deterministic amplitudes, φk are independent random phases uniformly distributed over [0; 2π[ (whose outcomes are fixed for a given experiment), and Ω = {2πk/N, k = 0, . . . , N − 1} is a set of frequency lines. Assuming the model plant of Fig. 1, the systematic error is then found to have the expression
Fig. 7. Noise assumption for the H2 estimator.
Interestingly, the variance due to leakage is now an O(N −3 ) for the coherence function instead of an O(N −1 ), as previously found for the FRF.
ˆ k) ≈ Bias G(ω
N −1 l=0
|Al |2 |W (ωl − ωk )|2 [G(ωl ) − G(ωk )] N −1 l=0
|Al |2 |W (ωl − ωk )|2 (28)
B. H2 Function The H2 function is the mean-square estimator of a FRF in the presence of input noise only, as illustrated in Fig. 7. Denoting by x(n) = u(n) + νx (n) the noisy measurement of the input signal, it is given by M −1 |Ywi (ωk )|2 ˆ G2 (ωk ) = M −1i=0 . (24) ∗ i=0 Xwi (ωk )Ywi (ωk ) Under Assumptions A1–A4 and provided that v(τ ) has a ˆ 2 (ωk ) due to bounded derivative, the bias and variance of G leakage are (see the Appendix) 2 2 C (1) (ω )| |G k ˆ 2 (ωk ) −−−−→ v ) G (ωk )+2 Bias G N τc 2N 2 G(ωk )∗ (25) 2 2 ˆ 2 (ωk ) = Cv (θ) |G (ωk )|2 + Dv (θ) · NSRin (ωk ), Var G M (θ)N 2 M (θ) ωk = mod(π)
(26)
where NSRin (ωk ) = S2νx (ωk )/S2u (ωk ) stands for the input NSR at frequency ωk , and where the correction factors Cv2 (θ) ˆ k ) (output-noiseand Dv2 (θ) are identical to those found for G(ω only case). C. Conclusion It is noteworthy that (21), (22), (25), and (26) all involve the same structure as (16) and (17). In conclusion, those settings ˆ k )—i.e., a half-sine or that were found to be optimal for G(ω a diff window with about 2/3 overlap—are also optimal for measuring the coherence function and the H2 function. VI. M ULTISINE E XCITATIONS So far, our analysis has focused on the case where the excitation signal is random stationary. A related question of interest is whether the same reasoning can be followed to find the optimal settings in the case of a multisine excitation. For that purpose, let us write the multisine excitation signal of period N as Ak ejωk n+φk (27) u(n) = u(n + N ) = ωk ∈Ω
with W (ω) as the DTF of the data window w(n). Clearly, the bias is exactly zero provided that W (ωl − ωk ) = 0 for any ωl = ωk , that is, provided that W (ω) has zeros at all frequencies ω = 2πk/N , k = 1, . . . , N − 1. For an N -long window, this can only be achieved by the rectangular shape. Similarly, the expression of the variance due to leakage noise can be found to depend on factors involving W (ωl − ωk ), so that those windows that cancel the systematic error also cancel the stochastic leakage errors. Finally, the variance due to measurement noise is found to be 2 ˆ k ) ≈ Fv (θ) · Bw · NSRout (ωk ) Var G(ω M (θ) ∼ θFv2 (θ) · O(L−1 )
(29)
where NSRout (ωk ) = S2νy (ωk )/|Ak |2 stands for the output NSR, Bw = N R2v (0)/|W (0)|2 stands for the bandwidth of the data window (e.g., Bw = 1/N for the rectangular window), and Fv2 (θ) = 1 + 2
M (θ)−1
i=1
1−
i M (θ)
κ(iθ)
(30)
is a positive correction factor very similar to (14) but without the power 2 on κ(iθ). Contrary to the previous cases, it can be checked that Fv2 (θ) is now minimum at θ = 1, that is, when there is no overlap. Moreover, the window shape that minimizes (29) is that with the minimum bandwidth Bw , which is again the rectangular window. This is illustrated in Fig. 8. It is also shown in this figure that overlapping the rectangular window has some unexpected consequences due to the “bouncing” behavior of Fv2 (θ): the worst choice is 1/4 overlap; suboptimal choices are 1/2, 2/3, 3/4, . . . overlap, but with an increased computational cost. This is summarized in Table III. In conclusion, the WOSA optimal settings for measuring FRFs in the case of multisine excitations with period N is a rectangular window of length N without overlap. Such a choice completely cancels the impact of leakage noise and minˆ k )} ≈ imizes that of measurement noise (in that case, Var{G(ω NSRout /L). As expected, this conclusion is fully coherent with that of the maximum-likelihood approach presented in [1].
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
ANTONI AND SCHOUKENS: OPTIMAL SETTINGS FOR MEASURING FRFs WITH WOSA
3283
Fig. 9. Simulated system. TABLE IV MEASUREMENT ERRORS AT THE RESONANCE FREQUENCY OF A ONE-DEGREE-OF-FREEDOM OSCILLATOR Fig. 8. Bw Fv2 (θ)/M(θ) as a function of the amount of overlap (1 − θ) for the (continuous line) rectangular and (dotted line) half-sine windows, with L = 10 000 and N = 100. Note that M (θ) = (L − N )/N θ + 1. TABLE III CORRECTION FACTORS FOR THE RECTANGULAR WINDOW
with i as the index of the experiment. In each case, the window length was set to N = 210 , which is about four times as great as the system time constant τc . A. Validation of Bias and Variance Formulas VII. N UMERICAL V ERIFICATIONS The results established in the preceding sections are now verified on simulations. The simulated system is a one-degreeof-freedom oscillator with the transfer function G(z) =
1 − 2z −1 + z −2 1 − 1.6079z −1 + 0.9875z −2
(31)
which shows a sharp resonance at ω0 = 0.2π. According to (9), its corresponding time constant τc is about 223 samples. The system was excited by a random Gaussian input of unit variance and of length L = 214 , as illustrated in Fig. 9. Note that significant bias and variance are expected at the resonance frequency ω0 according to (16) and (17), which involve the first and second derivatives of G(ω) (Table IV). To make the computation of experimental bias and variance possible, the same experiment was repeated 1000 times but with different (independent) inputs, so that, concerning the measurement of the FRF, for instance, its experimental bias and variance could be computed as follows: ˆ k) = G(ω Bias
ˆ k) = G(ω Var
1000 1 ˆ G(ωk ; i) − G(ωk ) 1000 i=1
2 1000 1000 1 ˆ 1 ˆ G(ωk ; i) G(ωk ; i) − 1000 i=1 1000 i=1
Our primary concern is to check the validity of the main results (16), (17), (21), (22), (25), and (26). No measurement noise was added to better emphasize the impact of leakage noise only on the systematic and stochastic errors. The parameters for the WOSA measurements were set with a half-sine window and 2/3 overlap. Fig. 10 compares the experimental bias and standard deviation with the corresponding formulas provided in this paper, for the coherence function, the classical FRF estimate, and the “H2 ” function. A very good match is generally noticeable. As expected from our formulas, the maximum errors occur in the vicinity of the resonance frequency where the leakage effect is the strongest. Concerning the FRFs, the bias error becomes so considerable in that frequency region that it overruns the stochastic errors (i.e., standard deviation; see Fig. 10(b) and (c)]. As for the coherence function [see Fig. 10(a)], its bias error largely dominates over the whole frequency range; this is because its variance was found to be very small: an O(N −3 ) instead of an O(N −1 ) for the FRFs [see (17), (22), and (26)]. Finally, it is noteworthy that the bias of the “H2 ” function suddenly drops at exactly ω0 . Again, this behavior can be well predicted by plugging (31) into (25), and it suggests that “H2 ” may be an excellent estimator of the resonance amplitude (this is true, in general, for any resonant system). B. Validation of the Optimal Settings Our second concern is to experimentally check the impact of overlap and of the window shape on the WOSA measurement.
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
3284
IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 58, NO. 9, SEPTEMBER 2009
Fig. 10. Comparison of (continuous lines) experimental and (dotted lines) theoretical biases and standard deviations in the case of a half-sine window ˆ 0 ). (c) G ˆ 2 (ω0 ). Bold lines depict the with 2/3 overlap. (a) γ ˆ 2 (ωk ). (b) G(ω actual quantities to be measured.
A first experiment was conducted to compare the estimation errors in the cases without overlap and with 2/3 overlap. No measurement noise was added, and a half-sine window was used. The experimental biases and standard deviations are shown in Fig. 11. It is obvious from these results that 2/3 overlap can reduce the standard deviation of stochastic errors by almost a factor of 2. Precisely, the case without overlap yields M (θ = 1) = 16 segments and, according to Table I, a correction factor Cv2 (1) = π 2 , whereas the case with 2/3 overlap yields M (θ = 1/3) = 46 and Cv2 (1/3) = 9.61(1 + 0.13/46); thus, the reduction in standard deviation is like the square root of Cv2 (1)/M (1) over Cv2 (1/3)/M (1/3), that is, 1.7. In addition, it is checked that overlap has no effect on the bias error, as predicted from (16), (21), and (25). A second experiment was conducted to obtain more quantitative results with respect to both overlap and the window shape. To summarize each simulation by a single couple of points, the experimental bias and variance results were integrated over the whole positive frequency axis, i.e.,
Bias2 =
L/2−1
2 ˆ k ) G(ω Bias
k=1 L/2−1
Var =
k=1
ˆ k) . G(ω Var
Fig. 11. Comparison of experimental biases and standard deviations (dotted lines) without overlap and (continuous lines) with 2/3 overlap. (a) γ ˆ 2 (ωk ). ˆ ˆ (b) G(ω0 ). (c) G2 (ω0 ). Bold lines depict the actual quantities to be measured.
Note that, according to these definitions, the total MSE is simply obtained as Bias2 + Var. The corresponding results are displayed in Fig. 12. As expected, it is shown that the squared bias is (statistically) constant against overlap, whereas the variance is a monotonically decreasing curve, with its minimum nearly reached at 2/3 overlap. This is independently of the window shape. These experimental curves fairly well reproduce the trends of Fig. 5, notwithstanding a different scaling factor relating to the system characteristics. Indeed, from (17), it theoretically holds that Var
L/2−1 Cv2 (θ) 2 |G (ωk )| M (θ)N 2 k=1
where Cv2 (θ)/M (θ) is the quantity displayed in Fig. 5. This again validates (17) and its implication about the optimality of the 2/3 rule. Finally, Fig. 12 demonstrates that the half-sine window is much better than the Hanning window as long as leakage errors (bias and variance) are of concern. A reduction of about 1.8 is observed on the squared bias, which is in perfect accordance with what can be predicted from (16) and Tables I and II, i.e., 2 2 C (1)Hanning 2 = 4 1.778. = v2 Cv (1)half -sine 3 Bias2 (half-sine) Bias2 (Hanning)
The reduction on the variance term is less significant, particularly at high levels of overlap. Equation (17) and Tables I and II,
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
ANTONI AND SCHOUKENS: OPTIMAL SETTINGS FOR MEASURING FRFs WITH WOSA
Fig. 12. Integrated squared bias and variance due to leakage with respect to the amount of overlap for the Hanning and half-sine windows.
in the case of no overlap, give Var(Hanning) C 2 (1)Hanning 4 = v2 = 1.33 Cv (1)half -sine 3 Var(half-sine) which fairly well compares with the reduction of about 1.4, as shown in Fig. 12, and in the case of 2/3 overlap, i.e., 0.27 10.50 1 + 2 M (1/3) Var(Hanning) Cv (1/3)Hanning = = Var(half-sine) Cv2 (1/3)half -sine 9.61 1 + 0.13 M (1/3)
1.10, which again is coherent with Fig. 12. Very similar results were obtained for the measurement of the coherence and the H2 functions, which are not displayed here due to the lack of space. C. Investigation of the Effect of Measurement Noise Our last concern is to check the scope of validity of our conclusions when the SNR becomes low. To do so, the total MSE Bias2 + Var was computed for different output SNRs (see Fig. 1), and the values were reported for both the Hanning and half-sine windows at various amounts of overlap. The experimental results are displayed in Fig. 13. It is shown that the half-sine window still outperforms the Hanning window, until the SNR reaches the unrealistically low value of −40 dB. At this stage, measurement noise completely overwhelms leakage noise, and the curves of Fig. 13(d) actually reproduce the trends of Fig. 6. This demonstrates that, at least for the very resonant system simulated in this section, leakage errors are very predominant (actually mainly the bias) that the half-sine window remains a very good choice, even in the presence of (realistic levels of) additive noise. VIII. C ONCLUSION The measurement of FRFs by means of the WOSA procedure suffers from three types of errors, namely, a systematic error
3285
Fig. 13. Overall MSE as a function of the SNR for the (· · · • · · ·) Hanning and (· · · ◦ · · ·) half-sine windows. (a) SNR = ∞ dB. (b) SNR = 0 dB. (c) SNR = −20 dB. (d) SNR = −40 dB.
due to leakage, a stochastic error due to leakage noise, and a stochastic error due to measurement noise. We have proven that the optimal real symmetric data-window shape that globally minimizes the leakage errors is the half-sine window, which incidentally achieves all the same properties as the (complexvalued and nonsymmetric) diff window recently proposed in [17]. Moreover, the stochastic errors can significantly be reduced by increasing the percentage of overlap between adjacent segments. The ultimate variance reduction that can achieved this way is about 3 on the leakage errors and 1.7 on the measurement noise errors with the half-sine window. An important practical result is that, in all cases, this is nearly achieved by setting about 2/3 overlap independently of the system FRF and of the data window. The “2/3 law” should replace the customary 1/2 law used in commercial data analyzers without increasing too much the computational demand. Moreover, it has been shown that the same recommendations apply as well to the measurement of other spectral quantities and, in particular, to the coherence function, which is often used in conjunction with FRF measurements. As a final note, we emphasize again that the optimality of these settings applies provided that the data-window length is set larger than the system time constant—which is in any case the usual recommendation with WOSA—and only with stationary random excitations. The special case of multisine excitations was verified to yield completely different conclusions. A PPENDIX The sketch of our proofs starts from the demonstration of relation (6). From Fig. 3, the leakage error due g(k)u(n − to segment i is twi (n) = y(n)w(n − iR) − ∞ k=0 k)w(n − k − iR), i.e.,
Twi (ωk ) =
Wi (ωk − ν) (G(ν) − G(ωk )) dU (ν)
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
(32)
3286
IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 58, NO. 9, SEPTEMBER 2009
where Wi (ωk − ν) = W (ωk − ν)e−jiR(ωk −ν) , and dU (ν) is the spectral increment of process u(n) at frequency ν [6]. Therefore, Ywi (ωk ) = G(ωk )Uwi (ωk ) + Twi (ωk ) + Nwi (ωk ), from which (6) immediately follows with M −1
i=0 βL (ωk ) = M −1 i=0
M −1
βN (ωk ) = i=0 M −1 i=0
Twi (ωk )Uwi (ωk )∗
(33)
Uwi (ωk )Uwi (ωk )∗
D. Equations (21) and (22)
∗
Nwi (ωk )Uwi (ωk )
Uwi (ωk )Uwi (ωk )∗
expected value of the denominator. Here, again, only the former quantity needs to be evaluated. From P2, P3 with p = 0, P4, P1, −1 ∗ and Assumptions A1–A3, Var{ M i=0 Nwi (ωk )Uwi (ωk ) } = 2 2 2 2 S2u (ωk )S2ν (ωk )(2π) M Dv (θ)R2v (0)N . Taking the ratio and using definition (12), (18) follows.
.
(34)
Now, from Assumption A4, Twi (ωk ) may be approximated as
1 2 Wi (ωk − ν) (ν − ωk )G (ν) − (ωk − ν) G (ν) dU (ν). 2 (35)
The evaluation of (21) and (22) starts with replacing Ywi (ωk ) in (20) by its expansion G(ωk )Uwi (ωk ) + Twi (ωk ). From Assumption A4, the magnitude of the transient term Twi (ωk ) is much smaller than Uwi (ωk ) on the average. Expanding ratio (20) to the second order with respect to Twi (ωk ) then yields the approximation M −1
2
γˆ (ωk ) 1 − i=0 M −1 i=0
All proofs proceed by making use of (33) and (34) and a set of properties: (P1): (P2): (P3): (P4):
E{dU (ν1 )dU (ν2 )} = S2u (ν1 )δ(ν1 + ν2 )dν1 dν2 . E{dU (ν1 )dN (ν2 )} = 0. (p) |W (ν)|2 ν p ≈ (−j)p 2πN (1−p) R2v (0), p = 0, 1, 2. The input signal u(n) is assumed to be Gaussian.
Property P1 follows from Assumption A1, Property P2 from Assumption A3, and Property P3 from the assumed differentiability of the window shape v(τ ). Property P4 is invoked here for simplicity, although not fundamentally necessary since the DFT makes (most) signal tend to Gaussianity in virtue of the central limit theorem [6].
From Assumption A4, the expected value of the ratio (33) tends to the ratio of the expected values. Namely, from P1, with p = 0, and Assumptions A1 and P3 −1 ∗ } = S2u (ωk )2πN R2v (0)M , A2, E{ M i=0 Uwi (ωk )Uwi (ωk ) −1 ∗ and from P3 with p = 2, E{ M i=0 Twi (ωk )Uwi (ωk ) } = −1 −S2u (ωk )πN R2v (0)M G (ω). Taking the ratio and using definition (12), (16) immediately follows. B. Equation (17) From Assumption A4, the variance of the ratio (33) tends to the ratio of the variance of the numerator to the square of the expected value of the denominator. Since the latter has just been evaluated in the preceding discussion, we only need to work out the former. Namely, from P1, P3−1with p = 0, 1, P4,∗ and Assumptions A1 and A2, Var{ M i=0 Twi (ωk )Uwi (ωk ) } = 2 2 (ωk )(2π)2 M Cv2 (θ)R2v (0)|G (ω)|2 . Taking the ratio and S2u using definition (12), (17) immediately follows.
C. Equation (18) From Assumption A4, the variance of the ratio (34) tends to the ratio of the variance of the numerator to the square of the
|Ywi (ωk )|2
.
(36)
The rest of the proof follows as for the bias and variance of ˆ k ) by taking the expected value and the variance of (36), G(ω simplifying the statistics of a ratio as the ratio of the statistics due to Assumption A4, and using Assumptions A1–A2 and Properties P1–P4.
E. Equations (25) and (26) The evaluation of (25) and (26) proceeds as for the H2 function; expansion of the ratio (24) to the second order with respect to Twi (ωk ) then yields the approximation
A. Equation (16)
|Twi (ωk )|2
G2 (ωk ) G(ωk ) ·
M −1
1+
∗ i=0 Twi (ωk )Ywi (ωk ) M −1 2 i=0 |Ywi (ωk )|
. (37)
The rest of the proof follows as before, by taking the expected value and the variance of (37), simplifying according to Assumption A1–A4, and using Properties P1–P4. R EFERENCES [1] R. Pintelon and J. Schoukens, System Identification—A Frequency Domain Approach. Piscataway, NJ: IEEE Press, 2001. [2] P. Welch, “The use of the fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms,” IEEE Trans. Audio Electroacoust., vol. AU-15, no. 2, pp. 70–73, Jun. 1967. [3] R. B. Blackman and J. W. Tukey, The Measurement of Power Spectra From the Point of View of Communications Engineering. New York: Dover, 1959. [4] G. M. Jenkins and D. G. Watts, Spectral Analysis and Its Applications. Boca Raton, FL: Emerson-Adams Press. [5] J. Bendat and A. Piersol, Random Data: Analysis and Measurement Procedures, 2nd ed. New York: Wiley-Interscience, 1986. [6] D. R. Brillinger, Time Series: Data Analysis and Theory. New York: SIAM, 1975. [7] A. H. Nuttall and G. C. Carter, “A generalized framework for power spectral estimation,” IEEE Trans. Acoust., Speech, Signal Process., vol. ASSP-28, no. 3, pp. 334–335, Jun. 1980. [8] A. H. Nuttall, “Spectral estimation of means of overlapped FFT processing of windowed data,” Naval Underwater Syst. Center, New London, CT, Rep. 4169. suppl. TR4169S.
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
ANTONI AND SCHOUKENS: OPTIMAL SETTINGS FOR MEASURING FRFs WITH WOSA
[9] A. H. Nuttall, “On the weighted overlapped segment averaging method for power spectral estimation,” Proc. IEEE, vol. 68, no. 10, pp. 1352–1354, Oct. 1980. [10] D. R. Brillinger, “The key role of tapering in spectrum estimation,” IEEE Trans. Acoust., Speech, Signal Process., vol. ASSP-29, no. 5, pp. 1075– 1076, Oct. 1981. [11] F. J. Harris, “On the use of windows for harmonic analysis with the discrete Fourier transform,” Proc. IEEE, vol. 66, no. 1, pp. 51–83, Jan. 1978. [12] A. Eberhard, “An optimal discrete window for the calculation of power spectra,” IEEE Trans. Audio Electroacoust., vol. AU-21, no. 1, pp. 37–43, Feb. 1973. [13] N. Geçkinli and D. Yavuz, “Some novel windows and a concise tutorial comparison of window families,” IEEE Trans. Acoust., Speech, Signal Process., vol. ASSP-26, no. 6, pp. 501–507, Dec. 1978. [14] A. Nuttall, “Some windows with very good sidelobe behavior,” IEEE Trans. Acosut., Speech, Signal Process., vol. ASSP-29, no. 6, pp. 84–91, Feb. 1981. [15] J. Le Roux and J. Ménez, “A cost minimization approach for optimal window design in spectral analysis of sampled signals,” IEEE Trans. Signal Process., vol. 40, no. 4, pp. 996–999, Apr. 1992. [16] J. L. Douce and L. Balmer, “Statistics of frequency-response estimates,” Proc. Inst. Elect. Eng., vol. 137, no. 5, pp. 290–296, Sep. 1990. [17] J. Schouken, Y. Rolain, and R. Pintelon, “Analysis of windowing/leakage effects in frequency response function measurements,” in Proc. 16th IFAC World Congr., Prague, Czech Republic, Jul. 4–8, 2005.
3287
Jérôme Antoni was born in Strasbourg, France, on November 11, 1972. He received the Ph.D. degree (with highest honors) in signal processing from the Polytechnic Institute of Grenoble, Grenoble, France, in 2000. He is currently a Lecturer with the University of Technology of Compiègne, Compiègne, France. His research interests include nonstationary signal analysis, system identification, and signal separation. His applications are concerned with noise and vibration analysis, modal analysis, and system diagnosis.
Johan Schoukens (M’90–SM’92–F’97) received the degree in engineering and the Ph.D. degree in applied sciences from the Vrije Universiteit Brussel (VUB), Brussels, Belgium, in 1980 and 1985, respectively. He is currently a Professor with the VUB. The prime factors of his research are in the field of system identification for linear and nonlinear systems. Dr. Schoukens was the recipient of the Best Paper Award and the Society Distinguished Service Award from the IEEE Instrumentation and Measurement Society in 2002 and 2003, respectively.
Authorized licensed use limited to: Rik Pintelon. Downloaded on September 29, 2009 at 03:22 from IEEE Xplore. Restrictions apply.
