Movingaverage decomposition
In the CarathodoryFejrPisarenko subspace methods, a Toeplitz structured covariance matrix is decomposed as
where is a singular covariance matrix attributed to the spectral lines and a covariance attributed to white noise. The estimation of spectral lines is then based on .
A possible generalization of this idea is to allow the noise component to have a movingaverage (MA) spectrum [1] . Thus, we decompose as
where is a covariance matrix corresponding to a MA process. Let be a companion matrix of size and
Then
is the covariance matrix of a MA process if and only if there exist a such that
This is a consequence of the KYP positivereal lemma [1]. Assuming that is the order of the MA process, it is natural to subtract the maximum amount of noise from . Thus we are called to solve
This decomposition is obtained using routine MA_decomp.m which in turn uses the convex optimization package cvx.
An alternative approach is proposed in [2].
Let
be the vector of the coefficients of a th order MA filter, and let
be the corresponding covariance sequence of the output of the filter driven by white noise with unit variance.
Then
where denotes the sum of the th subdiagonal of a matrix.
This was relaxed to the constraint used in [2] as follows
Example
Consider the following sequence
with and , , and the noise component generated by a MA filter driven by
white noise (complex with independent real and imaginary parts). The MA filter is rd with zeros at , , and respectively.
MA
N=100; n=15;
a0=1; a1=1.5; theta1=1.3; a2=2; theta2=1.35;
k=1:N; k=k(:);
y=a1*exp(1i*(theta1*k+2*pi*rand))+a2*exp(1i*(theta2*k+2*pi*rand));
b=[1.0000 0.4318 0.4496 0.5440];
u_real=randn(1,N+length(b)1);
u_imag=randn(1,N+length(b)1);
u=a0*(u_real+i*u_imag)/sqrt(2);
noise=b*toeplitz(u(length(b):1:1),u(length(b):end));
y=y'+noise;
Ac=compan(eye(1,n+1));
Bc=eye(n,1);
That=dlsim_complex(Ac,Bc,y);
NN=2048; th=linspace(0,2*pi,NN);
T=ToepProj(That);
[Ts, Tma]=MA_decomp(T,3);
figure(1);clf
h=freqresp(f,th);
subplot(3,1,1);
plot(theta1,a1,'r^','MarkerSize',12);hold on;
plot(th,vec(abs(h)).^2,'b','LineWidth',1.2);
legend('true spectral line','true MA spectrum','Location','North');
arrow([theta1 theta2],[a1 a2],12);
set(gca,'xlim',[0 2*pi]);
set(gca,'YTick',[]);
r_ma=[fliplr(Tma(1,:)) Tma(1,2:end)];
p_ma=abs(fft(r_ma,length(th)));
[w_s,r_s]=sm(Ts,Ac,Bc,rank(Ts,0.001));
r_s=r_s(w_s<pi); w_s=w_s(w_s<pi);
subplot(3,1,2);
plot(w_s(1),r_s(1),'r^','MarkerSize',12);hold on;
plot(th,p_ma,'b','LineWidth',1.2);
legend('estimated spectral line','estimated MA spectrum','Location','North');
arrow(w_s,r_s,12);
set(gca,'xlim',[0 2*pi]);
set(gca,'YTick',[]);
[Ts0, Tn]=MA_decomp(T,0);
[w_s0,r_s0]=sm(Ts0,Ac,Bc,rank(Ts0,0.001));
r_s0=r_s0(w_s0<pi); w_s0=w_s0(w_s0<pi);
subplot(3,1,3);
plot(w_s0(1),r_s0(1),'r^','MarkerSize',12);hold on;
plot(th,Tn(1,1)*ones(size(th)),'b','LineWidth',1.2);
legend('estimated spectral line','estimated noise variance','Location','North');
arrow(w_s0,r_s0,12);
set(gca,'xlim',[0 2*pi]);
set(gca,'YTick',[]);
The first figure below shows the true spectral lines marked with arrows superimposed on the MA spectrum. The figure in the middle shows the estimated spectral lines and the estimated rd order MA spectrum. There are two dominant spectral lines locate very closely
to the true values. The figure at the bottom shows the estimate when the noise component is taken to be white.
The blue line shows the noise variance. In this case, only one dominant spectral line is detected near the true ones.
Reference
[1] T.T. Georgiou, “Decomposition of Toeplitz matrices via convex optimization,” IEEE Signal Processing Letters, 13(9): 537540, September 2006.
[2] P. Stoica, L. Du, J. Li and T. Georgiou, “A new method for movingaverage parameter estimation”, Signals, Systems and Computers (ASILOMAR), 2010 Conference Record of the Forty Fourth Asilomar Conference on, 1817–1820, 2010.
