simulate a MIMO spatial multiplexing system using MATLAB

% Simulation parameters
Nt = 2;                 % Number of transmit antennas
Nr = 2;                 % Number of receive antennas
modulationOrder = 4;    % QPSK modulation
numBits = 1e6;          % Total number of bits
SNR_dB = 0:5:30;        % SNR values in dB
ber = zeros(size(SNR_dB)); % Store BER results

% Validate that numBits is divisible by Nt * log2(modulationOrder)
bitsPerSymbol = log2(modulationOrder);
numSymbolsPerStream = numBits / (Nt * bitsPerSymbol);
if mod(numBits, Nt * bitsPerSymbol) ~= 0
    error('numBits must be divisible by Nt * log2(modulationOrder).');
end

% Generate random bits
bits = randi([0 1], numBits, 1);

% Split bits into Nt streams
bitsStreams = reshape(bits, bitsPerSymbol, numSymbolsPerStream, Nt);
bitsStreams = permute(bitsStreams, [2 1 3]); % numSymbolsPerStream x bitsPerSymbol x Nt

% Modulate each stream
txSymbols = zeros(numSymbolsPerStream, Nt);
for i = 1:Nt
    txSymbols(:, i) = qammod(bitsStreams(:,:,i), modulationOrder, ...
        'InputType', 'bit', 'UnitAveragePower', true);
end

% Normalize transmit power by sqrt(Nt)
txSymbols = txSymbols / sqrt(Nt);

% Simulate for each SNR value
for snrIdx = 1:length(SNR_dB)
    SNR = SNR_dB(snrIdx);
    SNR_linear = 10^(SNR/10);
    N0 = 1 / SNR_linear;        % Noise variance
    numErr = 0;                 % Reset error counter
    
    % Process each symbol
    for symIdx = 1:numSymbolsPerStream
        % Transmitted symbols for current symbol time
        X = txSymbols(symIdx, :).'; % Nt x 1
        
        % Generate Rayleigh fading channel matrix
        H = (randn(Nr, Nt) + 1i*randn(Nr, Nt)) / sqrt(2); % CN(0,1)
        
        % Generate complex noise
        noise = sqrt(N0/2) * (randn(Nr, 1) + 1i*randn(Nr, 1));
        
        % Received signal
        Y = H * X + noise;
        
        % Zero Forcing detection
        H_inv = pinv(H);            % Moore-Penrose pseudo-inverse
        X_hat = H_inv * Y;          % Estimate symbols
        X_est = X_hat * sqrt(Nt);   % Rescale to original power
        
        % Demodulate each stream
        rxBits = zeros(bitsPerSymbol, Nt);
        for i = 1:Nt
            rxBits(:, i) = qamdemod(X_est(i), modulationOrder, ...
                'OutputType', 'bit', 'UnitAveragePower', true);
        end
        rxBits = rxBits(:); % Flatten to column vector
        
        % Original bits for this symbol
        originalBits = bitsStreams(symIdx, :, :);
        originalBits = originalBits(:);
        
        % Count errors
        numErr = numErr + sum(rxBits ~= originalBits);
    end
    
    % Calculate BER for current SNR
    ber(snrIdx) = numErr / numBits;
    fprintf('SNR = %d dB, BER = %.4f\n', SNR, ber(snrIdx));
end

% Plot BER vs SNR
figure;
semilogy(SNR_dB, ber, 'bo-', 'LineWidth', 2);
grid on;
xlabel('SNR (dB)');
ylabel('Bit Error Rate (BER)');
title('2x2 MIMO Spatial Multiplexing with Zero Forcing');