Separation of line spectra and system dynamics

Consider the following time-series

 y_t =sum_{k=1}^m a_ksin(omega_k t+phi_k)+{rm noise}_t, {rm ~for~} tin mathbb{Z},

where {{rm noise}_t} represents colored noise. The routine L1AR.m estimates spectral lines together with an autoregressive (AR) component corresponding to the noise. The routine and relevant theory are explained in [1] and rely on the ideas of compressive sensing (L1-optimization, sparse selection from a suitable dictionary).

help L1AR
  function [omega,mag,phi,a,sigma]=L1AR(y,freqNum,ARorder,freqrange,N)
               y: observation record signal+noise
         freqNum: number of signals
         ARorder: order of autoregressive filter (default value 0)
       freqrange: interested frequency interval (default value [0 pi])
               N: resolution=2*pi/N (default value 2*length(y))
          omega:  estimated frequencies of signals
            mag:  magnitudes of signals
            phi:  estimated phase angles
              a:  coefficients of AR filter
          sigma:  whitened noise variance

Needed software: convex optimization tool cvx.

Sample problem

Consider the following time-series

 y_t = a_1sin(omega_1 t+phi_1)+ a_2sin(omega_2 t+phi_2)+a_3sin(omega_3 t+phi_3)+{rm noise}_t, {rm ~for~} t=0, ldots, N-1,

where N=150, a_1=1, a_2=a_3=0.5, omega_1=frac{59pi}{200}, omega_2=frac{139pi}{200}, omega_3=frac{141pi}{200}, phi_1=frac{pi}{4}, phi_2=frac{pi}{3} and phi_3=-frac{pi}{3}. The resolution needed is frac{pi}{100} which is higher than Fourier analysis that is frac{pi}{75}. The noise component is generated as the output of an AR filter having transfer function


and driven by white noise with variance 2.25. The periodogram of the observation record is shown in the following figure.


There are multiple peaks in the plot. The “ground truth” in the form of noise spectral density and spectral lines are shown in the first subplot of the following figure. The noise spectral density and spectral lines which are approximated using L1AR.m are shown in the second subplot. All three frequencies are identified exactly at the chosen resolution scale. The code for this example is shown below.

% parameters for signals
o1   = 59*pi/200;  phi1 = pi/4;  a1   = 1;
o2   = 139*pi/200; phi2 = pi/3;  a2   = .5;
o3   = 141*pi/200; phi3 = -pi/3; a3   = .5;
% parameters for noise
order = 1; sigma = 1.5;
a     = [1.0000   -.85];
% generate data
n     = 150; t=0:n-1;
T     = tf([1 zeros(1,order)],a,-1);
randn('seed',101); noise = sigma*randn(n,1);
sig   = sin(o1*t+phi1)+.5*sin(o2*t+phi2)+.5*sin(o3*t+phi3);
y     = lsim(T,noise)+sig';
% estimate
freqNum = 3; N = 200;
[omega,mag,phi,a,sigma_est]=L1AR(y,freqNum,order,[0 pi],N);
% plot results
w     = linspace(0,pi,500);
ptrue = abs(freqresp(T,w)).^2;
T_est = tf([1 zeros(1,order)],a,-1);
pest  = abs(freqresp(T_est,w)).^2*sigma_est^2;
figure(1); periodogram(y);
        plot(w,vec(ptrue),'r','LineWidth',1.2);hold on;
        legend('true spectrum of colored noise','true spectral lines');
        set(gca,'YScale','log','XLim',[0 pi]);
        arrow([o1 o2 o3],(2*pi*[a1 a2 a3]),12);
        plot(w,vec(pest),'b','LineWidth',1.2);hold on;
        legend('estimaged spectrum of colored noise','estimated spectral lines');
        set(gca,'YScale','log','XLim',[0 pi]);


[1] L. Ning, T.T. Georgiou and A. Tannenbaum, “Separation of system dynamics and line spectra via sparse representation,” 49 th IEEE Conference on Decision and Control, 473-478, December 2010.