OSIC(排序的连续干扰消除)信号检测

1️⃣ OSIC信号检测介绍

前面的文章提到线性检测方法，例如ZF检测和MMSE检测的性能比非线性检测（LS检测）方法要差，然而线性检测方法的硬件实现复杂度更低。

通过排序的连续干扰消除（Ordered Successive Interference Cancellation，OSIC）方法可以改善线性检测方法的性能，而不会显著提高复杂度。

---

2️⃣ 一些铺垫

先看一下一些符号的定义：

信道矩阵
信道矩阵 $\mathbf{H}$ 是一个 $N_{\mathrm{Rx}}\times N_{\mathrm{Tx}}$ 的矩阵，其中每个元素 $h_{ji}$ 表示从第 $i$ 根发射天线到第 $j$ 根接收天线的信道增益：

$$\mathbf{H}=\left[\begin{array}{cccc}{h}_{11}& h_{12}& \cdots & h_{1N_{\mathrm{Tx}}}\\ h_{21}& h_{22}& \cdots & h_{2N_{\mathrm{Tx}}}\\ ⋮& ⋮& \ddots & ⋮\\ h_{N_{\mathrm{Rx}}1}& h_{N_{\mathrm{Rx}}2}& \cdots & h_{N_{\mathrm{Rx}}{N}_{\mathrm{Tx}}}\end{array}\right]$$

信道模型的矩阵形式：

给定发射信号 $\mathbf{x}$ ，接收信号 $\mathbf{y}$ ，信道矩阵 $\mathbf{H}$ ，以及噪声 $\mathbf{z}$ ，我们可以将整个 MIMO 系统的信号模型表示为：

$$\mathbf{y}=\mathbf{H}\mathbf{x}+\mathbf{z}$$

• $\mathbf{H}$ 是 $N_{\mathrm{Rx}}\times N_{\mathrm{Tx}}$ 的信道矩阵。
• $\mathbf{x}$ 是 $N_{\mathrm{Tx}}\times 1$ 的发射信号向量。
• $\mathbf{y}$ 是 $N_{\mathrm{Rx}}\times 1$ 的接收信号向量。
• $z$ 是 $N_{\mathrm{Rx}}\times 1$ 的噪声向量，表示每个接收天线的噪声。

假设 $\mathbf{y}$ 的每个元素 $y_j$ 是第 $j$ 根接收天线接收到的信号，我们可以将它展开为：

从每个接收天线的角度来看，接收信号$\mathbf{y}$可以分解为：
$$y_j=h_{j1}{x}_1+h_{j2}{x}_2+\cdots +h_{j{N}_{\text{Tx }}}{x}_{N_{\text{Tx }}}+z_j$$

$y_j$ 表示的是，每根接收天线上接收到的信号由多个发射信号加权（权重由信道增益 $h_{ji}$ 给出）并加上噪声 $z_j$ 得到。

从发射信号的角度来看，接收信号$\mathbf{y}$可以分解为：
$$\mathbf{y}={\mathbf{h}}_1{x}_1+{\mathbf{h}}_2{x}_2+\cdots +{\mathbf{h}}_{N_{\mathrm{Tx}}}{x}_{N_{\mathrm{Tx}}}+\mathbf{z}=\left[\begin{array}{c}{h}_{11}\\ h_{21}\\ ⋮\\ h_{N_{\mathrm{Rx}}1}\end{array}\right]x_1+\left[\begin{array}{c}{h}_{12}\\ h_{22}\\ ⋮\\ h_{N_{\mathrm{Rx}}2}\end{array}\right]x_2+\cdots +\left[\begin{array}{c}{h}_{1N_{\mathrm{Tx}}}\\ h_{2N_{\mathrm{Tx}}}\\ ⋮\\ h_{N_{\mathrm{Rx}}{N}_{\mathrm{Tx}}}\end{array}\right]x_{N_{\mathrm{Tx}}}+\mathbf{z}$$
举个例子，$\left[\begin{array}{c}{h}_{11}\\ h_{21}\\ ⋮\\ h_{N_{\mathrm{Rx}}1}\end{array}\right]x_1$表示的是第1根发射天线发送的信号 $x_1$在各个接收天线上的接收情况

---

3️⃣ OSIC原理介绍

1.数据流和符号的定义

• 令 $x_{(i)}$ 表示第 $i$ 个需要检测的符号，这些符号并不一定来自第 $i$ 根发射天线。因为在 OSIC 中，符号的检测是按照顺序设定的【后面介绍】，所以 $x_{(i)}$ 可能是任意数据流的符号。
• ${\hat{x}}_{(i)}$ 表示通过信号检测算法得到的第 $i$ 个符号的估计值。

