High resolution tools for spectral analysis

Research supported by AFSOR, NSF. Maintained by L. Ning and T.T. Georgiou

Introduction

Consider a discrete-time zero-mean stationary process ${ldots, y_0, y_1, {ldots}, y_{N-1}, ldots}$ , $y_tin mathbb{C}$ , and let

$r(k)=E(y_t y^*_{t-k}), kin mathbb{Z}$

be the covariance lags. The power spectral density (PSD) function

$f(theta)=sum_{k=-infty}^{infty} r(k) e^{-itheta k},; thetain [pi,; -pi),$

represents power density over frequency. Under suitable smoothness conditions, it is also

$f(theta)= lim_{Nrightarrowinfty} E left{frac{1}{N} |sum_{t=0}^{N-1} y(t)e^{-itheta t}|^2 right}.$

Spectral estimation refers to esimating f(theta) when only a finite-sized observation record ${y_0, y_1, {ldots}, y_{N-1}},; N<infty$ of the time series is available. Different schools of thought have evolved over the years based on varying assumptions and formalisms. Classical methods began first based on Fourier transform techniques and the periodogram, followed by the so called modern spectral estimation methods such as the Burg method, MUSIC, ESPRIT, etc. The mathematics relates to the theory of the Szeg {ddot{o}} orthogonal polynomials and the structure of Toeplitz covariances.

In contrast, recently, the analysis of state covariance matrices, see e.g. [1, 2], suggested a general framework which allows even higher resolution. A series of generalized spectral estimation tools have been developed generalizing Burg, Capon, MUSIC, ESPRIT, etc.. In the following, we will overview some of these high resolution methods and the relevant computer software. Their use is shown via a representative academic examples.

We want to resolve two sinusoids in based on N=100 noisy measurements of the time series

$y_{k}=a_1 e^{itheta_1k+phi_1}+a_2 e^{itheta_2k+phi_2}+mbox{noise}_k, mbox{~for~} k=0,ldots,{N-1},$

where theta_1=1.3 and theta_2=1.35 . These are closer than the theoretical resolution limit 2pi/N of Fourier method.

code1

%  SIGNAL = sinusoids + noise
%  Setting up the signal parameters and time history
    N=100;
    a0=1.8; a1=1.5; theta1=1.3; a2=2; theta2=1.35;
    k=1:N; k=k(:);
    y=a0*randn(N,1)+a1*exp(1i*(theta1*k+2*pi*rand))+a2*exp(1i*(theta2*k+2*pi*rand));

FFT-based methods

We begin with two classical formulae. The correlogram is defined as

$hat{f}_c =sum_{k=-(N-1)}^{N-1} hat{r}(k) e^{-itheta k},$

where hat{r}(k) denotes an estimate of the covariance lag r(k) , whereas the periodogram

$hat{f}_p(theta) = frac{1}{N} |sum_{t=0}^{N-1} y(t) e^{-itheta t}|^2.$

The two coincide when the correlation lag hat{r}(k) is estimated via

$hat{r}_k=frac{1}{N} sum_{t=-N+1}^{N-1}y(t)y_b(k-t):= frac{1}{N} (y*y_{b})(k),; k=-(N-1), ldots, 0, ldots, N-1,$

where y_b(t):=y^*(-t) . The periodogram can be efficiently computed using the fast Fourier transform (FFT). There is a variety of methods, such as Welch and Blackman-Tukey methods, designed to improve the performance using lag window functions either in the time domain or in the correlation domain. In situations when the data length is short, to get a smooth spectrum, we may increase the data length by padding zeros to the sequence.

Using the above example, we pad up 2048 by adding zeros. The corresponding FFT-based spectrum is shown in the following figures. Rudimentary code is displayed as a demonstration. A file to reproduce the following results can be downloaded. Alternative matlab build-in routines for periodograms are periodogram, pwelch, etc.

fftbased

% the fft-based spectra
NN=2048; th=linspace(0,2*pi,NN);
Y =abs(fft(y,NN))/sqrt(N);
Y = Y.^2;
plot(th,Y);

