基于QPSK调制的多天线无编码通信系统中ZF、MMSE、SIC、ML符号检测性能比较。

一、题目

        输入数据为猫图，并行分成5个数据流，每个都是QPSK调制，分别通过5个天线并行发送，单经瑞丽衰落信道（5x5信道中每个元素皆为CN(0,1)分布）外加AWGN，接收端为5个接收天线。比较ZF，MMSE，SIC，ML符号检测的BER性能。（无信道信源编码）

1. 画出0:10:30dB输出猫图
2. 在一个图中，画出0:10:30dB 4种算法的BER性能
3. 在一个图中，画出4种算法仿真时间。分析时间复杂度与BER性能的关系。

二、题目分析

1、算法原理

        （1）ZF（零离差）算法：

        ZF 算法是一种线性检测算法，其核心思想是通过自适应滤波的方式减少通道的影响。ZF 算法中，需要对传输通道进行估计和反演，从而得到接收信号的估计值，然后根据估计值来计算符号的似然概率，选择具有最大似然概率的符号作为输出。

        （2）MMSE（最小均方误差）算法：

        MMSE 算法是一种非线性检测算法，其核心思想是通过最小化符号估计值与实际值之间的均方误差，提高检测性能。MMSE 算法中，需要通过先验概率分布和后验概率分布来优化符号的估计过程，从而获得更准确的符号估计值。

        （3）SIC（逐步干扰消除）算法：

        SIC 算法是一种迭代检测算法，其核心思想是通过消除多径干扰，提高符号检测精度。SIC 算法中，首先需要对接收信号进行多次干扰消除，然后根据估计值得到符号的似然概率，选择具有最大似然概率的符号作为输出。

        （4）ML（最大似然）符号检测算法：

        ML 算法是一种基于统计估计的非线性检测算法，它通过最大化符号的似然概率来实现符号检测。ML 算法中，需要对所有可能的符号序列进行枚举，计算其似然概率，然后选择具有最大似然概率的符号序列作为输出。

2、实现步骤

        题目总体上要求我们通过仿真比较ZF、MMSE、SIC、ML四种符号检测算法在不同信噪比下的BER性能和运行时间。具体步骤如下：

        （1）定义仿真参数:

        我们需要定义仿真所需的各种参数，包括调制方式、信道、信噪比、数据流数、天线数等。对于该题目，调制方式为 QPSK，信道为复高斯随机数信道，天线数为 5，数据流数为 5。

​​​​​​​        （2）生成随机数据:

        我们需要生成测试所需的随机猫图数据，该题目要求将数据并行分成 5 个数据流，每个流再进行 QPSK调制。

​​​​​​​        （3）发送信号:

        我们需要将生成的数据流通过 5 个天线并行发送，并经过单经瑞利衰落信道 (信道系数为 5x5 CN(0,1) 噪声矩阵)，再加入 AWGN 噪声。

        （4）​​​​​​​接收信号:

        我们需要在接收端使用 5 个接收天线接收到发射的信号，并对接收到的信号进行复信道估计，以便后续的符号检测。（为了更简便地看到四种符号检测算法的性能区别，这一步假定信道已经做了精准的估计，所以直接使用原有的信道数据）

        （5）​​​​​​​符号检测:

         我们需要使用四种不同的符号检测算法 (ZF, MMSE, SIC, ML) 对接收到的数据进行解调和符号检测，以恢复原始猫图数据。

        （6）​​​​​​​计算误差率:

        我们需要将检测错误的比特数除以总比特数，计算各个算法在不同信噪比下的误比特率 (BER) 。   

        （7）​​​​​​​绘制输出猫图:

        我们需要将模拟的猫图数据绘制输出，以便观察模拟效果。

        （8）​​​​​​​绘制BER性能曲线:

        我们需要将四种符号检测算法在不同信噪比下的误比特率 (BER) 计算结果，绘制 BER 性能曲线，以方便比较四种算法的性能。

​​​​​​​        （9）计算仿真时间:

        我们需要在进行第6步的同时记录每种符号检测算法的仿真时间，并根据信噪比和数据流数的不同进行比较分析。

三、仿真结果 

​​​​​​​​​​​​​​​​​​​​​

(a) ZF                        (b) MMSE                       (c) SIC                  (d) ML

图 1 不同检测算法与不同信噪比下输出猫图

图 2 四种检测算法的输出误码率比较

图 3 四种检测算法的仿真时间对比