2.采用ZF或MMSE进行信号检测

• 假设在 OSIC 的第一阶段，我们使用MMSE对符号 $x_{(1)}$ 进行检测，得到它的估计值 ${\hat{x}}_{(1)}$ 。
  $${\hat{x}}_{(1)}={\mathbf{w}}_{(1)}\mathbf{y}$$
  其中 ，MMSE加权矩阵大小为$N_{\mathrm{Tx}}\times N_{\mathrm{Rx}}$。$\mathbf{y}$的大小为$N_{\mathrm{Rx}}\times 1$。${\mathbf{w}}_{(1)}$表示MMSE加权矩阵的第一行

• 然后，从接收信号 $\mathbf{y}$ 中減去 ${\mathbf{h}}_{(1)}{\hat{\mathbf{x}}}_{(1)}$ ，这相当于从接收到的信号中删除第一个数据流的贡献，公式表示为：
  $${\tilde{\mathbf{y}}}_{(1)}=\mathbf{y}-{\mathbf{h}}_{(1)}{\hat{x}}_{(1)}①$$

  • 其中，$\mathbf{y}$表示接收信号，大小$N_{\mathrm{Rx}}\times 1$；
  • ${\mathbf{h}}_{(1)}$是一个列向量，大小$N_{\mathrm{Rx}}\times 1$，第一个数据流不知道是哪个天线上发出的，如果是第一个天线，那就是$\left[\begin{array}{c}{h}_{11}\\ h_{21}\\ ⋮\\ h_{N_{\mathrm{Rx}}1}\end{array}\right]$
  • ${\hat{x}}_{(1)}$表示第一个数据流的估计值，大小1×1；
  • ${\mathbf{h}}_{(1)}{\hat{x}}_{(1)}$表示的是第一个数据流【这里不是真实值，采用的估计值${\hat{x}}_{(1)}$】在各个接受天线上的接收情况。
  • ${\tilde{\mathbf{y}}}_{(1)}$表示接收到的信号中删除第一个数据流的贡献后的信号

• 将下面的公式代入①：
  $$\mathbf{y}={\mathbf{h}}_{(1)}{x}_{(1)}+{\mathbf{h}}_{(2)}{x}_{(2)}+\cdots +{\mathbf{h}}_{\left(N_{\mathrm{TX}}\right)}{x}_{\left(N_{\mathrm{TX}}\right)}+\mathbf{z}$$
  $x_{(1)}$表示第一个数据流的真实值，${\mathbf{h}}_{(1)}$是列向量
  可以得到：
  $${\tilde{\mathbf{y}}}_{(1)}={\mathbf{h}}_{(1)}\left(x_{(1)}-{\hat{x}}_{(1)}\right)+{\mathbf{h}}_{(2)}{x}_{(2)}+\cdots +{\mathbf{h}}_{\left(N_{\mathrm{TX}}\right)}{x}_{\left(N_{\mathrm{TX}}\right)}+\mathbf{z}$$

3.两种情况分析

• 情况 1：$x_{(1)}={\hat{x}}_{(1)}$
  符号 $x_{(1)}$ 已经成功地估计，即 $x_{(1)}={\hat{x}}_{(1)}$ ，那么剩余信号中就不会有第一个数据流的干扰。此时，剩余信号中的干扰主要来自其他数据流，并且这些干扰在后续阶段可以得到有效消除。
• 情况 2：$x_{(1)}\ne {\hat{x}}_{(1)}$
  如果在第一阶段符号估计时产生了误差，即 $x_{(1)}\ne {\hat{x}}_{(1)}$ ，那么后续阶段的检测就会面临差错传播的问题

4.总结
在 OSIC 方法中，每个阶段都检测一个数据流并将其从接收信号中减去，以减小对后续信号的干扰。然而，这种方法的 性能依赖于前一个符号的估计是否准确。如果前一个符号的估计存在误差，后续的符号检测将受到影响，导致差错传播。

---

4️⃣ 检测排序方法

由于前一阶段的错误判决会引起误差传播，因此检测顺序会明显影响 OSIC 检测的整体性能。接下来我们将描述不同的检测排序方法：

4.1 基于SINR的排序

步骤一：
对于每个符号，首先计算 SINR，假设使用线性 MMSE (Minimum Mean Squared Error) 检测，每个符号的 SINR 计算公式为：
$${\mathrm{SINR}}_i=\frac{E_x{\left|{\mathbf{w}}_{i,\mathrm{MMSE}}{\mathbf{h}}_i\right|}^2}{E_x{\sum }_{l\ne i}{\left|{\mathbf{w}}_{i,\mathrm{MMSE}}{\mathbf{h}}_i\right|}^2+{\sigma }_z^2{\Vert {\mathbf{w}}_{i,\mathrm{MMSE}}\Vert }^2}②$$

