% Function: mimo_gen_table
% Author: Stephen Abraham
%
% This function will generate the capacity PDF and CDF of an MxN MIMO link.
% It will then use the inputs provided by the user to generate the SNR
% table of goodput optimal rates for payload ranges.
%
% The link is modeled as follows:
% - M transmit antennas, N receive antennas
% - Rayleigh fading
% - H: channel gain matrix with CN(0,1) iid entries
%
% Parameters:
% - num_transmit - number of transmit antennas
% - num_receive - number of receive antennas
% - num_streams - number of streams to use
% - snr_string - list of SNRs (in dB) to generate PDFs and CDFs for
% - rates_string - list of rates (in Mbps) in protocol
% - sifs - SIFS length in microseconds
% - slot_time - slot time in microseconds
% - basic_rate - transmission rate (in Mbps) for control bits (preamble and
%                ACK)
% - cw - contention window size (minimum)
% - plcp_bits - number of PLCP bits
% - plcp_rate - PLCP transmission rate (in Mbps)
% - preamble_bits - number of bits in preamble
%
% Note: Parameters are entered as strings (surrounded by ''). See sample
%       mimoconfig.txt for how parameters are formated.
%
%
%
% 
% -------------------------------------------------
% Running mimo_gen_table.
% -------------------------------------------------
% 
% If you wish to generate the goodput optimal rate table and not use it in 
% NS-2, open MATLAB and call the mimo_gen_table function as you would a 
% normal MATLAB function. Note that the parameters need to be in string 
% format. Example:
% 
% mimo_gen_table('4','4','4','10,20,30','400,450,500,550,600,650,700,750','10','9','24','16','384', '24', '0')
% 
% 
% 
% -------------------------------------------------
% File format
% -------------------------------------------------
% 
% snrvalues.txt:
% 
% ...
% node1 node2 SNR(indB) distance(in meters)
% ...
% 
% 
% RateProb.txt:
% 
% rate1 rate2 rate3 ... 
% snr1 rate1_probability_of_success rate2_probability_of_success rate3_probability_of_success ...
% snr2 rate1_probability_of_success rate2_probability_of_success rate3_probability_of_success ...
% ...
% 
%
% SNRTable.txt:
% 
% snr1 snr2 snr3 ...
% First row is goodput optimal rates for payload ranges for snr1
% Second row is goodput optimal rates for payload ranges for snr2
% Third row is goodput optimal rates for payload ranges for snr3
% ...
% 
% 
% The ranges are in this format: goodput_optimal_rate,lower_bound_of_payload_size-upper_bound_of_payload_size
% 
% 
% 
% 

function mimo_gen_table(num_transmit, num_receive, num_streams, snr_string, rates_string, sifs, slot_time, basic_rate, cw, plcp_bits, plcp_rate, preamble_bits)

% Convert parameters to their integer values
num_transmit = str2num(num_transmit);
num_receive = str2num(num_receive);
num_streams = str2num(num_streams);
sifs = str2num(sifs);
slot_time = str2num(slot_time);
basic_rate = str2num(basic_rate);
cw = str2num(cw);
plcp_bits = str2num(plcp_bits);
plcp_rate = str2num(plcp_rate);
preamble_bits = str2num(preamble_bits);

% For formatting purposes for SNR values
tokenized = regexp(snr_string,',', 'split');

% Take SNR values from parameter and place them in an array
set_of_snrs = [];
for i = 1:length(tokenized)
   char(tokenized(i))
   set_of_snrs(i) = str2num(char(tokenized(i)));
end

% For formatting purposes for rates
tokenized = regexp(rates_string,',', 'split');

% Take rates from parameter and place them in an array
set_of_rates = [];
for i = 1:length(tokenized)
   char(tokenized(i))
   set_of_rates(i) = str2num(char(tokenized(i)));
end

% Open SNRTable.txt file for appending
file = fopen('SNRTable.txt','a');
    
% Write the SNR values on the first line of the file
for i = 1:length(set_of_snrs)
    fprintf(file, '%d ', set_of_snrs(i));
end

fprintf(file, '\n');

% Open RateProb.txt file for appending
fileid = fopen('RateProb.txt','a');
    
% Write the rates on the first line of the file
for i = 1:length(set_of_rates)
    fprintf(fileid, '%d ', set_of_rates(i));
end    
    
fprintf(fileid, '\n');

% Define the number of iterations and the bandwidth
num_iterations = 1000000;
bandwidth = 20e6;

