function [tsignal,tnoise]=MA_decomposition(T,m)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [Tsignal,Tnoise]=MA_decomposition(T,m)
%
%  T = Tsignal + Tnoise
%
% where Tnoise corresponding to a MA(m) process
% with spectrum Q0+2*Q1*cos(theta)+...+2*Qm*cos(m*theta)
% i.e., with m+1 coefficients.
%
% T,tsignal,tnoise are all Toeplitz covariance matrices.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

narg=nargin;
if narg<2, disp('needs at least 2 arguments'), return, end
if m==0, tnoise=eye(size(T))*min(eig(T)); tsignal=T-tnoise; return, end
n1=length(T); Neps=1e-15; m1=m+1; n=n1-1;
kk=1:m1; 
if m>n,
    tstructure=toeplitz([kk]);
    M=eye(n1,m1);
    f=eye(1,m1); g=eye(m1); g(:,1)=[];
else,
    zz=zeros(1,n-m); tstructure=toeplitz([kk zz]);
    M=eye(n1);
    f=eye(1,n1); g=eye(n1); g=g(:,2:m+1);
end
setlmis([])
Tnoise=lmivar(3,tstructure);
P=lmivar(1,[m,1]);
A=compan(eye(1,m1)); B=eye(m,1);
lmiterm([1 1 1 Tnoise],M,M')
lmiterm([-1 1 1 0],T);
lmiterm([-2 1 1 P],1,1);lmiterm([-2 1 1 P],-A',A);
lmiterm([-2 1 2 Tnoise],g',f');lmiterm([-2 1 2 P],-A',B);
%lmiterm([-2 2 1 Tnoise],f,g);lmiterm([-2 2 1 P],-B',A);
lmiterm([-2 2 2 Tnoise],f,f');lmiterm([-2 2 2 P],-B',B);
lmiterm([2 1 1 0],0);
lmiterm([-3 1 1 P],1,1);
lmiterm([3 1 1 0],0);
lmisys=getlmis;
if m>n,
    weights=eye(m+1,m+1);
else,
    weights=eye(n+1,n+1);
end
c=mat2dec(lmisys,weights,0);
options=[1e-13,0,0,0,1];
[copt,xopt]=mincx(lmisys,-c,options);
trcq=-copt;

tnoise=dec2mat(lmisys,xopt,Tnoise);
tsignal=T-M*tnoise*M';
% Written by T.T. Georgiou (August 2005)