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 R 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

 R-ARA^*=BH+H^*B^*,

for some row-vector H. 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 A 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 R 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 A. In the above example, we set n=5 eigenvalues of A 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

 AA^*+BB^*=I.

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 3dB 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,

 int -p(x) log(p(x)) dx,

quantifies the uncertainty of the corresponding random variable. When p(x) is a multi-dimensional Gaussian distribution with zero mean and covariance matrix M, 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 N 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 Szegddot{o}-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 R 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 Carathacute{e}odory-Fejacute{e}r-Pisarenko result for Toeplitz matrices. Given a ntimes n positive definite Toeplitz matrix T, 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 Q be the unique solution to the Lyapunov equation

 Q-AQA^*=BB^*.

The matrix Q is the state covariance when the input is pure white noise. Let now R be an arbitrary state covariance, let c be the smallest eigenvalues of the matrix pencil R-rho Q, and assume that the rank of R-c Q is m, 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 U 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 m and the last n-m columns of U, 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 m 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.

  1. Spectral MUSIC:  identify theta_k, for k=1,ldots, m as the values on [-pi,; pi] where 1/hat{d}(e^{itheta},e^{-itheta}) achieves the m-highest local maxima.

  2. Root MUSIC:  identify theta_k, for k=1,ldots, m as the angle of the m-roots of hat{d}(z,  z^{-1}) which have amplitude <1 and are closest to unit circle.

Signal subspace analysis

A signal subspace method relies on the fact that for the pair A, B 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 mu 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 R is a Toeplitz matrix and A, B given as in (1) with A 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 R 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 R 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]