% Loop through each SNR value to generate distributions and determine
% goodput optimal rate for payload ranges
for j = 1:length(set_of_snrs)
    
    % Convert SNR to pure ratio (signal power/noise power)
    % Initialize a vector for capacity value
    snr = 10^(set_of_snrs(j)/10);
    C = zeros(1,num_iterations);

    % Loop for the number of iterations
    % Step 1: Generate H matrix
    % Step 2: Compute the SVD
    % Step 3: Generate new H matrix
    % Step 4: Calculate capacity and store value
    for i = 1:num_iterations
        H = (1/sqrt(2))*randn(num_receive, num_transmit) + (1/sqrt(2))*1i*randn(num_receive, num_transmit);
        [A,D,B] = svd(H);
        V = B(:,1:num_streams);
        H_eff = H*V;

        C(i) = real(bandwidth*log2(det(eye(num_receive) + (snr/num_streams)*H_eff*V'*H')))/1e6;
        disp(sprintf('%d SNR: Iteration %d of %d', snr, i, num_iterations));
    end


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Generate PDF and CDF
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
    % Note: Graphs will be plotted. If figures are not needed, the
    %       appropriate lines can be commented out.
    x_range = min(C):0.1:max(C); % bins originally width of 0.5
    figure;
    subplot(2,1,1)
    hist(C)
    num_elements = histc(C,x_range);
    pdf = num_elements/num_iterations; 
    plot(x_range,pdf)
    grid on
    string = sprintf('Probability Density Function for Capacity of %d x %d MIMO Link, %d streams, %d dB SNR', num_transmit, num_receive, num_streams, set_of_snrs(j));
    ylabel('Probability')
    xlabel('Capacity (Mbps)')
    title(string);
    subplot(2,1,2)
    cumulative_elements = cumsum(num_elements);
    cdf = cumulative_elements/num_iterations;
    plot(x_range, cdf)
    grid on
    string_two = sprintf('Cumulative Density Function for Capacity of %d x %d MIMO Link, %d streams, %d dB SNR', num_transmit, num_receive, num_streams, set_of_snrs(j));
    title(string_two);
    ylabel('Probability')
    xlabel('Capacity (Mbps)')

    % Generate rate vector with rates sorted and converted to bps units
    sorted_rate_vector = sort(set_of_rates)*10^6;
    num_of_rates = length(sorted_rate_vector);
    cont_time_rates = x_range*10^6;
    count = 0;
    index = 1;

    % Initialize vector for packet fail probabilities
    associated_prob_fail = [];

    % For each rate, loop through and determine the probability of a packet
    % transmission failing by using the distributions generated.
    while count ~= num_of_rates     
       matching_rate_index = find(cont_time_rates < sorted_rate_vector(index),1,'last');

       if length(matching_rate_index) ~= 0
         associated_prob_fail(index) = cdf(matching_rate_index);
       else
         associated_prob_fail(index) = 0;
       end

       if matching_rate_index == length(cdf)
           associated_prob_fail(index) = 1;
       end

       index = index + 1;
       count = count + 1;
    end

    % For debugging purposes, print out the associated probabilities of
    % failure and success. Comment out if not needed.
    format long
    associated_prob_fail
    success = 1-associated_prob_fail
    
    % Open RateProb.txt for appending. Print out SNR value this
    % distribution was generated for.
    fileid = fopen('RateProb.txt','a');
    fprintf(fileid, '%d ', set_of_snrs(j));
    
    % Print out the associated probabilities of success for each rate for
    % this distribution.
    for i = 1:length(success)
        fprintf(fileid, '%f ', success(i));
    end    
    
    fprintf(fileid, '\n');

    % Calculate the fixed cost transmission time
    fixed_cost = 0;
    plcp = plcp_bits/plcp_rate;
    difs = sifs + 2*slot_time;
    ack = (14 * 8)/basic_rate;
    backoff = (cw * slot_time)/2;
    preamble = preamble_bits/basic_rate;
    
    % Fixed cost. ACK is included in total transmission time calculation.
    fixed_cost = (difs + backoff + plcp + sifs + plcp)*10^-6;

       
    % MAC HEADER = 30 bytes
    % PAYLOAD = 0-2312 bytes
    % CRC = 4 bytes
    % Therefore, range of bytes for entire packet is 0 to 2346
    bits = 34*8:8:2346*8;
    bits_payload = 0:8:2312*8;
    gp = []; 
    
    
    % For each rate, loop through and determine the total transmission time
    % and the goodput for payload bits.
    for k = 1:length(sorted_rate_vector)
        
        % Total transmission time
        transmit_time = fixed_cost + bits/(sorted_rate_vector(k)) + ack*10^-6;
        
        if success(k) == 0
            gp = vertcat(gp, 0*bits_payload./transmit_time);
        else
            gp = vertcat(gp, bits_payload./(transmit_time*(1/success(k))));
        end
        
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    % Loop through and for each payload size determine which rate provides
    % the maximum gooodput. This allows us to determine which rate is
    % optimal for each range
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [maxgp previndex] = max(gp(:,[1]));
    row = 1;
    rate = sorted_rate_vector(previndex)/10^6;
    
    ranges = [];
    ranges(row,1) = rate;
    ranges(row,2) = bits_payload(1)/8;    
    [r,c] = size(gp);
    
    for i=2:c
        [maxgp, index]  = max(gp(:,[i]));
        
        if index ~= previndex
            ranges(row,3) = bits_payload(i)/8 - 1;
            row=row+1;

            rate = sorted_rate_vector(index)/10^6;

            ranges(row,1) = rate;
            ranges(row,2) = bits_payload(i)/8;
            previndex = index;
        end

    end

    [r,c] = size(ranges);
    ranges(r,3) = bits_payload(length(bits_payload))/8;
    
    if (ranges(1,2) == 0 && ranges(1,3) == 0)
        ranges(1,:)=[];
        [r,c] = size(ranges);
    end
    
    % Print ranges for debugging purposes. Can be uncommented if not needed.
    ranges
    
    % Open SNRTable.txt for appending
    file = fopen('SNRTable.txt','a');
    
    % Loop through for each range and print the rate and range values to
    % the fil.
    for i = 1:r
        fprintf(file, '%d,%d-%d ', ranges(i,1), ranges(i,2), ranges(i,3));
    end

    fprintf(file, '\n');
    

end


end