四、结果分析

        本次仿真主要比较ZF、MMSE、SIC、ML四种符号检测算法的性能。因为在信噪比SNR大于10dB之后，四种算法的误码率基本趋于0，故将信噪比范围限定在0至10dB,观察输出的各种指标。

        图1(a)(b)(c)(d)分别是ZF、MMSE、SIC、ML四种算法在0至10dB情况下输出的猫图。

        图2是四种检测算法的输出误码率比较。由图2可以看出，在任何信噪比情况下，ML算法都是四种算法中性能最好的，因为其牺牲了算法的复杂度来换取高性能，由图3可以看出ML是这四种算法中时间复杂度最高的一种。

        由图2还可以看出，当信噪比较低的时候，非线性检测算法 MMSE 比线性检测算法 ZF 和 SIC 更具优势。这是因为在低信噪比下，信号受到的噪声干扰非常大，而线性检测算法无法完全抵消这些噪声干扰。然而，非线性检测算法能够利用先验的信号统计信息估计符号的先验概率分布，从而更好地抵抗噪声干扰。

        而随着信噪比的提高，SIC 算法相对于 MMSE 算法在性能表现方面会出现逐渐超越的情况，这主要是因为SIC 算法可以通过多次干扰消除的方式来消除干扰（从图3也可以看出其正因为多次干扰消除的步骤，使得时间复杂度相对于ZF和MMSE稍高），从而最大化接收信号的质量。在高信噪比的情况下，SIC算法能够更好地将有效信号与噪声和干扰分离开来，减少信号受到噪声和干扰的影响，从而提高检测性能。

        综上所述，在使用不同的符号检测算法时，需要根据具体的情况进行选择和权衡，并进行适当的优化和调整，以实现最佳的性能和效果。

五、代码

clear ;
close all;

% 生成QPSK调制的发送数据
M = 4; % 调制阶数
k = log2(M); % 每个符号所需比特数
img = imread('coolLearn.PNG');              % 读取彩色图片
data = reshape(img, [], 1);                 % 将图像数据展开为列向量
bData = dec2bin(data,8);                %转换为8位二进制数
bData1 = reshape(bData, [], k);         %将二进制数分为4组，每组2位
decData = bin2dec(bData1);              %将二进制转成十进制0-3
% 序列补位
N = length(decData);                    % 获得序列元素个数
zero_tail = mod(5-mod(N,5),5);          % 获得需要对输入序列进行补位的个数
decData_0 = [decData;zeros(zero_tail,1)]; % 末尾以0补位

s_qam = qammod(decData_0, M, 'gray');       % QPSK调制发送数据   
X = reshape(s_qam, 5, []);                  % 将发送数据并行分成5个流

% 设置信道参数
N_TX = 5;        % 发送天线数
N_RX = 5;        % 接收天线数
H = (randn(N_RX, N_TX) + 1j * randn(N_RX, N_TX)) / sqrt(2); % 生成信道矩阵
SNRdB = 0:2:10;       % 信噪比范围（单位：dB）
% 初始化BER和time结果变量
berZF = zeros(length(SNRdB), 1);
berMMSE = zeros(length(SNRdB), 1);
berSIC = zeros(length(SNRdB), 1);
berML = zeros(length(SNRdB), 1);
timeZF = zeros(length(SNRdB), 1);
timeMMSE = zeros(length(SNRdB), 1);
timeSIC = zeros(length(SNRdB), 1);
timeML = zeros(length(SNRdB), 1);

