Periodogram

Massachusetts Institute of Technology Department of Electrical Engineering and Computer Science 6.341: Discrete-Time Signal Processing OpenCourseWare 2006 Lecture 18 Periodogram Reading: Sections 10.6 and 10.7 in Oppenheim, Schafer & Buck (OSB). We begin this lecture by introducing three common illusions in spectral analysis: THREE ILLUSIONS • If you can’t see it, it’s not there. (the picket fence effect) • The more zero padding, the better the spectral resolution. (resolution vs. sampling, spectral smearing) • For a random process, as the data record length → ∞, the magnitude-squared of the DTFT converges to the power spectral density. (the periodogram) In the last lecture, we discussed the first and the second illusions. The picket fence effect re­ sults from the spectral sampling imposed by the DFT and can be avoided using zero padding. However, zero padding does not improve the spectral resolution which depends on the shape and length of the window. In this lecture, we will see the third illusion which relates to spectral analysis of stochastic signals. For a deterministic signal, more data (i.e. longer window) results in better frequency resolution, but this does not hold for a stochastic signal. The power spectral density(PSD) of a random process is useful for the following purposes: 1. System Identification x(t) - h(t) Pyy (Ω) = Pxx (Ω)|H (jΩ)|2 - y(t) Pyx (Ω) = Pxx (Ω)H (jΩ) From the PSDs of input and output, we can estimate the frequency response of the system. 1 2. Noise Removal Consider a signal: x(t) = s(t) + η(t), where s(t) is a deterministic signal, and η(t) is a stochastic noise process. In order to filter out the noise using a frequency-selective filter, we need to determine the PSD of the noise. The Periodogram Let xc (t) be a bandlimited stationary random process, then it can be sampled without aliasing. xc (t) - C/D x[n] - 6 T The power spectrum of x[n] is proportional to that of xc (t): Pxx (ω) = 1 Px x (Ω)|Ω=ω/T , T c c |ω| < π Therefore, we can estimate Pxc xc (Ω) from a reasonable estimate of Pxx (ω). The periodogram, defined as below, can be used as an estimate of the PSD of x[n]: I(ω) � 1 |V (ejω )|2 , L where V (ejω ) is the Fourier transform of the windowed sequence v[n] = RL [n]x[n]. The process of computing the periodogram of x[n] is illustrated in the figure below: w[n] = RL [n] x[n] ? -� × - v[n] DFT V [k] - | · |2 - |V [k]|2 = |V (ejω )|2 |ω=2πk/L The PSD is the Fourier transform of the autocorrelation function of the signal if the signal can be treated as a stationary random process. Using the properties of Fourier transform, it can be shown that I(ω) = = 1 DT F T {v[n] ∗ v[−n]} L L−|m|−1 � 1 v[n]v[n + m]}. DT F T { L n=0 2 Therefore, the periodogram is in fact the Fourier transform of the autocorrelation of the win­ dowed data sequence. (a) Periodogram of White Noise (b) Periodogram of Colored Noise Figure (a) above shows a white noise process and its periodogram using the 512-point DFT and linear interpolation. The PSD of the noise process is indicated as the flat line in the peri­ odogram figure. Notice that the periodogram has many deviations from the actual PSD. Figure (b) shows a colored noise process, its periodogram, and PSD. Although the periodogram looks very different from the actual PSD, it is apparent that the process has significant content only in low frequency. It is clear from the figure above that the periodogram is not very good estimate of the PSD in this example. We have already learned in the last lecture that increasing the size of the DFT does not improve the frequency resolution. However, in this case, increasing the length of the window is not helpful either, as discussed in the following section. Properties of the Periodogram An estimator is unbiased if its expectation is equal to the quantity that is being estimated. A consistent estimator is an estimator that converges to the quantity being estimated as the data size grows. We can determine whether the periodogram is biased and whether it is consistent by computing its mean and variance. As developed in OSB Section 10.6.2, E{I(ω)} = ∞ � 1 RL [n]RL [n + m]} Pxx (ω) ∗ DT F T { L n=−∞ 3 � where Pxx (ω) is the PSD of the signal. The term ∞ n=−∞ RL [n]RL [n+ m] is the autocorrelation of the rectangular window, thus it has a shape of a triangle. Denoting the Fourier transform of the window as W (ejω ), let B(ω) = DT F T { ∞ � RL [n]RL [n + m]} = |W (ejω )|2 n=−∞ then, 1 Pxx (ω) ∗ B(ω). L This can be interpreted similarly as in the deterministic case: the desired quantity is smeared by the spectrum of the window (|W (jω)|2 in this case). Since E{I(ω)} is not equal to Pxx (ω), we see that the periodogram is a biased estimate of the power spectrum. E{I(ω)} = As L goes to infinity, W (ejω ) approaches an impulse, and thus B(ω) also approaches an impulse at the origin. In this case, E{I(ω)} ≈ Pxx (ω), so the periodogram is an asymptotically unbiased estimator. It has been shown that as the window length increases, 2 var{I(ω)} ≈ Pxx (ω). Therefore, the variance does not approach zero as L → ∞, and the periodogram is not a con­ sistent estimate of the power spectrum density. OSB Figure 10.20 shows the periodograms of a white-noise sequence with variance 1. The correct PSD for this process is a constant of 1. We can see that as the window length in­ creases from L = 16 in (a) to L = 1024 in (d), the variation from the mean does not go to zero, and the periodogram with longer window does not give a better estimate of the power spectrum. Periodogram Averaging In order to reduce the fluctuations and obtain a smooth spectrum estimate, we can average multiple measurements of periodogram estimates. Let x[n] be an ergodic random signal, then the expectation can be calculated by time averaging. Assume that we want to estimate the mean defined as follows: +∞ 1 � x[n] N →∞ 2N −∞ E{x[n]} = mx = lim Consider using the following estimator: K 1 � m̂x = x[n] K 1 4 Then, m ˆ x is an unbiased estimator since E{m ˆ x } = mx . If we assume that all the observation of x[n] are independent, and denote the variance of x[n] as σx2 , then var{m ˆ x } = K1 σx2 . Therefore, var{m ˆ x } → 0 as K → ∞, and m ˆ x is a consistent estimator. Similarly, we can construct a consistent estimator of power spectrum utilizing an arbitrarily long data record. If we have multiple measurements of I(ω), we can think of averaging them: K 1 � ¯ Ir (ω) Ixx (ω) = K r=1 Then, E{I¯xx (ω)} = E{Ir (ω)} = 1 Pxx (ω) ∗ B(ω) L and 1 1 2 (ω), var{Ir (ω)} = Pxx K K so as K → ∞, I¯xx (ω) converges to L1 Pxx (ω) ∗ B(ω). var{I¯xx (ω)} = Now, consider a fixed data record of length Q. If there is no overlap, Q = KL, where L is the length of the window, and K is the number of measurements of I(ω). We get more accurate estimate of I(ω) if K is larger, but we also need a longer window in order to increase spectral resolution. So, there is a tradeoff between K and L as illustrated in OSB Section 10.6.5 and summarized in the following example. Example: Consider the sequence x[n] = A cos(ω0 n + θ) + e[n], where θ is a random variable uniformly distributed between 0 and 2π, and e[n] is a zeromean white-noise sequence that has a constant power spectrum Pee (ω) = σe2 for all ω. It can be shown that E{I¯xx (ω)} = A2 L + σe2 . 4 OSB Figure 10.23 shows average periodogram estimates for A = 0.5, σe2 = 1, and with different values of L and K. Notice that as K increases (i.e. more sections are averaged), the periodogram becomes smoother because the variance of the estimate decreases. However, since the data length is fixed, the window length should become shorter to average more sections. As a result, the frequency resolution decreases because of the smearing effect, and in OSB Figure 10.23(d), the spectral peak due to the cosine is very broad and barely above the noise. 5 In OSB Figure 10.23(b), we used overlapped sections when averaging periodograms. The variance is reduced by almost a factor of 2 when the overlap is one-half the window length. Our last example illustrates periodogram averaging using a non-rectangular window and overlapping sections. Example: Consider the two sequences shown below: a noisy sequence and a ”clean” sequence after filtering noise using an elliptic filter. If we apply a Fourier transform to the entire sequences using the FFT (covered in the next lecture), the resulting figures are hard to interpret because of fluctuations. 6 The figures below show averaged periodograms using 4096- and 1029-point Hamming win­ dows, respectively. In both cases, neighboring segments are overlapped by one-half the window length. Notice that shorter window results in increased bias, but the estimate is smoother because more segments are averaged. Window: 4096 point, Hamming, Overlap: 50 % Window: 1024 point, Hamming, Overlap: 50 % 7