Untitled

clear all;
close all;
clc;
runid = 0;
rng(runid);

% Meta parameters
par.SNR_dB_list = 0:30;
par.T = 1e4;

% I am assuming that we use a = 1

par.mod = 'BPSK';
par.channel = 'SIMO_flat_rayleigh_optimal';
par.SIMO = [1,2,3,4];

res_SIMO_flat_rayleigh = channel_simulation(par);

par.mod = 'BPSK';
par.channel = 'SIMO_flat_rayleigh_selection_combining';
par.SIMO = [1,2,3,4];

res_SIMO_flat_rayleigh_selection_combining = channel_simulation(par);


for i=1:length(par.SIMO)
    semilogy(par.SNR_dB_list,res_SIMO_flat_rayleigh.P_e(:,i));
    hold on;
    semilogy(par.SNR_dB_list,res_SIMO_flat_rayleigh_selection_combining.P_e(:,i),'o-');
end

grid on
axis([0 20 1e-6 1]);
xlabel('average receive SNR [dB]');
ylabel('symbol error rate (SER)');
legend('M_r = 1', 'SC, M_r = 1', 'M_r = 2', 'SC, M_r = 2', 'M_r = 3', 'SC, M_r = 3', 'M_r = 4', 'SC, M_r = 4');
title('SIMO Optimal Detection and Selection Combining Comparison');

function symbols = constellation_set(par)
    if par.mod == "QPSK"
        symbols = [1+1i, -1+1i, 1-1i, -1-1i];
    elseif par.mod == "BPSK"
        symbols = [1, -1];
    end
end

function x_est = decoder(y, par, h)
    if par.mod == "QPSK"
        x_est = sign(real(y)) + sign(imag(y))*1i;
    elseif par.mod == "BPSK" && isfield(par, 'SIMO')
        if (~exist('h', 'var'))
            error("SIMO flat rayleigh channel estimator requires h");
        end
        x_est = SIMO_est_BPSK(y, h);
    elseif par.mod == "BPSK"
        x_est = sign(real(y));
    end
end


function res = channel_simulation(par)
    SNR = 10.^(par.SNR_dB_list / 10);
    symbols = constellation_set(par);
    E_s = mean(abs(symbols).^2);
    
    var_N0 = E_s ./ SNR;
    x = symbols(randi([1, length(symbols)], 1, par.T));
    x_est = zeros(length(par.SNR_dB_list), par.T);
    
    if par.channel == "SIMO_flat_rayleigh_optimal" || par.channel == "SIMO_flat_rayleigh_selection_combining"
        if ~isfield(par, 'SIMO')
            error("SIMO flat rayleigh channel requires M_r values in par.SIMO");
        end
        
        res.P_e = zeros(length(par.SNR_dB_list), length(par.SNR_dB_list));
        
        for kk=1:length(par.SNR_dB_list)
            for i=1:length(par.SIMO)
                h = sqrt(0.5)*( randn(par.SIMO(i), par.T) + 1i * randn(par.SIMO(i), par.T));
                n = sqrt((var_N0(kk) * 0.5)) * (randn(par.SIMO(i), par.T) + 1i * randn(par.SIMO(i), par.T));
                y = h .* kron(ones(par.SIMO(i), 1), x) + n;
                
                if par.channel == "SIMO_flat_rayleigh_selection_combining"
                    [~,idx_h_max] = max(abs(h),[],1);
                    y = y(sub2ind(size(y), idx_h_max, 1:par.T));
                    h = h(sub2ind(size(h), idx_h_max, 1:par.T));
                end
                
                x_est = decoder(y, par, h);
                res.P_e(kk,i) = sum(x ~= x_est) / par.T;
            end
        end
    elseif par.channel == "SIMO_flat_rayleigh_selection_combining"
        if ~isfield(par, 'SIMO')
            error("SIMO flat rayleigh channel requires M_r values in par.SIMO");
        end
        
        res.P_e = zeros(length(par.SNR_dB_list), length(par.SNR_dB_list));
        
        for kk=1:length(par.SNR_dB_list)
            for i=1:length(par.SIMO)
                h = sqrt(0.5)*( randn(par.SIMO(i), par.T) + 1i * randn(par.SIMO(i), par.T));
                n = sqrt((var_N0(kk) * 0.5)) * (randn(par.SIMO(i), par.T) + 1i * randn(par.SIMO(i), par.T));
                y = h .* kron(ones(par.SIMO(i), 1), x) + n;
                
                x_est = decoder(y, par, h);
                
                res.P_e(kk,i) = sum(x ~= x_est) / par.T;
            end
        end
    elseif par.channel == "AWGN"
        n = sqrt((var_N0 ./ 2))' * (randn(1, par.T) + 1i * randn(1, par.T));
        y = x + n;
        x_est = decoder(y, par);
        
        res.P_e = mean(x ~= x_est, 2); % average over trials (axis 2)
    end
end
    
function x_est = SIMO_est_BPSK(y, h)
    % y in R^M_r x T, h in R^M_r x T
    h_conj = conj(h);
    y_combined = sum(h_conj .* y, 1);  
    x_est = sign(real(y_combined));
end