% Compare exact versus approximate Theta
% M = 4;                                  % No. of users
K = 64;                                % No. of subcarriers
BW = 5e6;                               % 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)

SNRvec = [5:5:20];



% Generate channel
% wvec = 0;
C = zeros(M,length(wvec));
G = zeros(1,length(wvec));
I = G;
L = G;



Pbar = 1;

NumPoints = 50;
Thetavec = zeros(NumPoints,1);
Thetaestvec = Thetavec;
Mvec = [4,16,32];
for m = 1 : length(Mvec),
  figure
  M = Mvec(m);
  nidx = [-K/2:K/2-1]';
  E = exp(-sqrt(-1)*2*pi*fs/K*nidx*Tau);
  h = (gnorm'*ones(1,M)).*(randn(length(gnorm),M)+j*randn(length(gnorm),M))/sqrt(2);
  H = E*h;
  w = rand(M,1);
  w = w/sum(w);
  W = (ones(K,1)*w');
  for s = 1 : length(SNRvec),
    s
    sigma = 10^(-SNRvec(s)/20); % standard deviation
    G = abs(H).^2/sigma^2;
    lambdamin = K/log(2)*min(w)/(Pbar+sum(max(1./G.')));
    lambdamax = max(max(W.*G))/log(2);
    lambdavec = linspace(lambdamin,lambdamax,NumPoints);
    for i = 1:length(lambdavec),
      [Theta, maxP] = Dual(lambdavec(i),W,G,M,K,Pbar);
      Thetaest = DualEst(lambdavec(i),W,G,M,K,Pbar,K*log(M),K*log(2));
      Thetavec(i) = Theta;
      Thetaestvec(i) = Thetaest;
    end;
    [x,fval,exitflag,output] = fminbnd(@(x) Dual(x,W,G,M,K,Pbar),lambdamin,lambdamax);
    plot(lambdavec,Thetavec,'b-',lambdavec,Thetaestvec,'rx',x,fval,'ko')
    text(x,1.5*fval,)    
    hold on
  end;
end;