其中：

• $E_x$是发送符号的能量。
• ${\mathbf{w}}_{i,\mathrm{MMSE}}$是MMSE加权矩阵中的第$i$行。加权矩阵大小为$N_{\mathrm{Tx}}\times N_{\mathrm{Rx}}$
• ${\mathbf{h}}_i$是信道矩阵中的第$i$列，表示第i个符号在各个天线上的接收情况

步骤二：
然后选择SINR最大的符号，使用 MMSE检测或其他线性检测方法对其进行检测。假设选择的符号是 $x_{(l)}$，其估计值为 ${\hat{x}}_{(l)}$。然后，将该符号的贡献从接收信号中去除。那么新的接收信号为：
$${\tilde{\mathbf{y}}}_{(l)}=\mathbf{y}-{\mathbf{h}}_{(l)}{\hat{x}}_{(l)}①$$
步骤三：
在去除干扰后，信道矩阵 $\mathbf{H}$ 也需要更新。去掉第 $l$ 个符号后的新的信道矩阵为：
$${\mathbf{H}}^{(l)}=\left[{\mathbf{h}}_1,{\mathbf{h}}_2,\cdots ,{\mathbf{h}}_{l-1},{\mathbf{h}}_{l+1},\cdots ,{\mathbf{h}}_{N_{\mathrm{Tx}}}\right]$$
然后重新计算MMSE加权矩阵：
$${\mathbf{W}}_{\text{MMSE }}={\left({{\mathbf{H}}^{(\mathbf{l})}}^{\mathrm{H}}{\mathbf{H}}^{(\mathbf{l})}+{\sigma }_z^2\mathbf{I}\right)}^{-1}{{\mathbf{H}}^{(\mathbf{l})}}^{\mathrm{H}}$$
步骤四：
将加权矩阵代入到②中重新计算SINR，此时我们需要计算剩余的 $N_{\mathrm{Tx}}-1$ 个符号的 SINR，即：
{ SINR ⁡ i } i = 1 , i ≠ l N T x \left\{\operatorname{SINR}_i\right\}_{i=1, i \neq l}^{N_{\mathrm{Tx}}} {SINRi​}i=1,i=lNTx​​
计算所有剩余符号的 SINR 后，再次选择 SINR 最大的符号

最终，重复步骤二到步骤四，直到所有符号都被检测。

定量分析：
计算总的 SINR 次数为：

$$\sum _{j=1}^{N_{\mathrm{Tx}}}j=N_{\mathrm{Tx}}\left(N_{\mathrm{Tx}}+1\right)/2$$

4.2 基于SNR的排序

ZF 加权矩阵
ZF 方法通过设计加权矩阵，使得接收信号的不同符号之间的干扰被最大程度地消除，从而在接收端独立地接收每个符号。

SINR 的简化公式
对于符号 $x_i$，ZF 加权后的信号功率为 $|{\mathbf{w}}_i{\mathbf{h}}_l{|}^2=1$，即每个符号的信号功率被标准化。此时检测后的 SINR 简化为
$${\mathrm{SNR}}_i=\frac{E_x}{{\sigma }_z^2{\Vert {\mathbf{w}}_i\Vert }^2},i=1,2,\cdots ,N_{\mathrm{Tx}}$$

• $E_x$ 是发送信号的平均能量。
• ${\sigma }_z^2$ 是接收端的噪声功率。
• $|{\mathbf{w}}_i{|}^2$ 是加权矩阵 ${\mathbf{w}}_i$ 的范数的平方，表示了加权矩阵的功率增益

排序过程
和基于SINR的排序类似，但这里使用的是 SNR 来进行排序

定量分析：
和SINR一样，需要计算SNR值的次数为：
$$\sum _{j=1}^{N_{\mathrm{Tx}}}j=\frac{N_{\mathrm{Tx}}\left(N_{\mathrm{Tx}}+1\right)}{2}$$

4.3 基于列范数的排序

信道矩阵列向量的范数
基于SINR和SNR的方法需要计算每个符号的 SINR 或 SNR 来进行排序，在大规模MIMO系统中，计算可能非常庞大。
为了简化这一过程，提出了 基于信道矩阵的列范数 的排序方法。接收信号表示为：
$$\mathbf{y}=\mathbf{H}\mathbf{x}+\mathbf{z}={\mathbf{h}}_1{x}_1+{\mathbf{h}}_2{x}_2+\cdots +{\mathbf{h}}_{N_{\mathrm{Tx}}}{x}_{N_{\mathrm{Tx}}}+\mathbf{z}$$
其中，${\mathbf{h}}_i$ 是信道矩阵的第 $i$ 个列向量，它表示从第 $i$ 根发射天线到接收天线的信道增益。

