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

% function [Tsignal,Tnoise]=toprange4(T,m) : decomposes T into T_signal + T_{MA(m)-noise}
%
%  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.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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]);
    %tstructureones=toeplitz(ones(size(kk)));
    M=eye(n1,m1);
    f=eye(1,m1); g=eye(m1); g(:,1)=[]; % D=.5*f*Tnoise*f'; C=f*Tnoise*g;
else,
    zz=zeros(1,n-m); tstructure=toeplitz([kk zz]);
    %tstructureones=toeplitz([ones(size(kk)) zz]);
    M=eye(n1);
    f=eye(1,n1); g=eye(n1); g=g(:,2:m+1); % D=.5*f*Tnoise*f'; C=f*Tnoise*g;
end
      
setlmis([])
Tnoise=lmivar(3,tstructure);
P=lmivar(1,[m,1]);
A=compan(eye(1,m1)); B=eye(m,1);

% Tnoise<T
lmiterm([1 1 1 Tnoise],M,M')
lmiterm([-1 1 1 0],T);

% 0<[P-A'PA C'-A'PB; C-B'PA D+D'-B'PB]>0
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);

% 0<P
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 (revised December 2005)
