mimo_detection.m

%% 
% ML-ZF-MMSE 检测算法比较
% bpsk调制-瑞利衰落信道-2个发送端，2个接收端

%%
clc
clear
N = 10^5; % 发送的符号数目
Eb_N0_dB = 0:20; % 信噪比范围
nTx = 2;
nRx = 2;
%%
for ii = 1:length（Eb_N0_dB）

     % 发送端
     ip = rand（1N）>0.5; % 等概率产生0和1
     s = 2*ip-1; % BPSK 调制 0 -> -1; 1 -> 1
     sMod = kron（sones（nRx1））; % 
     sMod = reshape（sMod[nRxnTxN/nTx]）; % 将矩阵转换为[nRxnTxN/nTx ]形式 

     h = 1/sqrt（2）*[randn（nRxnTxN/nTx） + j*randn（nRxnTxN/nTx）]; % 瑞利衰落信道
     n = 1/sqrt（2）*[randn（nRxN/nTx） + j*randn（nRxN/nTx）]; % 0均值高斯白噪声


     % 加入噪声后在信道中传输

     y = squeeze（sum（h.*sMod2）） + 10^（-Eb_N0_dB（ii）/20）*n;
%%
%ML
     % 最大似然接收
     % ----------------------------
     % 当 [s1 s2 ] = [+1+1 ]
     sHat1 = [1 1];
     sHat1 = repmat（sHat1[1 N/2]）;
     sHat1Mod = kron（sHat1ones（nRx1））;
     sHat1Mod = reshape（sHat1Mod[nRxnTxN/nTx]）;%发送矩阵
     zHat1 = squeeze（sum（h.*sHat1Mod2）） ; %接收矩阵
     %e = abs（y - zHat1）; %取理论接收与实际接收的绝对值
     J11 = sum（abs（y - zHat1）1）;%将两行加为一行

     % 当 [s1 s2 ] = [+1-1 ]
     sHat2 = [1 -1];
     sHat2 = repmat（sHat2[1 N/2]）;
     sHat2Mod = kron（sHat2ones（nRx1））;
     sHat2Mod = reshape（sHat2Mod[nRxnTxN/nTx]）;
     zHat2 = squeeze（sum（h.*sHat2Mod2）） ;
     J10 = sum（abs（y - zHat2）1）;

     % 当 [s1 s2 ] = [-1+1 ]
     sHat3 = [-1 1];
     sHat3 = repmat（sHat3[1 N/2]）;
     sHat3Mod = kron（sHat3ones（nRx1））;
     sHat3Mod = reshape（sHat3Mod[nRxnTxN/nTx]）;
     zHat3 = squeeze（sum（h.*sHat3Mod2）） ;
     J01 = sum（abs（y - zHat3）1）;

     % 当 [s1 s2 ] = [-1-1 ]
     sHat4 = [-1 -1];
     sHat4 = repmat（sHat4[1 N/2]）;
     sHat4Mod = kron（sHat4ones（nRx1））;
     sHat4Mod = reshape（sHat4Mod[nRxnTxN/nTx]）;
     zHat4 = squeeze（sum（h.*sHat4Mod2）） ;
     J00 = sum（abs（y - zHat4）1）;

     % 从四组数值中找出最小值jj，所在位置dd
     rVec = [J11;J10;J01;J00];
     [jj dd] = min（rVec[]1）;

     % mapping the minima to bits
     ref = [1 1; 1 0; 0 1; 0 0 ];
     ipHat = zeros（1N）;
     ipHat（1:2:end） = ref（dd1）;
     ipHat（2:2:end） = ref（dd2）;

     % 统计错误
     f = find（[ip - ipHat]）;%发生错误所在位置
     nErrML（ii） = size（find（[ip - ipHat]）2）;%错误个数


%%
%ZF
     % 接收端

    % 求伪逆矩阵 G = inv（H^H*H）*H^H
    % H^H*H 的矩阵维数为 [nTx x nTx]. 即 [2 x 2] 
    % [2x2]矩阵的转置为 [a b; c d] = 1/（ad-bc）[d -b;-c a]

    hCof = zeros（22N/nTx）; 
    hCof（11:） = sum（h（:2:）.*conj（h（:2:））1）; % d 
    hCof（22:） = sum（h（:1:）.*conj（h（:1:））1）; % a 
    hCof（21:） = -sum（h（:2:）.*conj（h（:1:））1）; % c 
    hCof（12:） = -sum（h（:1:）.*conj（h（:2:））1）; % b 
    hDen = （（hCof（11:）.*hCof（22:）） - （hCof（12:）.*hCof（21:）））; % ad-bc
    hDen = reshape（kron（reshape（hDen1N/nTx）ones（22））22N/nTx）; % 
    hInv = hCof./hDen; % inv（H^H*H）

    hMod = reshape（conj（h）nRxN）; % H^H 

    yMod = kron（yones（1