x= DataV'; %Data collected by the ULA (31 Spatial samples-31Antennas), 501 time samples
N = size(x,1); % Number of antennas
M = size(x,2); % Number of samples in fast-time dimension
% window dimension for the smoothing dim[m1xm2]
m1 = 25;
m2 = 300;
p1 = N-m1+1; %positions in the spatial dimension
p2 = M-m2+1; %positions in the time dimension
N_submat = p1*p2; %number of all the possible scan
M_tot = m1*m2;
ind_tot=1;
for p1_ind = 1:p1
 for p2_ind = 1:p2
     x_win_new(:,ind_tot) = reshape (x(p1_ind:m1+(p1_ind-1),p2_ind:m2+(p2_ind-1)),[1,M_tot])';
     ind_tot=ind_tot+1;
 end
end
X_win_cor = x_win_new*x_win_new';
X_win_cor_conj = conj(X_win_cor);
J= fliplr(eye(M_tot));  %transition matrix
%data smoothed covariance matrix (forward+backward)
C=(1/2*N_submat)*(X_win_cor + J*X_win_cor_conj*J);

x x x .      . x x x      . . . .      . . . .
x x x .      . x x x      x x x .      . x x x
x x x .      . x x x      x x x .      . x x x
. . . .      . . . .      x x x .      . x x x