Input-to-state filter and state covariance

The first step to explain the high resolution spectral analysis tools is to consider the input-to-state filter below and the corresponding the state statistics. The process ${y_k}$ is the input and x_k is the state.

Then the filter transfer function is

$G(z)=(I-z^{-1}A)^{-1}B.$

A positive semi-definite matrix is the state covariance of the above filter, i.e. R:=E(x_kx_k^*) for some stationary input process ${y_t}$ , if and only if it satisfies the following equation

for some row-vector . Starting from a finite number of samples, we denote

$hat{R}=frac{1}{N} sum_{k=1}^N x_k x_k^*$

the sample state covariance matrix.

If the matrix is the companion matrix

A = left[ begin{array}{ccccc} 0 & 0 & ldots & 0 & 0 1 & 0 & ldots & 0 & 0 vdots & vdots & & vdots& vdots 0 & 0 & ldots & 1 & 0 end{array} right], ~~ B= left[ begin{array}{c} 1 0 vdots 0 end{array} right], hspace*{3cm} (1)

the filter is a tapped delay line:

The corresponding state covariance is Toeplitz matrix. Thus the state-covariance formalism subsumes the theory of Toeplitz matrices.

The input-to-state filter works as a “magnifying glass” or, as type of bandpass filter, amplifying the harmonics in a particular frequency interval. Shaping of the filter is accomplished via selection of the eigenvalues of . In the above example, we set n=5 eigenvalues of at $0.88e^{1.325i}$ , and the phase angle 1.325 is the middle of the interval [theta_1,;theta_2] where resolution is desirable. Then we build the pair A, B using routine cjordan.m, (and rjordan.m for real valued problem). The pair is normalized to satisfy

The routine dlsim_complex.m generate the state covariance matrix (dlsim_real.m for real valued problem). The following figure plots $|G(e^{itheta})|$ versus theta . The gain $|G(e^{itheta})|$ at theta equal 1.2 and 1.45 is approximately dB below the peak value. The window is marked in the following figure. The figure to the right shows the detail within the window of interest.

code2

