% ==================================================
%
%   This file statistically computes the CDL-A channel covariance. 
%
%   The covariance will be used for MMSE estimation of CDL-A channel. 
%
%
%
%
% ==================================================
clear; close all; clc; 

addpath('algs/'); 
addpath('utils/');
fc          = 3.5e9; 
c           = physconst('lightspeed'); 
lambda      = c/fc; 
k0          = 2*pi/lambda; 

rng(0); 

% The Array Aperture should be set to match that of the sim .m files. 
L           = 1;
% L           = 3; 
N_test      = 1e4; 
P           = 100;
sigma2      = 1; 
Ncombiner   = 120; 

% Setup bandwidths. 
reg.Gamma       = 0.05; 
RegionsTx       = cell(1,1);
RegionsTx{1}    = [0, 0, 0.15, reg.Gamma];   
RegionsRx       = cell(1,1); 
RegionsRx{1}    = [0, 0, 0.15, reg.Gamma];   
reg.RegionsTx   = RegionsTx; 
reg.RegionsRx   = RegionsRx; 

% MIMO = setup("SIMOULAwithFixedK", fc, L, reg);     
MIMO = setup("SymmetricULAwithFixedK", fc, L, reg);  

NRx = MIMO.RxArray.N; NTx = MIMO.TxArray.N;
Nh = NRx*NTx; 


cdl_model = nrCDLChannel('ChannelFiltering', false);
cdl_model.DelayProfile              = 'CDL-A'; 
cdl_model.MaximumDopplerShift       = 300;
cdl_model.CarrierFrequency          = MIMO.fc;
cdl_model.SampleRate                = 10e2;     
cdl_model.SampleDensity             = Inf;
cdl_model.NumTimeSamples            = 1;
cdl_model.TransmitAntennaArray.Size = [MIMO.TxArray.N, 1, 1, 1, 1];   % [M N P Mg Ng]
cdl_model.ReceiveAntennaArray.Size  = [MIMO.RxArray.N, 1, 1, 1, 1]; 
cdl_model.Seed                      = randi(10000);

corrMat = zeros(Nh); 

% figure(1);
FR = dftmtx(NRx); FT = dftmtx(NTx); 
is_disp = false; 

for idx = 1:N_test
    [pathGains, ~] = cdl_model();
    
    H = reshape(sum(pathGains(1, :, :, :), 2), [NTx, NRx]).';
    
    % H = H/norm(H, 'fro'); 
    h = H(:); 
    corrMat = corrMat + (h*h'); 
    
    % Get statistical correlation of H. 
    % Correlation.CovH = eye(NRx*NTx)/(NRx*NTx); 
    fprintf('%d\n', idx); 

    if is_disp
        Htilde = (FR')*H*(FT); 
    
        imagesc(mag2db(abs(Htilde))); 
        colormap('turbo'); 
        xlabel('Tx Beam Index'); ylabel('Rx Beam Index'); colorbar; 
        clim([-30, 30]); 
    
        
        pause(0.01); 
    end

end

corrMat = corrMat/(N_test-1); 

%%
save('results/CDL_corrMat.mat', 'corrMat', 'MIMO'); 
fprintf('MAT file saved.\n'); 

