%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Experiments with OFDMA Resource Allocation %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
warning('off','MATLAB:divideByZero')
CSI_Type = 0; % 0 - Perfect, 1 - Predicted
PLOT_CHAN = 0;
% close all
M = 2;                                  % No. of users
% Nfft = 2;                             % No. of subcarriers
% Nfftvec = [4,16,64,256];                  % No. of subcarriers
Nfftvec = 64;
Nu = 33;
CPfrac = 1/16;                          % Fraction of Nfft for CP
colorvec = ['m','b','r','g'];

% Nchan = 1000;
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(2);                          % Velocity (m/s)
fd = v/c*fc;                            % Doppler Frequency (Hz)
tcoh = 1/(2*fd);                        % Coherence time (s)
lambda = 3e8/fc;                        % Wavelength (m)
fs = floor(8/7*BW/8000)*8000;           % Sampling Frequency (Hz)


% Prediction params
P = 8;
Dt = 21;

% SNR = 0;    % dB
sigma = 10^(-SNR/20); % standard deviation

Eta = 
% Generate channel
wvec = [0.1:0.2:0.9];
% wvec = 0;
for n = 1 :length(Nfftvec),
  Nfft = Nfftvec(n);
  tsym = Nfft*(1+CPfrac)/fs;
  C = zeros(M,length(wvec));
  G = zeros(1,length(wvec));
  I = G;
  L = G;

  C2 = zeros(M,length(wvec));
  G2 = zeros(1,length(wvec));
  I2 = G;
  L2 = G;

  C3 = zeros(M,length(wvec));
  G3 = zeros(1,length(wvec));
  I3 = G;
  L3 = G;
  %   if CSI_Type == 0
  H = channel(M, Nchan, fs, Nu, tau, gnorm);
  %   else
  [Hhat,sigma2ehat] = channelpred(M, Nchan, fs, tsym, Nu, tau, gnorm, fd, P, Dt, CPfrac, sigma^2);
  E = (randn(size(H))/sqrt(2)+sqrt(-1)*randn(size(H))/sqrt(2)).*repmat(sigma2ehat,[1,M,Nchan]);
  %   end
  for c = 1 : Nchan
    c
    Gamma = abs(H(:,:,c)).^2./(sigma^2);
    Gammahat = abs(H(:,:,c)-E(:,:,c)).^2./(sigma^2);

    for i = 1:length(wvec),
      %         i = 6;
      w = [wvec(i),1-wvec(i)];
      [Gtmp, Itmp, Ltmp, Ctmp, P]=ExpectedWeightedPerfectCSI(Gamma,Gamma,w,1);
      [Gtmp2, Itmp2, Ltmp2, Ctmp2, P2]=ExpectedWeightedImperfectCSI(Gammahat,Gamma,sigma^2,sigma2ehat,w,1);
      [Gtmp3, Itmp3, Ltmp3, Ctmp3, P3]=ExpectedWeightedPerfectCSI(Gammahat,Gamma,w,1);
      f = sum(w.*Ctmp);
      C(:,i) = C(:,i) + Ctmp'/Nchan;
      G(i) = G(i) + Gtmp/Nchan;
      I(i) = I(i) + Itmp/Nchan;
      L(i) = L(i) + Ltmp/Nchan;

      f2 = sum(w.*Ctmp2);
      C2(:,i) = C2(:,i) + Ctmp2'/Nchan;
      G2(i) = G2(i) + Gtmp2/Nchan;
      I2(i) = I2(i) + Itmp2/Nchan;
      L2(i) = L2(i) + Ltmp2/Nchan;

      f3 = sum(w.*Ctmp3);
      C3(:,i) = C3(:,i) + Ctmp3'/Nchan;
      G3(i) = G3(i) + Gtmp3/Nchan;
      I3(i) = I3(i) + Itmp3/Nchan;
      L3(i) = L3(i) + Ltmp3/Nchan;
    end;
  end;
end;
figure(1)
plot(C(1,:),C(2,:),'-',C2(1,:),C2(2,:),'--',C3(1,:),C3(2,:),'-.')
grid on
hold on
title('Capacity Region for 2 Users')
xlabel('Capacity for user 1 (bps/Hz)')
ylabel('Capacity for user 2 (bps/Hz)')
legend('P-CSI: WF','I-CSI: Opt','I-CSI: WF')
figure(2)
semilogy(wvec,G,[colorvec(n),'x-'])
grid on
hold on
figure(3)
plot(wvec,I,[colorvec(n),'x-'])
hold on
grid on
figure(4)
plot(wvec,L,[colorvec(n),'x-'])
hold on
grid on
% end;
figure
plot(wvec,G)
grid on
title('Duality Gap')
xlabel('User 1 Weight')
figure
plot(wvec,10*log10(I))
title('No. of Iterations')

