%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Experiments with OFDMA Resource Allocation %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
warning('off','MATLAB:divideByZero')
CSI_Type = 0; % 0 - Perfect, 1 - Predicted
PLOT_CHAN = 0;
% close all
Mvec = (4:4:20);                                  % No. of users
Nfftvec = 128;
Nu = 76;
CPfrac = 6/128;                          % Fraction of Nfft for CP
colorvec = ['m','b','r','g'];

Nframes = 10;
Nsymperframe = 100;
BW = 1.25e6;                               % Bandwidth (Hz)
fc = 2.6e9;                             % Carrier frequency (Hz)
[tau,gnorm] = Pdp('vehicularA');        % Power Delay Profile
c = 3e8;                                % Speed of light (m/s)
Velvec = [30, 75, 120, 180, 250]*1000/3600;
v = Velvec(5);                          % Velocity (m/s)
fd = v/c*fc;                            % Doppler Frequency (Hz)
tcoh = 1/(2*fd);                        % Coherence time (s)
lambda = 3e8/fc;                        % Wavelength (m)
fs = 1.92e6;           % Sampling Frequency (Hz)


% Prediction params
Lpred = 4;
Dt = 1;

SNRvec = 10;

BER = 1e-3;
R = 0:2:6;
Eta = -log(5*BER)/1.6*(2.^R-1);
Nfft = Nfftvec(1);
tsym = Nfft*(1+CPfrac)/fs;

C0 = zeros(length(Mvec),length(SNRvec));
P0 = zeros(length(Mvec),length(SNRvec));
G0 = zeros(length(Mvec),length(SNRvec));
I0 = G0;
Iint0 = G0;

C = zeros(length(Mvec),length(SNRvec));
G = zeros(length(Mvec),length(SNRvec));
I = G;

C2 = zeros(length(Mvec),length(SNRvec));
G2 = zeros(length(Mvec),length(SNRvec));
I2 = G;

C3 = zeros(length(Mvec),length(SNRvec));

for m = 1 : length(Mvec)
  M = Mvec(m);
  gbarvec = ones(1,M);
  for nf = 1 : Nframes
    w = unifrnd(.1,.9,1,M);
    w = w/sum(w);
    for s = 1 : length(SNRvec)
      SNR = SNRvec(s);    % dB
      sigma = 10^(-SNR/20); % standard deviation
      for ns = 1 : Nsymperframe
        [M,nf,s,ns]
        H = channel(M, 1, fs, Nu, Nfft, tau, gnorm);
        Gamma = abs(H).^2/(sigma^2);        
        if ns == 1
          [Gtmp0, Itmp0, Itmpint0, Ltmp0, Ctmp0, Ptmp0]=ExpectedDiscretePerfectCSI2(gbarvec/sigma^2,Nu,w,Nu,R,Eta);
          if Gtmp0>1
            keyboard
          end
          C0(m,s) = C0(m,s) + sum(Ctmp0)/Nu/Nframes;
          G0(m,s) = G0(m,s) + Gtmp0/Nframes;
          I0(m,s) = I0(m,s) + Itmp0/Nframes;
          Iint0(m,s) = Iint0(m,s) + Itmpint0/Nframes;
        end
        [Gtmp, Itmp, Ltmp, Ctmp, Ptmp] = ExpectedDiscretePerfectCSI(Gamma,Gamma,w,Nu,R,Eta,Ltmp0);
        [Gtmp2, Itmp2, Ltmp2, Ctmp2, Ptmp2] = ExpectedDiscretePerfectCSI(Gamma,Gamma,w,Nu,R,Eta,0);
        Ctmp3 = WeightedDiscretePerfectCSIConst(Gamma,w,Eta,R);

        C(m,s) = C(m,s) + sum(Ctmp)/Nu/Nframes/Nsymperframe;
        G(m,s) = G(m,s) + Gtmp/Nframes/Nsymperframe;
        I(m,s) = I(m,s) + Itmp/Nframes/Nsymperframe;

        C2(m,s) = C2(m,s) + sum(Ctmp2)/Nu/Nframes/Nsymperframe;
        G2(m,s) = G2(m,s) + Gtmp2/Nframes/Nsymperframe;
        I2(m,s) = I2(m,s) + Itmp2/Nframes/Nsymperframe;

        C3(m,s) = C3(m,s) + sum(Ctmp3)/Nu/Nframes/Nsymperframe;
      end; % SNR
    end; % Nsym
  end; % Nframes
end; % No. of users
figure
bar([C0,C,C2,C3])
Diff = [(C-C3)./C3, (C-C2)./C2];
figure
bar(Diff)
str = ['save DiscRatePerfCSIMuser_SNR_',num2str(SNRvec),'_Nf_',num2str(Nframes),'_Nsym_',num2str(Nsymperframe)];
eval(str)