可以看出，第 $i$ 个发射符号 $x_i$ 的接收信号强度与信道矩阵的第 $i$ 个列向量 ${\mathbf{h}}_i$ 的范数 $|{\mathbf{h}}_i|$ 成正比，列范数越大的信道，其对应的信号的接收强度越大，也就是该信号在接收端的贡献越强。

因此，信道矩阵的列向量范数 $|{\mathbf{h}}_i|$ 可以作为判断符号检测优先级的一种指标

排序过程
基于列范数的排序方法可以通过以下步骤进行：

• 计算列范数：对于信道矩阵 $\mathbf{H}$，计算每个列向量 ${\mathbf{h}}_i$ 的范数 $|{\mathbf{h}}_i|$，即计算每个信道分量的强度。
• 排序：按照列范数的大小对信号进行排序。列范数大的信号对应的符号应该优先检测，因为它们的接收强度较大，检测误差较小。
• 检测顺序：根据排序结果，从列范数最大的信号开始检测，逐个去除已检测的信号对接收信号的干扰。此时，剩下的信号会变得更加独立，从而减少干扰。

定量分析：
假设系统有 $N_{\mathrm{Tx}}$ 个符号（或数据流），则只需计算 $N_{\mathrm{Tx}}$ 个列范数，并进行一次排序。

---

5️⃣ 代码

5.1 创建名为OSIC_plot.m的文件

%%%%%%%%%%%%%%%%%%%实例仿真说明OSIC性能%%%%%%%%%%%%%%%%%%%%
clear all;
close all;

%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%
%信噪比（SNR）范围，从 0 dB 到 20 dB
SNR = 0:1:20;
%10.^(SNR/10)将信噪比转成线性形式
%信噪比是信号功率比上噪声功率，很多情况下噪声功率就是噪声的方差sigma2
%假设信号的功率为 1（单位功率），所以噪声方差就是SNR的倒数
sigma2 = 1./10.^(SNR/10);
%每个 SNR 点要进行 1000 次独立仿真
num_iteration = 1000;
%发送天线4
Ntx = 4;
%接收电线;
Nrx = 4;
%每次仿真中每个发射天线发送的符号数
num_symbol = 1000;
%16QAM调制
modulation = 16;
errors =zeros(length(SNR),3);
run_times = zeros(length(SNR), 3);  % 用来存储运行时间
%%%%%%%%%%%%%%%%%%%%%%%主循环%%%%%%%%%%%%%%%%%%%%%%%
for index_snr = 1:length(SNR)
    for iteration = 1:num_iteration
        %信道矩阵的大小是Nrx*Ntx
        %生成一个Nrx × Ntx的随机信道矩阵H，该矩阵的元素服从复高斯分布，即每个元素的实部和虚部都服从标准正态分布
        H = sqrt(0.5)*(randn(Nrx,Ntx)+1i*randn(Nrx,Ntx));
        %生成一个Ntx × num_symbol的矩阵x，每个元素是一个从0到15的随机整数
        x = randi([0,modulation-1],Ntx,num_symbol);
        %对符号x进行16QAM调制，得到Ntx × num_symbol的复数矩阵x_modulation
        x_modulation = qammod(x,modulation);
        %生成接收噪声N，大小是Nrx × num_symbol，其每个元素服从复高斯分布
        N = sqrt(sigma2(index_snr)/2).*(randn(Nrx,num_symbol)+1i*randn(Nrx,num_symbol));
        %接收信号，大小是Nrx × num_symbol
        y = H*x_modulation+N;
        % 对每种OSIC类型进行处理
        for OSIC_type = 1:3
            % 调用OSIC函数进行符号检测，得到估计符号x_hat
            x_hat = OSIC_detector(y, H, Ntx, sigma2(index_snr), OSIC_type);
            % 对估计符号x_hat进行解调，得到x_hat_demod
            x_hat_demod = qamdemod(x_hat, modulation);
            errors(index_snr, OSIC_type) = errors(index_snr, OSIC_type) + sum(sum(x_hat_demod ~= x));
        end
    end
end
error_rate = errors/(num_symbol*Ntx*num_iteration);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%作图比较%%%%%%%%%%%%%%%%%%
plot(SNR,error_rate(:,1),'r-o',SNR,error_rate(:,2),'b-o',SNR,error_rate(:,3),'g-s');
xlabel('SNR(dB)');
ylabel('BER');
title('OSIC三种类型的误码率比较');
legend('SINR排序','列范数排序','SNR排序');

