correlation matrix for complex white noise does not equal to identity matrix?
Why don't you simply compute $X^{H}  X / N$ ?
clear all;
close all;

fs = 4000;
N = 100*fs;
R = 0.5;% overlap percent
Nseg = fs;

M = 8;

I noticed your data was not complex, since you metioned complex
in the question here I made it complex.
X = 0.5*randn(N,M) + 0.5j *randn(N,M);

This will be close to the identity matrix there will be small complex perturbations
the expected perturbations are inversely proportional to $\sqrt{N}$. The $X^H$ corresponds to X', it is the conjugate transpose of $X$.
R_hat = X' * X / N

Edit
Here I compute RR such that RR(:,:,k) is the accumulated covariance matrix of the frequency bin $k$ for the different inputs.
clear all;
close all;

fs = 4000;
N = 100*fs;
R = 0.5;% overlap percent
Nseg = 256;

M = 8;

for m = 1:M
   x(:,m) = randn(N,1) ;    
end  

segNum = floor((N-R*Nseg)/((1-R)*Nseg));

RR = zeros(M, M, Nseg);

for n = 1:segNum
   y = x([1:Nseg]+(n-1)*Nseg*(1-R),:) ;
   Y = fft(y);
   for k = 1:Nseg
     %% the correlation for for one bin one slice.
     RR(:,:,k) = RR(:,:,k) + Y(k, :)' * Y(k, :);
   end
end
%% Normalize
RR = RR * (1/segNum/Nseg);

%% Acculate over all bins
R_hat = 0;
for n = 1:Nseg
   R_hat = R_hat + RR(:,:,k); 
end

R_hat = R_hat/Nseg