SpaceTimeCovMatMSE

SpaceTimeCovMatMSE(R,N,T,mode);

  SpaceTimeCovMatMSE(R,N,T) determines the mean square estimation error of 
  a space time covariance matrix R over a support length of (-T ... +T) 
  given a data set containing N snap shots, based on the analysis in [1].

  Via [xi_1,xi_2] = SpaceTimeCovMatMSE(R,N,T), this function returns two
  arguments: the estimation error due to the finite sample set in xi_1, 
  and a truncation error in xi_2. If 2T+1 exceeds the support of R, then 
  xi_2 will be zero.

  [xi_1,xi_2,VarE] = SpaceTimeCovMatMSE(R,N,T) additionally returns the
  variance of the estimation error for each element of the estimate.
     
  [xi_1,xi_2,VarE] = SpaceTimeCovMatMSE(R,N,T,'vec') returns vectors of 
  dimension T for xi_1 and xi_2, containing errors for lags tau, 0<tau<=T.

  Input parameters:
    R         MxMxL space-time covariance matrix
    N         number of snapshots
    T         maximum lag range over which estimate is calculated
    mode      optional, 'vec' for vectorial rather than total error terms  
     
  Output parameters
    xi_1      estimation error due to finite sample size N
    xi_2      truncation error, xi_2>0 in case T<L
    VarE      MxMx(2T+1) variance of estimated space time covariance
                matrix, element by element

  Reference:
    [1] C. Delaosa, J. Pestana, N. Goddard, S. Somasundaram, S. Weiss:
        "Sample Space-Time Covariance Matrix Estimation", submitted to 
        IEEE ICASSP, Brighton, UK, May 2019.

0001 function [xi_1,xi_2,VarE] = SpaceTimeCovMatMSE(R,N,T,VecMode);
0002 %SpaceTimeCovMatMSE(R,N,T,mode);
0003 %
0004 %  SpaceTimeCovMatMSE(R,N,T) determines the mean square estimation error of
0005 %  a space time covariance matrix R over a support length of (-T ... +T)
0006 %  given a data set containing N snap shots, based on the analysis in [1].
0007 %
0008 %  Via [xi_1,xi_2] = SpaceTimeCovMatMSE(R,N,T), this function returns two
0009 %  arguments: the estimation error due to the finite sample set in xi_1,
0010 %  and a truncation error in xi_2. If 2T+1 exceeds the support of R, then
0011 %  xi_2 will be zero.
0012 %
0013 %  [xi_1,xi_2,VarE] = SpaceTimeCovMatMSE(R,N,T) additionally returns the
0014 %  variance of the estimation error for each element of the estimate.
0015 %
0016 %  [xi_1,xi_2,VarE] = SpaceTimeCovMatMSE(R,N,T,'vec') returns vectors of
0017 %  dimension T for xi_1 and xi_2, containing errors for lags tau, 0<tau<=T.
0018 %
0019 %  Input parameters:
0020 %    R         MxMxL space-time covariance matrix
0021 %    N         number of snapshots
0022 %    T         maximum lag range over which estimate is calculated
0023 %    mode      optional, 'vec' for vectorial rather than total error terms
0024 %
0025 %  Output parameters
0026 %    xi_1      estimation error due to finite sample size N
0027 %    xi_2      truncation error, xi_2>0 in case T<L
0028 %    VarE      MxMx(2T+1) variance of estimated space time covariance
0029 %                matrix, element by element
0030 %
0031 %  Reference:
0032 %    [1] C. Delaosa, J. Pestana, N. Goddard, S. Somasundaram, S. Weiss:
0033 %        "Sample Space-Time Covariance Matrix Estimation", submitted to
0034 %        IEEE ICASSP, Brighton, UK, May 2019.
0035   
0036 % S. Weiss, UoS, 28/10/18
0037 %    return a vector of errors       06/05/2022
0038 
0039 %------------------------------------------------------------------------------
0040 %  parameters
0041 %------------------------------------------------------------------------------
0042 [M,~,L] = size(R);
0043 Tmax = (L-1)/2;            % lag range of R is (-Tmax ... +Tmax)
0044 ComplexMode = isreal(R);   % 0: complex;  1: real
0045 VectorMode=0;
0046 if nargin>3,
0047    if VecMode=='vec',
0048       VectorMode=1;
0049    end;
0050 end;      
0051    
0052 %------------------------------------------------------------------------------
0053 %  check for truncation
0054 %------------------------------------------------------------------------------
0055 if VectorMode==0,
0056   xi_2 = 0.0;
0057   if T<Tmax,                 % only evaluate positive lag range
0058     %   disp('truncation error incurred');
0059     for t = T+1:Tmax,   
0060       xi_2 = xi_2 + 2*norm(R(:,:,Tmax+1+t),'fro')^2;
0061     end;
0062   end;
0063 else
0064   dummy = zeros(Tmax,1);
0065   for t=1:Tmax,
0066      dummy(t)= 2*norm(R(:,:,Tmax+1+t),'fro')^2;
0067   end;
0068   dummy = cumsum(dummy,'reverse'); 
0069   xi_2 = zeros(T,1);
0070   if T<=Tmax,
0071      xi_2 = dummy(1:T);
0072   else; 
0073      xi_2(1:Tmax) = dummy;
0074   end;        
0075 end;
0076 
0077 %------------------------------------------------------------------------------
0078 %  determine estimation error
0079 %------------------------------------------------------------------------------
0080 if T > Tmax,
0081    L = T-Tmax;
0082    dummy = zeros(M,M,2*T+1);
0083    dummy(:,:,L+1:L+2*Tmax+1)=R;
0084    R = dummy; Tmax = T; 
0085 end;
0086 VarE = zeros(M,M,2*T+1);
0087 for m = 1:M,
0088   r_xx = squeeze(R(m,m,:)).'; 
0089   for mu = 1:M,
0090     r_yy = squeeze(R(mu,mu,:)).';
0091     r_xx_yy_ext = [zeros(1,T) r_xx.*conj(r_yy) zeros(1,T)];   
0092     Var_rh_xy2 = zeros(1,2*T+1);   
0093     rb_xy = squeeze(R(m,mu,:)).';
0094     if ComplexMode==0,
0095       rb_xy = zeros(size(rb_xy));
0096     end;
0097     rb_xy_ext = [zeros(1,T) rb_xy zeros(1,T)];      
0098     for tau = -T:T,
0099       W = [1:(T+1), T:-1:1] + (N-abs(tau)-T-1);    % weighting function
0100       W1 = r_xx_yy_ext((-T:T)+Tmax+T+1);
0101       W2 = rb_xy_ext((-T:T)+Tmax+T+1+tau);
0102       Var_rh_xy2(tau+T+1) = (W1 + W2.*fliplr(conj(W2)) )*W'/(N-abs(tau)).^2;
0103     end;
0104     VarE(m,mu,:) =  Var_rh_xy2;
0105   end;
0106 end;
0107 if VectorMode==0,
0108   xi_1 = sum(sum(sum(VarE)));
0109 else
0110   xi_1 = zeros(T,1);
0111   xi_1(1) = sum(sum(VarE(:,:,T+1))); 
0112   xi_1(2:T) = 2*cumsum(squeeze(sum(sum(VarE(:,:,T+2:2*T+1)))));  
0113 end;