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 thoughts 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 stucture od Toeplitz covariances.

In contrast, recently, the analysis of state covariance matrices, see e.g. [1, 2], suggested a general framework which alows 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));