Journal ArticleThis paper introduces novel applications of the short-time unbiased sprectral estimation (STUSE) algorithm, which adds biased estimates to yield unbiased spectral estimates. It is shown that STUSE algorithm is an effective tool for estimating time delays and magnitude-squared coherence (MSC) functions between two stationary signals received at spatially separated sensors, especially when one of the signals is a delayed version of the other. Computer simulation results are presented to compare the performances of the STUSE algorithm and the conventional weighted overlapped segment averaging (WOSA) method for spectrum estimation

Mathews, V. John

Youn, D. H.

English

The University of Utah: J. Willard Marriott Digital Library

ON USING THE SHORT-TIME UNBIASED SPECTRAL ESTIMATION ALGORITHM FOR ESTIMATING TIME DELAYS AND MAGNITUDE SQUARED COHERENCE FUNCTIONS D.H. YOUN and V.J. MATHEWS Department of Electrical and Computer Engineering The University of Iowa Iowa City, Iowa 52242 ABSTRACT This paper introduces novel applications of the short-time unbiased sprectral estimation (STUSE) algorithm, which adds biased estimates to yield unbiased spectral estimates. It is shown that STUSE algorithm is an effective tool for estimating time delays and magnitude-squared coherence (MSC) functions between two stationary signals received at spatially separated sensors, especially when one of the signals is a delayed version of the other. Computer simulation results are presented to compare the performances of the STUSE algorithm and the conventional weighted overlapped segment averaging (WOSA) method for spectrum estimation. 1. INTRODUCTION The purpose of this paper is to introduce the use of the STUSE algorithm [1,2] for estimating time delays and MSC functions between stationary signals, received at spatially separated sensors. Estimating the time delay between the arrival times of the same signal at spatially separated sensors has applications in sonar, radar, acoustics, biomedical engineering and other areas [9]. The MSC function between two signals is a measure of linearity between them [10] and finds applications in time delay estimation (TDE) [3,11], system analysis [8], and measuring signal-to-noise ratio [4]. A common approach to time delay estimation is to find the time lag at which the cross-correlation function between the received signals x 1(k) and x2(k) is maximum [3]. The cross-correlation function, (la) is computed as (lb) where E{o} denotes statistical expectation of {o}, F-1{o} denotes the inverse Fourier transform (1FT) of {o}, and G12(f) is the cross-power density spectrum (cross-PDS) of xl(k) and x2(k). CHI944-8/83/0000-0060 $1.00 © 1983 IEEE 60 The MSC function of two signals xl(k) and x 2(k) is defined as [4] where G11 (f) and G22 (f) are the auto-PDS's of x1(k) and x2(k), respectively. (2) From (1b) and (2) it can be seen that computing the auto- and cross-PDS's plays a key role in measuring time delays and MSC functions. One of the widely used spectral estimation method is the WOSA method [5,6], in which the cross-PDS between two signals xl(k) and x2(k) is computed as where carets denote estimated quantities, * denotes complex conjugates, and ~ n(f) represents the 2M-point discrete Fourier transform (DFT) of the n-th weighted segment of xi(k) defined as Xi,n(f) = F{xi(nR + k)w(k)} (3b) for 0 ~ k ~ L - 1, 0 ~ n ~ N - 1, i = 1 or 2, and M ~ L, where w(k) is a window function of length L, R denotes the number of samples between successive segments (i.e., L - R samples are overlapped) and N is the number of segments. In (3a) rww(O) is the autocorrelation function of the window function at zero lag. Here rww(m) is defined as L-l rww(m) - 2M L w(k) w(k-m); - M ~ m ~ M - 1 k=O (3c) It has been shown [1,2,6] that the expected values of G12 (f) and C12 (m) computed using the WOSA algoriffim are given by ( 4a) and (4b) where 0 denotes complex convolution and Gww( f) is the auto-PDS of the window function. From (4a) and (4b) we, can see that, for G12 (f) and ~12(m) to be unbiased IW(f)1 2 should be a delta function. However, for all finite length linear windows, IW(f)1 2 exhibit low-pass characteristics and the estimates in (4a) and (4b) are biased. The implications of these biases are two-fold; in the frequently domain, convolution of the cross-PDS with the auto-PDS of the window results in reduced spectral resolution, and in the time domain, multiplication of the cross-correlation function with the autocorrelation function of the window results in a correlation function estimate which is gradually tapered to zero as the lag approaches the window length in either direction. Similarly, when one of the signals is a delayed version of the other, it has been shown [12] that the MSC function estimates via WOSA method are biased. This bias, due to misalignment of data or rapidly changing phase, has been shown to be [7] B(0,;{-~+(Q.)2}h (01 2 (5) Y L L 12 where D is the time delay between two signals. The bias in (5) can be reduced by realigning the data [7] to yield zero delay (i.e., D = 0) and/or increasing the analysis window length L. Recently, a new spectrum estimation technique has been proposed, in which the influences of the finite window length on a spectral estimate can be removed by linearly combining biased estimates [1,2]. This method has been referred to as the STUSE algorithm [1]. After discussing the STUSE algorithm in the next section, the performance of the STUSE algorithm for estimating the MSC and cross correlation functions will be demonstrated in Section 3 via computer simulations. 2. STUSE Algorithm The STUSE algorithm for computing the cross-PDS of xl(k) and x2(k) can be summarized as follows; G12(0 - N ___ ~2 [N~l rww(O) q=ql n=O ] . 21T f R * J N q Xl,n(f) X2,n+q(f) e ( 6a) where ~ (m) ww (6b) In (6a) the exponential term takes care of the time delay between the n-th and (n+q)-th segments. In (6a), it is assumed that the cross-correlation function of xl(k) and x2(k) is zero for Iml > M. i.e., 61 We will later on see, how the algorithm can be modified if (7) does not hold. Comparing (6a) and (3a), we can see that the STUSE algorithm becomes the WOSA method if ql = q2 = O. As for the case of the WOSA method, we can show that the expected values of the cross-correlation function and cross-PDS estimates using the STUSE algorithm are given by G (0 E{G12 (0} G (00~ ( Sa) 12 r (0) and ww r (m) E{c 12 (m)} = c 12(m) • ww Iml .. M (Sb) ~ (0) ww where G (f) denotes the Fourier Transform ww of r (m) in (6b). Thus, from (6b) and (Sb) we can ~e that, if ql' q2 ~nd R are carefully chosen to have constant r (m) for Iml .. M, the A ww estimates G12(f) in (6a) a~d the corresponding crosscorrelation function c 12(m) obtained using the STUSE algorithm are unblased. It has been shown [1] that any band limited window function should be applied at or above the Nyquist decimation rate of the window to yield unbiased estimates. For example, 0 < R .. L/4 for Hamming window [1]. Also, for the estimated cross correlation function to be unbiased in -M .. m .. M, we can see that ql and q2 should be chosen such that - q 1 = q 2 = l2~ J ( 9 ) where ~.~ denotes the largest integer smaller than (0). If the above conditions are met, we have r (m) and ww ~ (0) ww Iml .. M (lOa) G (0 r (0) 6(0 (lOb) whereww ww 1 f = 0 6(f) = {O: elsewhere , (lOc) thus yielding unbiased estimates of cross-PDS and crosscorrelation function, as can be seen by substituting (lOb) and (lOa) in (Sa) and (Sb), respectively. If (7) is not satisfied, the algorithm should be modified in the following manner. For q = ql' ••• , q2' let us define G (q) (0 12 - N N-l * ~ Xl,n(f) X2,n+q(f), n=O and Now compute G12(f) as (lla) (lIb) (l2) Note that Iml in (12) may be larger than M and thus the FFT length in (lIb) may be more than 2M. In (12), the time shift qR corresponds to the exponential term in (6a). In the next section, we demonstrate the effectiveness of the STUSE algorithm for TDE and MSC function estimates via computer simulation results. 3. SIMULATION RESULTS Time Delay Estimation TDE-1 and Consider a multipath TDE problem defined by 3 1: s(k - D.) i=-3 1 x 2(k) = s(k), for k = 0,1, ••• ,499 (13) where s(k) is zero mean white Gaussian signal with unit variance. There are seven delay parameters, Di = 10i due to seven different paths, to be estimated here. Figure 1a displays a crosscorrelation function estimate via the WOSA method, where 75% overlapping and a 32 sample Hamming window were used and a 64 point FFT was taken for each segment. The sequence of arrows in Fig. 1a is the theoretical cross correlation function between x1(k) and x2(k). From Fig. 1a we can see that the WOSA method yields biased estimates of the cross correlation function, missing altogether the peaks at ±30 samples delay and almost missing the ones for ±20 samples delay. Figure 1b shows a cross correlation function estimate using the STUSE algorithm, where R = 8, -q1 = q2 = 3 and the other parameters are the same as those used for the WOSA method. Comparing Figs. 1a and 1b it can be observed that using the STUSE algorithm significantly reduces the bias in the estimated crosscorrelation function. TDE-2 To demonstrate the use of the ST~SE algorithm when (7) is not satisfied, consider the following set of signals; s(k + 40) + s(k + 10) + s(k) + s(k - 10) + s(k - 45) s(k), for k = 0,1, ••• ,499 (14) where s(k) is, as before, white Gaussian signal with zero mean and unit variance. Figure 2a displays the crosscorrelation function between x1(k) and x2(k) obtained using (6a) where a 32 point Hamming window and 64 point FFT's were used with -q1 = q2 = 3 and R = 16; i.e., taking care of the shifts in the frequency domain. It can be seen that the crosscorrelation function in Fig. 2a has spurious peaks at -24 and +19 sample delays. These peaks are due to the wrapping around of the cross correlation function for time delays greater than 32 samples. 62 Figure 2b displays the crosscorrelation function estimate of x1(k) and x2(k) obtained using (12). All the parameters used were the same as before. From Fig. 2b it can be seen that one can obtain correct estimates of cross-correlation functions using the STUSE algorithm, even when the cross correlation functions are nonzero for time lags greater than the window lengths in either direction. A Figures 2(c) - (i) show plots of c1~q) (m - qR); Iml ~ 50 for q = 3, 2, ••• , -3 respectively. The sequence of plots heuristically demonstrates how the STUSE algorithm obtains unbiased crosscorrelation function estimates, using a sequence of biased estimates. MSC Function Estimation The set up in Fig. 3 was used to generate signal sequences x1(k) and x2(k) used to compute IY12(f)12. In Fig. 3, s(k) is white Gaussian signal with zero mean and unit variance. Two different simulation examples were run. The results presented in each case are the averages over 60 independent runs using 500 data points each. MSC-1 H1(z) = H2(z) = 1 and w1(k) and w2(k) are mutually independent additive white Gaussian noises with zero mean and unit variance. The actual and estimated MSC functions using the STUSE and WOSA algorithms are displayed in Fig. 4a where 32 point Hamming window, 64 point FFT's, 75% overlap and -q2 = q1 = 2 were used. It can be seen that the WOSA algorithm yields biased estimates, while almost unbiased MSC function estimates were obtained using the STUSE algorithm using only 5 shifts (different values of q) when one of the signals is a delayed version of the other. Figure 4b shows plots of the quality ratio of the estimates, defined as (15) where var (.) and av (.) represent the variance and average value of (.), respectively. Figure 4b suggests that the MSC function estimates using the STUSE algorithm are more stable than those obtained using the WOSA algorithm, at least when the delays involved are reasonably large. MSC-2 This is the same example that was considered in [8] for comparing seismometers with MSC functions. H1(z) 1, -1 z H2 (z) = -----:-----;;-1 - z-l + 0.8z-2 and w2(k) = 0 for all k and w1(k) was white Gaussian noise with zero mean and unit variance. The other parameters were the same as those used for the last experiment. Figure Sa displays the true MSC function along with the estimates using the STUSE and WOSA algorithms. The quality ratios of the two estimates are plotted in Fig. Sb. Once again, it can be seen that the estimate using the STUSE algorithm is less biased and has a better quality ratio than that using the WOSA algorithm. 4. CONCLUSIONS Preliminary results of using the STUSE algorithm for estimating time delays and MSC functions between two stationary signals are presented. The effectiveness of the algorithm was demonstrated specifically (1) when the delay parameter to be estimated is relatively large compared to the window length and (2) when one of the signals, whose MSC function is to be estimated, is a delayed (possibly noisy) version of the other. It was pointed out that to use (6a), the crosscorrelation function between xl(k) and x2(k) should be zero for lags of magnitude greater than half the FFT length. If this condition does not hold, or if one does not have prior knowledge of the cross correlation function, equations (10a)-(12) should be used in computing the PDS using the STUSE algorithm. Simulation results show that when large delays are involved, the quality ratio of the MSC function estimates using the STUSE algorithm is smaller than those using the WOSA algorithm. However, this will not be the case if the delays are small, since the improvement in quality ratio was mainly contributed by reduction in the bias due to misalignment. While it may be necessary to use large values of q (shift) to obtain completely unbiased estimates, it was demonstrated that, smaller values of q will provide considerably unbiased estimates. Further theoretical aspects of the STUSE method are being investigated. REFERENCES 1. L.R. Rabiner and J.B. Allen, "Short-time Fourier analysis techniques for FIR system identification and power spectrum estimation," IEEE Trans. Acoustics, Speech, and Signal Proc., Vol. ASSP-27, No.2, pp. 182-192, April 1979. 2. , "On the implementation of a short-time spectral analysis method for system identification", IEEE Trans. Acoustics, Speech, and Signal Proc., Vol. ASSP-28, No.1, pp. 69-78, Feb. 1980. 3. C. H. Knapp and G.C. Carter, "The generalized correlation method for estimation of time delay", IEEE Trans. Acoustics, Speech and Signal Proc., Vol. ASSP-24, pp. 320-327, Aug. 1976. 4. G.C. Carter, C.H. Knapp, and A.H. Nuttal, "Estimation of magnitude-squared coherence function via overlapped fast Fourier transform processing", IEEE Trans. Audio Electrocoust., Vol. AU-21, pp. 337-344, Aug.19 73. 5. P.D. Welch, "The use of FFT for estimation of power spectra: A method based on time averaging of short, modified periodograms", IEEE Trans. Audio Electroacoust., Vol. AU-IS, No.2, pp. 70-73, June 1967. 6. A.H. Nuttal, "Spectral estimation by means of overlapped FFT processing of windowed data", New London, CT, NUSC Rep. 4169, Oct. 13, 1971. 7. G.C. Carter, "Bias in magnitude-squared coherence estimation due to misalignment", IEEE Trans. Acoustics, Speech, and Signal Proc., Vol. ASSP-28, No.1, pp. 97-99, Feb. 1980. 8. S.D. Stearns, Applications of the Coherence Function in Comparing Seismometers, Sandia National Laboratories, SAND79-1633, Dec. 1979. 9. G. C • Carter, Ed. ,Special Issue Time Delay Estimation, IEEE Trans. Acoustics, Speech, Signal Proc., Vol. ASSP-29, No.3, June 1981. 10. G.C. Carter, "Estimation of the magnitude-squared coherence function", New London CT, NUSC Rep. 4343 May 19, 1972. 11. J.S. Bendat and A.G. Piersol, Engineering Applications of Correlation and Spectral Analysis, New York: Wiley, 1980. 12. L.H. Koopmans, The Spectral Analysis of Time Series, New York: Academic, 1974. 13. F.J. Harris,: "On the use of merged, overlapped, and windowed FFT's to generate synthetic time series data with a specified power spectrum", Proc. Sixteenth Asilomar Conference on Circuits, Systems and Computers, Nov. 8-11, Pacific Grove, California, pp. 316-321. 63 t -,uj (a) "i I ~J I -32 0 32 (b) Fig. 1. Cross correlation function estimates for example TDE-1 via m (a) WOSA method and (b) STUSE algorithm q~ 2 e I I ,~~, A, J q=- 3 ~ L.c-. -50 50 Fig. 2. (a) - (b): Cross correlation function estimate for example TDE-2 (a) using (6a) and (b) using (12). (c) - (i) : Plots of ~i~)(m-qR) for several values of q I (,2( I) 12 05r---------------------~--------------~~ --~-------~-~-~------,~~--~--~-Fig. 4a. MSC function estimates for example MSC-l. (1) True MSC function (2) Estimate using STUSE algorithm (3) Estimate using WOSA method. Q(!) Fig. 4b. Quality ratio of the estimates in example MSC-l using (1) STUSE algorithm and (2) WOSA method 64 w1(k) 1----... ~)----40) Xl (k) ----' s(k) t---~(+}---~ x2 (k) ----' Fig. 3. Block diagram for generating signals used for MSC function estimation NORMALIZED FREQUENCY 1.0 Fig. 5a. MSC function estimates for example MSC-2. (1) True MSC function (2) Estimate usi~!o STUSE algorithm (3) Est imate using WOSA method. Q(f) 1·5;-· ----------------------------, NORMALIZED FREQUENCY Fig. 5b. Quality ratio of the estimates in example MSC-2 using (1) STUSE algorithm and (2) using WOSA method 

On using the short-time unbiased spectrum estimation algorithm for estimating time delays and magnitude squared coherence functions

https://collections.lib.utah.edu/dl_files/67/87/6787c8353b4b51d5e99a7a597fb45ad13f3d92b6.pdf

On using the short-time unbiased spectrum estimation algorithm for estimating time delays and magnitude squared coherence functions

Abstract

Similar works

Full text

Available Versions

The University of Utah: J. Willard Marriott Digital Library