5.2 创建名为OSIC_detector.m的函数

function [X_hat]=OSIC_detector(y,H,nT,sigma2,OSIC_type)
%y是输出信号，H是信道参数，nT是发射天线数，sigma2噪声分布的方差
%OSIC_type是排序的类型，1表示SINR排序，2表示列范数排序,3表示SNR
if OSIC_type==1  % Post_detection_SINR
  Order=[];  % detection order
  index_array=1:nT; % yet to be detected signal index

  X_hat = zeros(nT,size(y,2));
  for stage = 1:nT
     %MMSE的加权矩阵
     Wmmse=inv(H'*H+sigma2*eye(nT+1-stage))*H';
     WmmseH=Wmmse*H;
     SINR=[];
     %计算每个符号的 SINR
     for i =1:nT-(stage-1)
        tmp= norm(WmmseH(i,[1:i-1 i+1:nT-(stage-1)]))^2 ...
              + sigma2*norm(Wmmse(i,:))^2;
        SINR(i)=abs(WmmseH(i,i))^2/tmp; 
     end
     %index_temp是具有最大SINR的符号的索引
     [~,index_temp] = max(SINR);
     %将具有最大SINR的符号索引index_max加入到Order中，记录符号的检测顺序
     Order = [Order index_array(index_temp)]; 
     %从index_array中移除已经选择的符号，确保下一次选择时不会再选择已经检测的符号
     index_array = index_array([1:index_temp-1 index_temp+1:end]); 
     %通过MMSE信号检测估计的符号
     X_hat(stage,:) = Wmmse(index_temp,:)*y;
     %删除第一个数据流的贡献
     %   H的大小是Nrx*Ntx
     %   H(:,index_max)是一个列向量，大小为Nrx*1，选哪一列，取决于当前检测的是哪个发射天线上的符号
     y_tilde = y - H(:,index_temp)*X_hat(stage); % interference subtraction
     H_tilde = H(:,[1:index_temp-1 index_temp+1:nT-(stage-1)]); % new H
     %更新信道矩阵H，把符号对应的那列H删除掉
     H = H_tilde;
     %更新y
     y = y_tilde;
  end
  X_hat(Order,:) = X_hat;


elseif OSIC_type==2 
  G = inv(H);
  %计算每个符号的列向量范数
  for i=1:nT
     norm_array(i) = norm(H(:,i));
  end 
  %按照信道列范数的大小对Order_temp进行排序，得到排序后的索引 index_sort
  [~,Order_temp] = sort(norm_array); 
  %使用 wrev() 函数将排序后的索引反转，使得索引从最大列范数的信号开始排序
  Order = wrev(Order_temp);
    
  for stage=1:nT
      %估计当前符号
      X_hat(stage,:)=G(Order(stage),:)*y; 
      %删除已检测符号的干扰
      y_tilde = y-H(:,Order(stage))*X_hat(Order(stage)); 
  end


else 
  Order=[];   
  index_array=1:nT; 
  for stage=1:nT
     % 计算ZF加权矩阵
     G = inv(H'*H)*H'; 
     norm_array=[];
     % 计算每个符号的SNR
     for i=1:nT-(stage-1)
        norm_array(i) = norm(G(i,:));
     end
     [~,index_max]=max(norm_array);
     % 将具有最大SNR的符号加入到Order数组中
     Order=[Order index_array(index_max)];
     % 从index_array中移除已经选择的符号
     index_array = index_array([1:index_max-1 index_max+1:end]); 
     % 通过ZF加权矩阵估计当前符号
     X_hat(stage,:) = G(index_max,:)*y;
     % 删除已检测符号的干扰
     y_tilde = y-H(:,index_max)*X_hat(stage,:);
     % 更新信道矩阵H，删除已检测符号对应的列
     H_tilde = H(:,[1:index_max-1 index_max+1:nT-(stage-1)]);
     H = H_tilde;   
     y = y_tilde;
  end
  % 根据Order数组中的顺序重新排列X_hat
  X_hat(Order,:) = X_hat;
end

结果分析：

• SNR与误码率（BER）关系：随着SNR的增加，误码率会逐渐下降，因为信号噪声比提高，传输质量更好
• 这里关于三种方法的对比，性能可能有问题，后续需要用到再回来看，理论上SINR应该最好的

---

6️⃣ 总结

• 排序连续干扰消除（OSIC）就是采用一定的顺序，逐阶段检测信号。排序方法包括基于SINR、SNR和列向量范数