%% 开始仿真过程
for i = 1:length(SNRdB)
    % 生成高斯白噪声
    sigma = sqrt(10^(-SNRdB(i)/10)); % 信噪比计算公式
    n = sigma*(randn(N_RX,length(X))+1i*randn(N_RX,length(X)))/sqrt(2);   
    % 生成接收信号
    Y = H*X + n;
    %% ZF符号检测
    tstart = tic; % 计算运行时间   
    tZF = pinv(H.'*H)* H.' *Y;      % ZF检测
    tZF1 = reshape(tZF, [], 1);     % 重构为向量
    sZF = qamdemod(tZF1, M, 'gray');% QPSK解调
    toc(tstart);
    timeZF(i) = toc(tstart); % 记录运行时间

    [~,berZF(i)] = biterr(decData_0 , sZF)% 计算ZF误码率
    
    sZF = sZF(1: end-zero_tail);    %去掉补位的零
    ZFb1 = dec2bin(sZF,2);
    ZFb2 = reshape(ZFb1,[],8);
    ZFd1 = bin2dec(ZFb2);           %二进制转十进制
    ZFd2 = reshape(ZFd1,size(img))/255;
    
    figure(1)
    subplot(length(SNRdB),1,i);
    imshow(im2uint8(ZFd2));%画猫图
    title(['SNRdB=',num2str(SNRdB(i)),'dB']);

    %% MMSE符号检测
    tstart = tic; % 计算运行时间
    gamma = 10 ^ (- SNRdB(i) / 10); %设置信噪比参数
    LAMBDA = gamma * eye(N_RX); % 计算正比例参数
    wMMSE = H' * pinv(H * H' + LAMBDA);%计算MMSE检测权重
    tMMSE = wMMSE * Y;%进行MMSE检测
    tMMSE1 = reshape(tMMSE, [], 1);%重构为向量
    sMMSE = qamdemod(tMMSE1, M, 'gray');%QAM解调
    toc(tstart);
    timeMMSE(i) = toc(tstart); % 记录运行时间

    [~,berMMSE(i)] = biterr(decData_0 , sMMSE)%计算MMSE误码率
   
    sMMSE = sMMSE(1: end-zero_tail);%去掉补位的零
    MMSEb1 = dec2bin(sMMSE,2);
    MMSEb2 = reshape(MMSEb1,[],8);
    MMSEd1 = bin2dec(MMSEb2);%二进制转十进制
    MMSEd2 = reshape(MMSEd1,size(img))/255; 
    
    figure(2)
    subplot(length(SNRdB),1,i);
    imshow(im2uint8(MMSEd2));%画猫图
    title(['SNRdB=',num2str(SNRdB(i)),'dB']);

    %% SIC检测
    tstart = tic; % 计算运行时间
    squares_sum = sum(abs(H) .^2, 1);% 计算H每列的平方和
    [~, sort_idx] = sort(squares_sum);% 对平方和进行排序，并返回排序后的索引
    H_sorted = H(:,sort_idx);  % 按索引顺序排序H每列
    Y_temp = Y ;              
    tSIC = zeros(size(X));     %给还原的数据预分配一个数组空间
    for sic_i = 1:N_TX
        X_temp = pinv(H_sorted.'*H_sorted)* H_sorted.' * Y_temp;   % ZF
        tSIC(end-sic_i+1,:) = qamdemod(X_temp(end,:),4); % 取ZF功率最大的一行，解调
        Y_temp = Y_temp - H_sorted(:,end)*qammod(tSIC(end-sic_i+1,:),4); % 更新Y_temp
        H_sorted = H_sorted(:,1:end-1);% 删除已估计的列
    end
    % 恢复原始矩阵行数顺序
    [~, original_idx] = sort(sort_idx);
    tSIC = tSIC(original_idx, :);
    toc(tstart);
    timeSIC(i) = toc(tstart); % 记录运行时间
    sSIC = reshape(tSIC, [], 1);%sic过程中已经解调,直接重构为向量

    [~,berSIC(i)] = biterr(decData_0 , sSIC)%计算SIC误码率

    sSIC = sSIC(1: end-zero_tail);%去掉补位的零
    SICb1 = dec2bin(sSIC,2);
    SICb2 = reshape(SICb1,[],8);
    SICd1 = bin2dec(SICb2);%二进制转十进制
    SICd2 = reshape(SICd1,size(img))/255;
    
    figure(3)
    subplot(length(SNRdB),1,i);
    imshow(im2uint8(SICd2));%画猫图
    title(['SNRdB=',num2str(SNRdB(i)),'dB']);

    %% ML检测
    tstart = tic; % 计算运行时间
    permutation_list = (dec2base(0:M^N_TX-1, 4) - '0').'; % 生成所有可能的符号排列
    permutation_list_qam = qammod(permutation_list,4);  % 对可能的符号进行QPSK调制
    num_permutations = size(permutation_list,2);        % 统计有多少种符号排序
    distance = zeros(1,num_permutations);       %初始化一个距离矩阵
    tML = zeros(size(X));                       %给还原的数据预分配一个数组空间
    for ml_i = 1:length(Y)
        ml_xk_matrix = H*permutation_list_qam;
        distance = sum(abs(Y(:,ml_i) - ml_xk_matrix).^2);
        [~, ml_index] = min(distance);
        tML(:,ml_i) = permutation_list(:,ml_index);
    end
    toc(tstart);
    timeML(i) = toc(tstart); % 记录运行时间
    sML = reshape(tML, [], 1);              %ML过程已解调，直接重构为向量
    
    [~,berML(i)] = biterr(decData_0 , sML)  %计算ML误码率

    sML = sML(1: end-zero_tail);            %去掉补位的零
    MLb1 = dec2bin(sML,2);
    MLb2 = reshape(MLb1,[],8);
    MLd1 = bin2dec(MLb2);%二进制转十进制
    MLd2 = reshape(MLd1,size(img))/255;
    
    figure(4)
    subplot(length(SNRdB),1,i);
    imshow(im2uint8(MLd2));
    title(['SNRdB=',num2str(SNRdB(i)),'dB']);
end

% 绘制结果图像
figure(5);
semilogy(SNRdB, berZF, 'b-o',...
    SNRdB, berMMSE, 'r-s',...
    SNRdB, berSIC, 'g-d',...
    SNRdB, berML, 'm-p', 'LineWidth', 2);
grid on;
xlabel('SNR(dB)');
ylabel('BER');
legend('ZF', 'MMSE', 'SIC', 'ML');

figure(6);
semilogy(SNRdB, timeZF, 'b-o',...
    SNRdB, timeMMSE, 'r-s',...
    SNRdB, timeSIC, 'g-d',...
    SNRdB, timeML, 'm-p','LineWidth', 2);
grid on;
xlabel('SNR(dB)');
ylabel('Time (s)');
legend('ZF', 'MMSE', 'SIC', 'ML');