本次对参考文献[1]中介绍的 V-BLAST，QRD 和 SQRD 算法进行了实现，并对参考文献 [1] 中的实验结果使用Matlab进行了仿真复现；同时根据参考文献 [2] 中给出的基于 MMSE 的线性检测，BLAST 检测，无排序的 QR 检测，有排序的 QR 检测这四种不同的信号检测方法原理使用C++进行了算法实现，并且与参考文献[1]中的算法进行了性能对比。
贝尔实验室分层时空码(BLAST)
Bell Labs Layered Space-Time，即贝尔实验室分层时空，是贝尔实验室提出的一种能够达到极高数据传输速率的无线通信技术。这种技术在一个带限无线信道内，通过多天线技术充分利用空间复用，使得信道容量随着天线数的增加呈线性增长，有着极高的频谱效率。
分层时空(LST)码是一种特殊的时空码，具有易于解码的优点。贝尔实验室最初提出了D-BLAST (对角贝尔实验室分层时空)体系结构，该体系使用对角分层编码结构，其中代码块在时空中分散在对角线上，由此实现了所有层相同的“平均”信道，并且降低了深度衰减的可能性。但是由于代码块的对角线排列，D-BLAST是不能实时实现的。而后提出的V-BLAST(垂直BLAST)，将每个层与特定的发送天线相关联，从而实现了更为容易的检测和解码过程。[3]
多输入多输出系统(MIMO)
多输入多输出系统（Multi-input Multi-output ; MIMO）是一种用来描述多天线无线通信系统的抽象数学模型，能利用发射端的多个天线各自独立发送信号，同时在接收端用多个天线接收并恢复原信息。该技术最早是由马可尼于1908年提出的，他利用多天线来抑制信道衰落（fading）。根据收发两端天线数量，相对于普通的单输入单输出系统（Single-Input Single-Output，SISO），MIMO此类多天线技术尚包含早期所谓的“智能天线”，亦即单输入多输出系统（Single-Input Multi-Output，SIMO）和多输入单输出系统（Multiple-Input Single-Output，MISO）
由于MIMO可以在不需要增加带宽或总发送功率耗损（transmit power expenditure）的情况下大幅地增加系统的资料吞吐量（throughput）及发送距离，使得此技术于近几年受到许多瞩目。MIMO的核心概念为利用多根发射天线与多根接收天线所提供之空间自由度来有效提升无线通信系统之频谱效率，以提升传输速率并改善通信质量。
下图2.1所示为一个有$n_{T}$个发射端天线和$n_{R}$个接收端天线的MIMO系统。
设有发射端信号$\bm c=[c_{1},c_{2},...,c_{n_{T}}]^{T}$，接收端信号$\bm x=[x_{1},x_{2},...,x_{n_{R}}]^{T}$，接受端信号可由下式给出：
$${x}={H}\cdot {c}+{\bm\nu}$$
其中$\bm\nu=[\nu_{1},\nu_{2},...,\nu_{n_{R}}]^{T}$为接收天线处的复高斯白噪音，$H$为$n_{R}\times n_{T}$的矩阵，作为信道增益矩阵，其中$h_{j,i}$为第$i$根发射天线到第$j$根接收天线的IID复衰落增益。
$$
{\bm H} = \Bigg[ \begin{array}{ccc}
h_{1,1} & \dots & h_{1,n_{T}} \
\vdots & \ddots & \vdots \
h_{n_{R},1} & \dots & h_{n_{R},n_{T}}
\end{array} \Bigg]
$$
最小均方误差算法(MMSE)
MMSE探测是以最小化均方误差MSE(Mean Square Error)为标准的信号检测优化方法。[4]是对于无法观察的参数$\theta$ 的一个估计函数$T$；其定义为:
$$\operatorname {MSE}(T)=\operatorname {E}((T-\theta )^{2})$$
即，它是“误差”的平方的期望值。误差就是估计值与被估计量的差。均方误差满足等式
$$\operatorname {MSE}(T)=\operatorname {var}(T)+(\operatorname {bias}(T))^{2}$$
其中
$$\operatorname {bias}(T)=\operatorname {E}(T)-\theta$$
也就是说，偏差$\operatorname {bias}(T)$是估计函数的期望值与那个无法观察的参数的差。
通过最小化实际传输符号和线性探测器输出之间的均方误差，得到滤波矩阵
$$
\mathbf{G}{\mathrm{MMSE}}=\left(\mathbf{H}^{H} \mathbf{H}+\sigma{n}^{2} \mathbf{I}{n{T}}\right)^{-1} \mathbf{H}^{H}
$$
由此得到的滤波器输出
$$
\tilde{\mathrm{s}}{\mathrm{MMSE}}=\mathbf{G}{\mathrm{MMSE}} \mathbf{x}=\left(\mathbf{H}^{H} \mathbf{H}+\sigma_{n}^{2} \mathbf{I}{n{T}}\right)^{-1} \mathbf{H}^{H} \mathbf{x}
$$
不同层的估计误差对应于误差协方差矩阵的对角线元素，得到
$$
\begin{aligned} \Phi_{\mathrm{MMSE}} &=\mathrm{E}\left(\left(\tilde{\mathbf{s}}{\mathrm{MMSE}}-\mathbf{s}\right)\left(\tilde{\mathbf{s}}{\mathrm{MMSE}}-\mathbf{s}\right)^{H}\right) \ &=\sigma_{n}^{2}\left(\mathbf{H}^{H} \mathbf{H}+\sigma_{n}^{2} \mathbf{I}{n{T}}\right)^{-1} \end{aligned}
$$
在QR分解中，为了使得MMSE可以应用于一般的豪斯霍尔德QR分解，可以通过定义$(n_{T}+n_{R})\times n_{T}$扩展信道矩阵$\underline{\mathbf{H}}$和$(n_{T}+n_{R})\times 1$扩展接收向量$\underline{\mathbf{x}}$
$$
\underline{\mathbf{H}}=\left[\begin{array}{c}{\mathbf{H}} \ {\sigma_{n} \mathbf{I}{n{T}}}\end{array}\right] \quad \text { and } \quad \underline{\mathbf{x}}=\left[\begin{array}{c}{\mathbf{x}} \ {\mathbf{0}{n{T}, 1}}\end{array}\right]
$$
MMSE滤波器的输出可以改写为
$$
\tilde{\mathrm{s}}_{\mathrm{MMSE}}=\left(\underline{\mathbf{H}}^{H} \underline{\mathbf{H}}\right)^{-1} \underline{\mathbf{H}}^{H} \underline{\mathbf{x}}=\underline{\mathbf{H}}^{+} \underline{\mathbf{x}}
$$
则误差协方差矩阵变为
$$
\Phi_{\mathrm{MMSE}}=\sigma_{n}^{2}\left(\underline{\mathbf{H}}^{H} \underline{\mathbf{H}}\right)^{-1}=\sigma_{n}^{2} \underline{\mathbf{H}}^{+} \underline{\mathbf{H}}^{+H}
$$
在上述过程中信道矩阵$\mathbf{H}$被替换为$\underline{\mathbf{H}}$，MMSE准则由此引入SQRD算法中。
算法实现
矩阵类
算法的实现使用C++，在算法设计中，采用了封装好的矩阵类，可以进行矩阵加减乘，转置，求逆的操作，为方便阅读，先给出矩阵类中的定义：
class Complex
{
public:
Complex() //构造函数
{
for (int i = 0; i <= 20; i++)
for (int j = 0; j <= 20; j++)
a[i][j] = 0;
n = 0; m = 0;
}
friend Complex operator +(Complex& c1, Complex& c2); //矩阵加法
friend Complex operator -(Complex& c1, Complex& c2); //矩阵减法
friend Complex operator *(Complex& c1, Complex& c2); //矩阵乘法
friend Complex operator !(Complex& c1); //矩阵转置
friend ostream& operator <<(ostream&, Complex&); //矩阵输出
friend istream& operator >>(istream&, Complex&); //矩阵输入
Complex inverse(Complex& c);//矩阵求逆
double a[21][21]; //定义矩阵
int n, m; //标识矩阵的大小,由用户输入,m为矩阵行大小，n为矩阵列大小
};
特别的，给出矩阵转置算法的求解方式，方案为根据$A$和$B$互逆，则$AB=BA=I$，通过初等行变换，可以得到$A$的逆矩阵$B$，代码如下：
Complex inverse(Complex mat)
{
Complex x;
x.m = x.n = mat.n;
for (int i = 1; i <= mat.n; i++) //矩阵扩展
mat.a[i][mat.n + i] = 1;
for (int i = 1; i <= mat.n; i++) //进行初等行变换
{
int max_row = i;
double max = fabs(mat.a[i][i]);
for (int k = i; k <= mat.n; k++)
if (fabs(mat.a[k][i]) > max)
{
max = fabs(mat.a[k][i]);
max_row = k;
}
for (int j = 1; j <= 2 * mat.n; j++)
std::swap(mat.a[i][j],mat.a[max_row][j]);
if (mat.a[i][i] == 0)
return x;
for (int k = 1; k <= mat.n; k++)
for (int j = 2 * mat.n; j >= i; j--)
if (k != i)
mat.a[k][j] = mat.a[k][j] - mat.a[i][j] * mat.a[k][i] / mat.a[i][i];
}
if (mat.a[mat.n][mat.n] == 0)
return x;
for (int i = 1; i <= mat.n; i++) //求解得到逆矩阵
for (int j = 1; j <= mat.n; j++)
x.a[i][j] = mat.a[i][j + mat.n] / mat.a[i][i];
return x;
}
高斯噪声生成算法
高斯噪声是一种具有正态分布（也称作高斯分布）概率密度函数的噪声。换句话说，高斯噪声的值遵循高斯分布或者它在各个频率分量上的能量具有高斯分布。它被极其普遍地应用为用以产生加成性高斯白噪声（AWGN）的迭代白噪声。
高斯噪声满足正态分布，其概率密度函数为
$$
{\displaystyle f(x)={\frac {1}{\sigma {\sqrt {2\pi }}}}e^{-{\frac {\left(x-\mu \right)^{2}}{2\sigma ^{2}}}}!}
$$
生成高斯噪声的算法代码如下：
double gaussrand()
{
static double V1, V2, S;
static int phase = 0;
double X;
if (phase == 0)
{
do
{
double U1 = (double)rand() / RAND_MAX;
double U2 = (double)rand() / RAND_MAX;
V1 = 2 * U1 - 1;
V2 = 2 * U2 - 1;
S = V1 * V1 + V2 * V2;
} while (S >= 1 || S == 0);
X = V1 * sqrt(-2 * log(S) / S);
}
else
X = V2 * sqrt(-2 * log(S) / S);
phase = 1 - phase;
return X;
}
传统V-BLAST算法
接收信号$x$是发射信号$c$的线性组合, 由式子$x=H\cdot{c}+v$给出。接收信号是发射信号的线性组合。由式子$x=H\cdot{c}+v$可以得到，在接收端恢复信号$n_T$的最佳方法是最大似然检测，但由于其复杂性过大，这种方法并不可行。V-BLAST算法通过用零强迫(ZF)零向量线性加权接收信号向量来消除干扰[5]。在每个检测步骤中，除一个信号外的所有信号都被视为干扰信号。通过将消零向量应用于干扰抵消，消除这些信号的影响，检测目标信号
选出2-范数最小的一行$i$
$${w}_i^T = {g}^{(i)}$$
接收信号{x}被${w}_i^T$线性加权，得到
$$y_i = {w}_i^T\cdot {x} ={g}^{(i)}\cdot ({H}\cdot {c}+ {\bm\nu})= c_i + \hat{\nu}_i$$
以此作为第$i$个子信号传输的最终统计量，对其进行调制处理，同时把检测出的信号$\hat{c}_{i}$引起的干扰从接收信号{x}中消去
$${x}{i+1}={x}{i}-{h}{i} \cdot \hat{c}{i}$$
其中${h}_{i}$满足
$$\mathbf{w}{i}^{T} \cdot \mathbf{h}{l}=\left ( \begin{array}{ll}{1} & {l=i} \ {0} & {l \neq i}\end{array} \right.$$
上述过程迭代进行，直到所有的发射信号分量被检测出。
实现该算法的代码如下
void V_BLHST(Complex H, Complex x, int nr, int nt, double* c)
{
Complex HT;
Complex G;
for (int i = 1; i <= nt; i++)
s[i] = i;
for (int i = 1; i <= nt; i++)
{
HT = !H;
G = (HT * H);
G = inverse(G); //矩阵转置
G = G * HT; //对于每一次，计算出矩阵G
double norm2 = 0, min = 0;
for (int j = 1; j <= G.n; j++) //计算范数
min += G.a[1][j] * G.a[1][j];
int m = 1;
for (int j = 1; j <= nt - i + 1; j++)
{
norm2 = 0;
for (int k = 1; k <= G.n; k++)
norm2 += G.a[j][k] * G.a[j][k];
if (norm2 < min) //取最大范数
{
min = norm2;
m = j;
}
}
double y = 0;
for (int j = 1; j <= G.n; j++)
y += G.a[m][j] * x.a[j][1]; //计算取值y
if (y < 0)//使用y值对接受矩阵进行修正
c[m] = -1;
if (y > 0)
c[m] = 1;
for (int j = 1; j <= G.n; j++) //计算迭代的x
x.a[j][1] = x.a[j][1] - c[m] * H.a[j][m];
for (int j = 1; j <= H.m; j++) //将已经计算好的列换到最后
{
H.a[j][m] = 0;
std::swap(H.a[j][m], H.a[j][nt - i + 1]);
}
H.n--;
std::swap(c[m], c[nt - i + 1]);
std::swap(s[m], s[nt - i + 1]);
}
for (int i = 1; i <= nt; i++) //恢复原始顺序
{
for (int j = 1; j <= nt - 1; j++)
{
if (s[j] > s[j + 1])
{
std::swap(c[j], c[j + 1]);
std::swap(s[j], s[j + 1]);
}
}
}
}
未排序QRD算法
QR分解算法在V-BLAST算法基础上更进一步优化。[6]信道增益矩阵{H}可进行QR分解${H}={Q}\cdot {R}$
将其左乘{Q}的埃尔米特矩阵，得到
$$\mathbf{y}=\mathbf{Q}^{H} \cdot \mathbf{x}=\mathbf{R} \cdot \mathbf{c}+\boldsymbol{\eta}$$
由于{Q}是幺正矩阵，所以噪声$\bm\eta={Q}^H\cdot \bm\nu$的统计特性不变。y的各个分量为
$$y_k = r_{k,k}\cdot c_k + \eta_k + d_k$$
其中干扰项为
$$d_k = \sum\limits_{i=k+1}^{n_T} r_{k,i}\cdot c_i$$
因为{R}是一个$n_{T}\times n_{T}$的上三角矩阵，所以该干扰项与上层信号$c_1,c_2,\cdots ,c_{k-1}$无关，故下层信号$c_{n_T}$表示如下
$$y_{n_T} = r_{n_T,n_T}\cdot c_{n_T} + \eta_{n_T}$$
且其独立于其他发射信号，因此，估计值可表示为
$$\hat{c}{n_T} = \mathcal{Q}\Big[\frac{y{n_T}}{r_{n_T,n_T}}\Big]$$
对于$n_{T}-1$层信号的检测，可以表示为
$$y_{n_{T}-1}=r_{n_{T}-1, n_{T}-1} \cdot C_{n_{T}-1}+r_{n_{T}-1, n_{T}} \cdot c_{n_{T}}+\eta_{n_{T}-1}$$
在认为$\hat{c}{n{T}}=c_{n_{T}}$ 的情况下，得到了用于估计$c_{n_{T-1}}$的无干扰决定统计量。检测层$k=n_{T-1},...1$同理。根据先前得到的估计值，可以计算干扰项$\hat{d_k}$并在接收信号$y_k$中消去，计算式如下
$$z_k = y_k -\hat{d_k} = r_{k,k}\cdot c_k +\eta_k$$
假定$\hat{d_k}=d_k$，那么$z_k$可视为零干扰，$c_k$可以被检测
$$\hat{c_k}=\mathcal{Q}\Big[\frac{z_{k}}{r_{k,k}}\Big]$$
实现该算法的代码如下
void QRD(Complex H, Complex x, int nr, int nt, double* c)
{
Complex R, Q;
Q.n = H.m;
Q.m = H.m;
for (int i = 1; i <= H.m; i++)
for (int j = 1; j <= H.m; j++)
{
Q.a[i][j] = R.a[i][j] = 0;
}
R.m = H.m;
R.n = H.n;
double a[20], b[20];
//QR分解
for (int j = 1; j <= H.m; j++)
{
for (int i = 1; i <= H.m; i++)
a[i] = b[i] = H.a[i][j];
for (int i = 1; i < j; i++)
{
R.a[i][j] = 0;
for (int m = 1; m <= H.m; m++)
R.a[i][j] += a[m] * Q.a[m][i];
for (int m = 1; m <= H.m; m++)
b[m] -= R.a[i][j] * Q.a[m][i];
}
double f2 = 0;
for (int i = 1; i <= H.m; i++)
f2 += b[i] * b[i];
f2 = sqrt(f2);
R.a[j][j] = f2;
if (f2 == 0) break;
for (int i = 1; i <= H.m; i++)
Q.a[i][j] = b[i] / f2;
}
double y[20], z[20];
for (int j = 1; j <= Q.m; j++) //计算参数y
for (int i = 1; i <= Q.m; i++)
y[j] += Q.a[i][j] * x.a[i][1];
for (int k = Q.m; k >= 1; k--)
{
double d = 0;
for (int i = k + 1; i <= Q.m; i++)
d += R.a[k][i] * c[i];
z[k] = y[k] - d;
//量化
if (z[k] / R.a[k][k] < 0)
c[k] = -1;
else
c[k] = 1;
}
}
SQRD算法
SQRD的性能接近于VBLAST。它基本上是对改进的Gram-Schmidt算法[7]的扩展，在每个正交化步骤中将H的列进行排序。本文提出的算法使用搜索实现最小$SNR_k$的检测序列。因此，其对角线元素的绝对值在三角形矩阵的左上方区域的绝对值很小。因此，由小信噪比引起的检测故障$SNR_k$的误差传播只影响$1，…，k−1$这几个层，有效控制了误差的传播。
QR分解算法的信噪比SNR如下
$$SNR_k=\frac{{E}{|c_k|^2}|r_{k,k}|^2}{{E}{|n_{k_i}|^2}}\sim |r_{k,k}|^2$$
可以看到，SNR与$|r_{k,k}|^2$成正比，意味着$r_{k,k}$越大，信噪比越大，检测中收到的干扰越小，故我们首先找到最大的$r_{k,k}$优先处理，以此降低整个算法的误码率。SQRD算法相较于QRD算法，在加入了对于最小范数的检测后，实现了对于误差传播的有效控制，大大降低了误码率。
void SQRD(Complex H, Complex x, int nr, int nt, double* c)
{
Complex R;
Complex Q;
Q = H;
Q.n = H.m;
Q.m = H.m;
R.m = H.m;
R.n = H.n;
for (int i = 1; i <= nt; i++)
s[i] = i;
for (int j = 1; j <= nt; j++)
{
double norm = 0, min = 0; //寻找最大的范数
int min_col = j;
for (int i = 1; i <= H.m; i++)
min += Q.a[i][j] * Q.a[i][j];
for (int k = j + 1; k <= nt; k++)
{
norm = 0;
for (int i = 1; i <= H.m; i++)
norm += Q.a[i][k] * Q.a[i][k];
if (norm < min)
{
min = norm;
min_col = k;
}
}
for (int i = 1; i <= H.m; i++)
std::swap(Q.a[i][j], Q.a[i][min_col]);
for (int i = 1; i <= nt; i++)
std::swap(R.a[i][j], R.a[i][min_col]);
std::swap(s[j], s[min_col]);
R.a[j][j] = sqrt(min);
for (int i = 1; i <= Q.m; i++)
Q.a[i][j] /= R.a[j][j];
for (int k = j + 1; k <= nt; k++)
{
for (int i = 1; i <= Q.m; i++)
R.a[j][k] += Q.a[i][j] * Q.a[i][k];
for (int i = 1; i <= Q.m; i++)
Q.a[i][k] -= R.a[j][k] * Q.a[i][j];
}
}
double y[20], z[20];
for (int j = 1; j <= Q.m; j++) //计算参数y
for (int i = 1; i <= Q.m; i++)
y[j] += Q.a[i][j] * x.a[i][1];
for (int k = Q.m; k >= 1; k--)
{
double d = 0;
for (int i = k + 1; i <= Q.m; i++)
d += R.a[k][i] * c[i];
z[k] = y[k] - d;
//量化
if (z[k] / R.a[k][k] < 0)
c[k] = -1;
else
c[k] = 1;
}
for (int i = 1; i <= nt; i++) //恢复原始顺序
{
for (int j = 1; j <= nt - 1; j++)
{
if (s[j] > s[j + 1])
{
std::swap(c[j], c[j + 1]);
std::swap(s[j], s[j + 1]);
}
}
}
}
MMSE线性探测算法
MMSE准则调整信道矩阵$\mathbf{H}$，并得到其伪逆矩阵
$$
\mathbf{G}{\mathrm{MMSE}}=\left(\mathbf{H}^{H} \mathbf{H}+\sigma{n}^{2} \mathbf{I}{n{T}}\right)^{-1} \mathbf{H}^{H}
$$
以伪逆矩阵左乘接收信号向量，其中复高斯白噪音矢量${\bm\nu}$被忽略，则检测到的发射信号直接由下表达式得到
$$\hat c=\mathbf{G}{\mathrm{MMSE}}{x} = \mathbf{G}{\mathrm{MMSE}}({H}\cdot {c}+{\bm\nu})={c}+\mathbf{G}_{\mathrm{MMSE}}\cdot {\bm\nu}$$
在上式中，$\sigma$为生成的噪声向量的功率，通过将该参数加入到MMSE准则中，实现对信号检测的优化。
void MMSE(Complex H, Complex x, int nr, int nt, double* c, double sigma)
{
Complex HT;
Complex G;
HT = !H;
G = (HT * H);
Complex I;
I.m = I.n = G.m;
for (int i = 1; i <= I.m; i++)
I.a[i][i] = sigma;
G = G + I;
G = inverse(G);
G = G * HT;
for (int i = 1; i <= nt; i++)
for (int j = 1; j <= nr; j++)
{
c[i] += G.a[i][j] * x.a[j][1];
if (c[i] > 0)
c[i] = 1;
else
c[i] = -1;
}
}
MMSE-BLAST算法
MMSE-BLAST与V-BLAST一样的是，在每个检测步骤中，除一个信号外的所有信号都被视为干扰信号。通过将消零向量应用于干扰抵消，消除这些信号的影响，检测目标信号。区别在于ZF准则更改为MMSE准则，应用MMSE准则得到矩阵$\mathbf{G}$
$$
\mathbf{G}{\mathrm{MMSE}}=\left(\mathbf{H}^{H} \mathbf{H}+\sigma{n}^{2} \mathbf{I}{n{T}}\right)^{-1} \mathbf{H}^{H}
$$
选出2-范数最小的一行$i$
$${w}_i^T = {g}^{(i)}$$
接收信号{x}被${w}_i^T$线性加权，得到
$$y_i = {w}_i^T\cdot {x} ={g}^{(i)}\cdot ({H}\cdot {c}+ {\bm\nu})= c_i + \hat{\nu}_i$$
以此作为第$i$个子信号传输的最终统计量，对其进行调制处理，同时把检测出的信号$\hat{c}_{i}$引起的干扰从接收信号{x}中消去
$${x}{i+1}={x}{i}-{h}{i} \cdot \hat{c}{i}$$
上述过程迭代进行，直到所有的发射信号分量被检测出。
void V_BLHST(Complex H, Complex x, int nr, int nt, double* c, double sigma)
{
Complex HT;
Complex G;
for (int i = 1; i <= nt; i++)
s[i] = i;
for (int i = 1; i <= nt; i++)
{
HT = !H;
G = (HT * H);
Complex I;
I.m = I.n = G.m;
for (int i = 1; i <= I.m; i++)
I.a[i][i] = sigma;
G = G + I;
G = inverse(G);
G = G * HT;
double norm2 = 0, min = 0;
for (int j = 1; j <= G.n; j++) //计算范数
min += G.a[1][j] * G.a[1][j];
int m = 1;
for (int j = 1; j <= nt - i + 1; j++)
{
norm2 = 0;
for (int k = 1; k <= G.n; k++)
norm2 += G.a[j][k] * G.a[j][k];
if (norm2 < min) //取最大范数
{
min = norm2;
m = j;
}
}
double y = 0;
for (int j = 1; j <= G.n; j++)
y += G.a[m][j] * x.a[j][1]; //计算取值y
if (y < 0)//使用y值对接受矩阵进行修正
c[m] = -1;
if (y > 0)
c[m] = 1;
for (int j = 1; j <= G.n; j++) //计算迭代的x
x.a[j][1] = x.a[j][1] - c[m] * H.a[j][m];
for (int j = 1; j <= H.m; j++) //将已经计算好的列换到最后
{
H.a[j][m] = 0;
std::swap(H.a[j][m], H.a[j][nt - i + 1]);
}
H.n--;
std::swap(c[m], c[nt - i + 1]);
std::swap(s[m], s[nt - i + 1]);
}
for (int i = 1; i <= nt; i++) //恢复原始顺序
{
for (int j = 1; j <= nt - 1; j++)
{
if (s[j] > s[j + 1])
{
std::swap(c[j], c[j + 1]);
std::swap(s[j], s[j + 1]);
}
}
}
}
无排序MMSE-QRD算法
在比较V-BLAST算法与ZF-QRD算法的时候提到，V-BLAST算法的时间复杂度过大，不适合计算大量的数据，同样的，对于MMSE准则下的MMSE-BLAST算法，其时间复杂度同样很大，所以与前文中提到的类似，可以使用QRD算法来对信号进行检测。
区别于普通QR分解算法，MMSE准则的QR算法使用拓展的信道矩阵$\mathbf{H}$和扩展的接收信号向量$\mathbf{x}$进行计算。
扩展信道矩阵和接收信号向量为
$$
\underline{\mathbf{H}}=\left[\begin{array}{c}
{\mathbf{H}} \
{\sigma_{n} \mathbf{I}{n{T}}}
\end{array}\right] \quad \text { and } \quad \underline{\mathbf{x}}=\left[\begin{array}{c}
{\mathbf{x}} \
{\mathbf{0}{n{T}, 1}}
\end{array}\right]
$$
对扩展信道增益矩阵{H}进行QR分解
$$\underline{\mathbf{H}}=\left[\begin{array}{c}{\mathbf{H}} \ {\sigma_{n} \mathbf{I}{n{T}}}\end{array}\right]=\underline{\mathbf{Q}} \mathbf{R}=\left[\begin{array}{l}{\mathbf{Q}{1}} \ {\mathbf{Q}{2}}\end{array}\right] \underline{\mathbf{R}}=\left[\begin{array}{c}{\mathbf{Q}{1} \mathbf{R}} \ {\mathbf{Q}{2} \mathbf{R}}\end{array}\right]$$
将左乘$\underline{\mathbf{Q}}$的埃尔米特矩阵，得到
$$\mathbf{y}=\underline{\mathbf{Q}}^{H} \cdot \underline{\mathbf{x}}=\underline{\mathbf{R}} \cdot \mathbf{c}+\boldsymbol{\eta}$$
由于$\underline{\mathbf{Q}}$是幺正矩阵，所以噪声$\bm\eta={Q}^H\cdot \bm\nu$的统计特性不变。y的各个分量为
$$y_k = r_{k,k}\cdot c_k + \eta_k + d_k$$
其中干扰项为
$$d_k = \sum\limits_{i=k+1}^{n_T} r_{k,i}\cdot c_i$$
因为$\underline{\mathbf{R}}$是上三角矩阵，所以该干扰项与上层信号$c_1,c_2,\cdots ,c_{k-1}$无关，故下层信号$c_{n_T}$表示如下
$$y_{n_T} = r_{n_T,n_T}\cdot c_{n_T} + \eta_{n_T}$$
且其独立于其他发射信号，因此，估计值可表示为
$$\hat{c}{n_T} = \mathcal{Q}\Big[\frac{y{n_T}}{r_{n_T,n_T}}\Big]$$
对于$n_{T}-1$层信号的检测，可以表示为
$$y_{n_{T}-1}=r_{n_{T}-1, n_{T}-1} \cdot C_{n_{T}-1}+r_{n_{T}-1, n_{T}} \cdot c_{n_{T}}+\eta_{n_{T}-1}$$
在认为$\hat{c}{n{T}}=c_{n_{T}}$的情况下，得到了用于估计$c_{n_{T-1}}$的无干扰决定统计量。检测层$k=n_{T-1},...1$同理。根据先前得到的估计值，可以计算干扰项$\hat{d_k}$并在接收信号$y_k$中消去，计算式如下
$$z_k = y_k -\hat{d_k} = r_{k,k}\cdot c_k +\eta_k$$
假定$\hat{d_k}=d_k$，那么$z_k$可视为零干扰，$c_k$可以被检测
$$\hat{c_k}=\mathcal{Q}\Big[\frac{z_{k}}{r_{k,k}}\Big]$$
void QRD(Complex H, Complex x, int nr, int nt, double* c, double sigma)
{
int row = H.m;
H.m = H.m + H.n;
for (int i = 1; i <= H.m; i++)
{
for (int j = 1; j <= H.n; j++)
{
if (j == i - row)
H.a[i][j] = sigma;
}
}
Complex R, Q;
Q.n = H.m;
Q.m = H.m;
for (int i = 1; i <= H.m; i++)
for (int j = 1; j <= H.m; j++)
{
Q.a[i][j] = R.a[i][j] = 0;
}
R.m = H.m;
R.n = H.n;
double a[20], b[20];
//QR分解
for (int j = 1; j <= H.m; j++)
{
for (int i = 1; i <= H.m; i++)
a[i] = b[i] = H.a[i][j];
for (int i = 1; i < j; i++)
{
R.a[i][j] = 0;
for (int m = 1; m <= H.m; m++)
R.a[i][j] += a[m] * Q.a[m][i];
for (int m = 1; m <= H.m; m++)
b[m] -= R.a[i][j] * Q.a[m][i];
}
double f2 = 0;
for (int i = 1; i <= H.m; i++)
f2 += b[i] * b[i];
f2 = sqrt(f2);
R.a[j][j] = f2;
if (f2 == 0) break;
for (int i = 1; i <= H.m; i++)
Q.a[i][j] = b[i] / f2;
}
double y[20], z[20];
for (int j = 1; j <= Q.m; j++) //计算参数y
for (int i = 1; i <= Q.m; i++)
y[j] += Q.a[i][j] * x.a[i][1];
for (int k = Q.m; k >= 1; k--)
{
double d = 0;
for (int i = k + 1; i <= Q.m; i++)
d += R.a[k][i] * c[i];
z[k] = y[k] - d;
//量化
if (z[k] / R.a[k][k] < 0)
c[k] = -1;
else
c[k] = 1;
}
}
MMSE-SQRD算法
实际上采用MMSE-QRD算法的时候可以发现，其误码率相比MMSE-BLAST算法更高，我们考虑对MMSE准则下的QRD算法进行同样的优化。
QR分解算法的信噪比SNR如下
$$SNR_k=\frac{{E}{|c_k|^2}|r_{k,k}|^2}{{E}{|n_{k_i}|^2}}\sim |r_{k,k}|^2$$
可以看到，SNR与$|r_{k,k}|^2$成正比，意味着$r_{k,k}$越大，信噪比越大，检测中收到的干扰越小，故我们首先找到最大的$r_{k,k}$优先处理，以此降低整个算法的误码率。
为了得到最佳的检测顺序，第一个$|r_{n_{T}, n_{T}}|$必须在扩展信道矩阵$\underline{\mathbf{H}}$的所有可能排列上最大化，其次是$| r_{n_{T}-1, n_{T}-1}|$，以此类推，后面也同样。然而在标准的QR算法中，$\mathbf{R}$的对角线元素的计算顺序是相反的，这使得检测过程变得复杂，因此我们提出基于MMSE原理的改进排序的QR算法。其基本思想是按照(1,...,$n_T$)的检测顺序使$\left|r_{k}, k\right|$最小化。由于最后探测的层通过传播误差只影响少数层，因此可能具有很小的信号与干扰加噪声比，这增加了第一层较大信号与干扰加噪声比的概率。
将信道矩阵$\mathbf{H}$按照这样的顺序进行处理，并由此得到$\mathbf{R}$，得到QR分解的结果。
void SQRD(Complex H, Complex x, int nr, int nt, double* c, double sigma)
{
int row = H.m;
H.m = H.m + H.n;
for (int i = 1; i <= H.m; i++)
{
for (int j = 1; j <= H.n; j++)
{
if (j == i - row)
H.a[i][j] = sigma;
}
}
Complex R;
Complex Q;
Q = H;
Q.n = H.m;
Q.m = H.m;
R.m = H.m;
R.n = H.n;
for (int i = 1; i <= nt; i++)
s[i] = i;
for (int j = 1; j <= nt; j++) //寻找最小范数
{
double norm = 0, min = 0;
int min_col = j;
for (int i = 1; i <= H.m; i++)
min += Q.a[i][j] * Q.a[i][j];
for (int k = j + 1; k <= nt; k++)
{
norm = 0;
for (int i = 1; i <= H.m; i++)
norm += Q.a[i][k] * Q.a[i][k];
if (norm < min)
{
min = norm;
min_col = k;
}
}
for (int i = 1; i <= H.m; i++) //将最小范数提前
std::swap(Q.a[i][j], Q.a[i][min_col]);
for (int i = 1; i <= nt; i++)
std::swap(R.a[i][j], R.a[i][min_col]);
std::swap(s[j], s[min_col]);
R.a[j][j] = sqrt(min);
for (int i = 1; i <= Q.m; i++)
Q.a[i][j] /= R.a[j][j];
for (int k = j + 1; k <= nt; k++)
{
for (int i = 1; i <= Q.m; i++)
R.a[j][k] += Q.a[i][j] * Q.a[i][k];
for (int i = 1; i <= Q.m; i++)
Q.a[i][k] -= R.a[j][k] * Q.a[i][j];
}
}
double y[20], z[20];
for (int j = 1; j <= Q.m; j++) //计算参数y
for (int i = 1; i <= Q.m; i++)
y[j] += Q.a[i][j] * x.a[i][1];
for (int k = Q.m; k >= 1; k--)
{
double d = 0;
for (int i = k + 1; i <= Q.m; i++)
d += R.a[k][i] * c[i];
z[k] = y[k] - d;
//量化
if (z[k] / R.a[k][k] < 0)
c[k] = -1;
else
c[k] = 1;
}
for (int i = 1; i <= nt; i++) //恢复原始顺序
{
for (int j = 1; j <= nt - 1; j++)
{
if (s[j] > s[j + 1])
{
std::swap(c[j], c[j + 1]);
std::swap(s[j], s[j + 1]);
}
}
}
}
Comments 3 条评论
这是一条私密评论
@function 这都被你发现了，这是我们当年的期中作业
@lzusa 这是一条私密评论