detection.m

MIMO2by2_16QAM_SoIC+Turbo迭代检测simulink链路

function [ Output ] = Detection( Input,snr,ITERATION_SoIC )
% MIMO Detection
% 

warning off all

TX = 2;
M = 16;% 16QAM

% 初始符号概率
% Yita = [0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625;
%         0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625];
% Yita = [0.25 0.25 0.25 0.25;     % 第一根天线上的发送00,01,10,11的概率
%         0.25 0.25 0.25 0.25];    % 第二根天线上的发送00,01,10,11的概率

Constellation = [-3 + 3i;-3 + 1i;-3 - 3i;-3 - 1i;-1 + 3i;-1 + 1i;-1 - 3i;-1 - 1i;  % 0000,0001,0010,0011,0100,0101,0110,0111
                  3 + 3i;3 + 1i;3 - 3i;3 - 1i;1 + 3i;1 + 1i;1 - 3i;1 - 1i]/sqrt(10)/sqrt(2);        % 1000,1001,1010,1011,1100,1101,1110,1111
% Constellation = [-1 + 1i;-1 - 1i;1 + 1i;1 - 1i];  % 00,01,10,11
Rvv = diag([1/(10^(snr/10)),1/(10^(snr/10))]);

[m,n] = size(Input);

for counter = 1:m/2
    % Divide input signals
    H = Input((2*(counter-1)+1):(2*(counter-1)+2),1:TX);
    ReceivedSignal = Input((2*(counter-1)+1):(2*(counter-1)+2),TX+1);
    Yita = Input((2*(counter-1)+1):(2*(counter-1)+2),(TX+2):(TX+1+M));

%     if 1%sum(sum(Yita)) < 1.5
%         Yita = [0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625;
%                 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625];
%     end

    for ITERATION_SoICNum = 1:ITERATION_SoIC
        % 计算均值、功率、方差
        Mean = Yita*Constellation;
        Power = Yita*(Constellation.*conj(Constellation));
        Var = Power - Mean.*conj(Mean);

        % 抵消干扰
        for i = 1:TX
            TempH = H;
            TempH(:,i) = 0;
            TempSignal = ReceivedSignal - TempH*Mean;
            Signal(:,i) = TempSignal;
        end

        W = diag(Var);
        prb = zeros(TX,M);
        DataforCompare = zeros(1,M);
        for i = 1:TX
            % MMSE滤波
            TempW = W;
            TempW(i,i) = Power(i);
            w_MMSE = Power(i)*inv(H*TempW*H'+Rvv)*H(:,i);
            
            % 计算符号概率
            % 求Yita的简化算法（此简化算法不损失性能）
            Miu = real(w_MMSE'*H(:,i));
            Delta = (Miu-Miu*Miu)*Power(i);
            for j = 1:M
                DataTemp = abs(w_MMSE'*(Signal(:,i)-H(:,i)*Constellation(j)));
                DataforCompare(j) = DataTemp*DataTemp;
            end
            Min_DataforCompare = min(DataforCompare);
            for j = 1:M
                prb(i,j) = exp((Min_DataforCompare-DataforCompare(j))/Delta);
            end
            prbsum = sum(prb(i,:));
            if isnan(prbsum) == zeros(TX)
                Yita(i,:) = prb(i,:)/prbsum;
            end
%             % 求Yita的非简化算法
%             S(:,i) = w_MMSE'*Signal(:,i);
%             Miu = real(w_MMSE'*H(:,i));
%             Delta = (Miu-Miu*Miu)*Power(i);
%             for j = 1:M
%                 Yita(i,j) = (1/pi/Delta)*exp(-abs(Miu*Constellation(j)-S(:,i))^2/Delta);
%             end
        end
    end

    % 计算比特似然比
    blr = [(Yita(:,9)+Yita(:,10)+Yita(:,11)+Yita(:,12)+Yita(:,13)+Yita(:,14)+Yita(:,15)+Yita(:,16))./(Yita(:,1)+Yita(:,2)+Yita(:,3)+Yita(:,4)+Yita(:,5)+Yita(:,6)+Yita(:,7)+Yita(:,8)),...
           (Yita(:,5)+Yita(:,6)+Yita(:,7)+Yita(:,8)+Yita(:,13)+Yita(:,14)+Yita(:,15)+Yita(:,16))./(Yita(:,1)+Yita(:,2)+Yita(:,3)+Yita(:,4)+Yita(:,9)+Yita(:,10)+Yita(:,11)+Yita(:,12)),...
           (Yita(:,3)+Yita(:,4)+Yita(:,7)+Yita(:,8)+Yita(:,11)+Yita(:,12)+Yita(:,15)+Yita(:,16))./(Yita(:,1)+Yita(:,2)+Yita(:,5)+Yita(:,6)+Yita(:,9)+Yita(:,10)+Yita(:,13)+Yita(:,14)),...
           (Yita(:,2)+Yita(:,4)+Yita(:,6)+Yita(:,8)+Yita(:,10)+Yita(:,12)+Yita(:,14)+Yita(:,16))./(Yita(:,1)+Yita(:,3)+Yita(:,5)+Yita(:,7)+Yita(:,9)+Yita(:,11)+Yita(:,13)+Yita(:,15))];     % 列：第几bit，行：第几天线
    % 0000,0001,0010,0011,0100,0101,0110,0111,1000,1001,1010,1011,1100,1101,1110,1111

    out = zeros(8,1);
    for i = 1:TX
        for j = 1:log2(M)
            if isinf(blr(i,j)) == 1     % 输出限幅
                blr(i,j) = 1e50;
            elseif isnan(blr(i,j)) == 1
                blr(i,j) = 1e-50;
            elseif blr(i,j) > 1e50
                blr(i,j) = 1e50;
            elseif blr(i,j) < 1e-50
                blr(i,j) = 1e-50;
            end
            blr(i,j) = log(blr(i,j));
            out((i-1)*log2(M)+j,1) = blr(i,j);% 将输出比特似然比的维数变为((TX-1)*2+log2(M))*1
        end
    end
    
    outtemp((8*(counter-1)+1):(8*(counter-1)+8),1) = out;
end

% 输出信号
Output = outtemp;