function QuantH = quantiz_channels(H, regions, P, QuantizType, CQIType)
% QuantH = quantiz_channels(H, regions, P, QuantizType, CQIType)
%
% Quantizes channels H using a codebook of size "regions", using:
%   QuantizType quantization (0 - No Quantization, 1- Random Vector Quantization, 2 - Quantization Upper Bound, 3 - Scalar Quantization), and
%   scaled by the CQI value depending on CQIType (0 - Channel Norm, 1 - SINR computed using quantization error and SNR P)
% H and QuantH are of dimensions (Number of users) X (Number of base station antennas)
%
%   Part of mud_vs_csi package.
%
%   Reference:  "Multi-User Diversity vs. Accurate Channel State Information in MIMO Downlink Channels"
%                       by Niranjay Ravindran and Nihar Jindal, Submitted to IEEE Transactions on Communications, 2009.
%

K = size(H, 1);
Nt = size(H, 2);
QuantH = zeros(K, Nt);

if QuantizType ~= 3
    %Choose quantized channels by generating random quantization errors
    temp = rand(K,1);

    if QuantizType == 0
        cos_sq = ones(1, K); %Perfect CSI
    elseif QuantizType == 1 %RVQ
        if log2(regions) <= 80
            cos_sq = 1 - (1-temp.^(1/regions)).^(1/(Nt-1));
        else %Very large codebook, use symbolic toolbox for higher precision (slow). Comment out if symbolic toolbox is not installed.
            cos_sq = double(vpa(1 - (1-temp.^(1/sym(regions))).^(1/(Nt-1))));
        end
    else %QUB
        if log2(regions) <= 80
            cos_sq = 1 - (temp./regions).^(1/(Nt-1));
        else  %Very large codebook, use symbolic toolbox for higher precision (slow). Comment out if symbolic toolbox is not installed.
            cos_sq = double(vpa(1 - (temp./sym(regions)).^(1/(Nt-1))));
        end
    end
    for k = 1:K
        channel_gain = norm(H(k, :));
        rv = H(k, :)/channel_gain;
        temp = randn(Nt,1) + i * randn(Nt,1);
        temp = (temp - rv * temp * rv')';
        if CQIType == 0
            CQI = channel_gain^2; %Norm CQI
        else
            CQI = channel_gain^2*cos_sq(k) / (1 + P*channel_gain^2*(1-cos_sq(k))); %SINR CQI
        end
        QuantH(k,:) = sqrt(CQI)* (sqrt(cos_sq(k)) * rv + sqrt(1-cos_sq(k)) * temp/norm(temp));
    end
else
    %Scalar Quantization
    for k = 1:K
        %Split bits evenly between magnitude and phase
        rpc = floor(regions^(1/(2*(Nt-1)))); %Quantization regeons per component

        mag_bin_size = 1/2* pi / rpc;
        phase_bin_size = 1/2* 4*pi / rpc;

        temp = H(k, 2:Nt) / H(k, 1);
        mag_vector = abs(temp);
        phase_vector = angle(temp)+pi;

        mag_vector = (ceil(mag_vector/mag_bin_size) * mag_bin_size) - (mag_bin_size/2); %Quantize relative magnitudes
        phase_vector = (ceil(phase_vector/phase_bin_size) * phase_bin_size) - (phase_bin_size/2); %Quantize relative phases
        phase_vector  = phase_vector - pi;

        QuantH(k,:) = [1 mag_vector.*exp(i*phase_vector)];
        QuantH(k,:) = QuantH(k, :)/norm(QuantH(k,:));
        channel_gain = norm(H(k, :));
        cos_sq = abs(QuantH(k, :)*H(k, :)')^2/channel_gain^2;
        if CQIType == 0
            CQI = channel_gain^2; %Norm CQI
        else
            CQI = channel_gain^2*cos_sq / (1 + P*channel_gain^2*(1-cos_sq)); %SINR CQI
        end
        QuantH(k,:) = sqrt(CQI) * QuantH(k, :);
    end

end %End if doscalar

%__EOF__