% setting up filter parameters and the svd of the input-to-state response
thetamid=1.325;
[A,B]=cjordan(5,0.88*exp(thetamid*1i));
% obtaining state statistics
R=dlsim_complex(A,B,y');
sv=Rsigma(A,B,th);
plot(th(:),sv(:));

Maximum entropy spectrum

The entropy of a probability distribution function p(x) with $xin mathbb{R}^n$ ,

quantifies the uncertainty of the corresponding random variable. When p(x) is a multi-dimensional Gaussian distribution with zero mean and covariance matrix , the entropy becomes

$frac{1}{2}log |M|+frac{n}{2}-frac{n}{2}log (2pi).$

For a zero-mean stationary Gaussian process {y(t), tin mathbb{Z}} , corresponding to infinite-sized Toeplitz structured covariance, the entropy “diverges”. Thus, one uses instead the ‘‘entropy rate". Let $T_{N-1}$ be a dimensional principle sub-matrix of the covariance matrix

$T_{N-1}=left[ begin{array}{cccc} r(0) & r(1) & ldots & r(N-1) r(-1) & r(0) & ldots & r(N-2) vdots & vdots & ddots & cdots r(-N+1) & r(-N+2) & r(-1) & r(0) end{array} right],$

Then the entropy rate is

$lim_{Nrightarrow infty} frac{1}{N} left( frac{1}{2}log |T_{N-1}|+frac{N}{2}-frac{N}{2}log (2pi) right)= lim_{Nrightarrow infty} frac{1}{2} log |T_{N-1}|^{1/N}+frac{1}{2}-frac{1}{2}log (2pi).$

The limit of $|T_{N-1}|^{1/N}$ converges to the optimal one-step prediction error variance, and from the Szeg-Kolmogorov formula

$lim_{Nrightarrow infty} |T_{N-1}|^{1/N} = frac{1}{2pi} e^{ int_{-pi}^{pi} f(theta)frac{dtheta}{2pi}}.$

Given state statistics, the maximum entropy power spectrum is

$f_{rm me}(theta)=argmax_{f(theta)} left{ int_{-pi}^{pi} log f(theta) dtheta ~mid~ int_{-pi}^{pi} G(e^{itheta}) f(theta)G(e^{itheta})^*dtheta = R right}.$

The solution to this problem (see [2]) is

$f_{rm me}(theta)= P(e^{itheta})^{-1}Omega (P(e^{itheta})^{-1})^*$

where

$begin{array}{rl}Gamma:=& (B^*R^{-1} B)^{-1} B^* R^{-1} Omega:=& Gamma RGamma^* P(z):=& Gamma G(z) end{array}$

and G(z) is the input-to-state filter defined above. This formula subsumes the classical Burg method/AR modeling where the is a Toeplitz matrix and G(z) is tapped delay line filter.

The maximum entropy spectrum is obtained using the routine me.m. For the example discussed above, the maximum entropy spectrum is shown in blue. There are two peaks detected inside the window. Burg's spectrum is shown in green. The resolution of Burg's solution is not sufficient to distinguish the two peaks.

mespectrum

% maximum entropy spectrum and Burg spectrum
me_spectrum=me(R,A,B,th);
plot(th,me_spectrum);
me_burg=pburg(y,5,th);
hold on;
plot(th,me_burg);

Subspace methods

We start with the Carathodory-Fejr-Pisarenko result for Toeplitz matrices. Given a ntimes n positive definite Toeplitz matrix , let c>0 be the smallest number such that T-cI is singular and has rank m<n , then T-cI has the formE

$T-cI=sum_{k=1}^m rho_k F(e^{itheta_k})F(e^{itheta_k})^*$

where

$F(e^{itheta_k})= left[ begin{array}{c} 1 e^{itheta_k} vdots e^{i(n-1)theta_k} end{array} right].$

Accordingly, the power spectrum f(theta) decomposes as

$f(theta)= sum_{k=1}^m 2pirho_k delta(theta-theta_k)+c$

where delta(cdot) is the Dirac function. The decomposition corresponds to a set of white noise. See MA decomposition for a decomposition corresponds moving-average (MA) noise.

The above result can be generalized to the case of state covariances [1]. More specifically, let be the unique solution to the Lyapunov equation

The matrix is the state covariance when the input is pure white noise. Let now be an arbitrary state covariance, let be the smallest eigenvalues of the matrix pencil R-rho Q , and assume that the rank of R-c Q is , then

$R-c Q=sum_{k=1}^m rho_k G(e^{itheta_k}) G(e^{itheta_k})^*,$

where $G(e^{itheta})$ is generalization of $F(e^{itheta})$ .

The subspace spectral analysis methods rely on the singular value decomposition

$R-cQ= U mbox{diag}{lambda_1, lambda_2, ldots, lambda_m, 0,ldots, 0 } U^*,$

where is a unitary matrix and lambda_k>0 , k=1,ldots, m . Partition

$U=[U_{1:m},; U_{m+1:n}]$

where $U_{1:m}$ and $U_{m+1:n}$ are the first and the last n-m columns of , respectively. Based on this decomposition, there are two ways we can proceed generalizing the MUSIC and ESPRIT methods, respectively [P. Stoica, R.L. Moses, 1997].

Noise subspace analysis

The columns of $U_{m+1:n}$ span the null space of R-cQ , while the signal $G(e^{itheta_k})$ , k=1,ldots, m is in the span of the columns of $U_{1:m}$ . So the nonnegative trigonometric polynomial

$d(e^{itheta},e^{-itheta})=G(e^{itheta})^* U_{m+1:n}U_{m+1:n}^*G(e^{itheta})$

has roots at ${theta_1, ldots, theta_m }$ .

Given a sample state covariance matrix hat{R} and an estimate on the number of signals m<n , we let $hat{U}_{m+1:n}$ be the matrix of singular vectors of the smallest singular values, and $hat{d}(e^{itheta},e^{-itheta})$ be the corresponding trigonometric polynomial for $hat{U}_{m+1:n}$ . Two possible generalization of the MUSIC method are as follows.

Spectral MUSIC: identify , for as the values on where $1/hat{d}(e^{itheta},e^{-itheta})$ achieves the -highest local maxima.
Root MUSIC: identify , for as the angle of the -roots of $hat{d}(z, z^{-1})$ which have amplitude and are closest to unit circle.

Signal subspace analysis

A signal subspace method relies on the fact that for the pair , and $U_{1:m}$ given above, there is a unique solution to the following matrix equation

$U_{1:m} =[B ~ AU_{1:m}] left[begin{array}{c}muPhi end{array} right],$

where is 1times m row vector and Phi a mtimes m matrix. The eigenvalues of Phi are precisely $e^{itheta_k}$ for k=1,ldots, m . If is a Toeplitz matrix and , given as in (1) with a companion matrix, then we recover the ESPRIT result [P. Stoica, R.L. Moses, 1997].

The signal subspace estimation is computed using sm.m, whereas music.m and esprit.m implement the MUSIC and ESPRIT methods, respectively. For the example discussed above, the estimated spectral lines are shown in the following figure. For an extensive comparison of the high resolution method sm.m with MUSIC and ESPRIT methods, see the example.

smspectrum

% subspace method
[thetas,residues]=sm(R,A,B,2);
arrowb(thetas,residues); hold on;
Ac=compan(eye(1,6));
Bc=eye(5,1);
That=dlsim_complex(Ac,Bc,y'); 
[th_esprit,r_esprit]=sm(That,Ac,Bc,2); % ESPRIT spectral lines
arrowg(th_esprit,r_esprit);

Spectral envelop

Let sigma(theta) denote a spectral measure of the stochastic process y_t , the envelop of maximal spectral power is defined as

$rho(theta) := sup_{sigma} left{lim_{epsilon rightarrow 0} (sigma(theta+epsilon)-sigma(theta-epsilon)) ~mid~ dsigma geq0, R=frac{1}{2pi} int_{-pi}^{pi} G(e^{itheta})dsigma(theta) G(e^{itheta})^* right},$

where represents the state covariance and G(z) is the input-to-state filter. In other words, rho(theta) represents the maximal spectral “mass” located at theta which is consistent with the covariance matrix. It can also be shown that

$rho(theta) =maxleft{ rho ~mid~ R-rho G(e^{itheta})G(e^{itheta})^* geq 0right},$

where $G(e^{itheta})G(e^{itheta})^*$ represents the state covariance for the sinusoidal input $left{e^{iktheta}, kin mathbb{Z}right}$ . The optimal solution is

$rho(theta)=frac{1}{G(e^{itheta})^* R^{-1}G(e^{itheta})},$

and implemented in envlp.m. This generalizes the Capon spectral estimation method which applies to being Toeplitz matrix and G(z) the tapped delay line transfer function. The Capon method may be motivated by a noise suppression formulation aimed at antenna array applications [P. Stoica, R.L. Moses, 1997].

For real valued processes with a symmetric spectrum with respect to the origin, a symmetrized version of the spectral envelop [1] can be similarly obtained:

$rho_{mbox{real}}(theta)=maxleft{rho ~mid~ R-frac{rho}{2} left(G(e^{itheta})G(e^{itheta})^*+G(e^{-itheta})G(e^{-itheta})^* right)right}.$

specenvlp

% spectral envelop
rhohalf = envlp(R,A,B,th);
rho = rhohalf.^2;
plot(th,rho);

Reference

[1] T.T. Georgiou, “Spectral Estimation via Selective Harmonic Amplification,” IEEE Trans. on Automatic Control, 46(1): 29-42, January 2001. [BIBTEX]

[2] T.T. Georgiou, “Spectral analysis based on the state covariance: the maximum entropy spectrum and linear fractional parameterization,” IEEE Trans. on Automatic Control, 47(11): 1811-1823, November 2002. [BIBTEX]