原创 关于LDPC的仿真和编码器的LU实现

<?xml:namespace prefix = o ns = "urn:schemas-microsoft-com:office:office" />

目前我了解的仿真环境主要有3个：C环境、Matlab环境和CCSS环境。其中CCSS环境运行在Linux下，主要用来搭建链路。C和Matlab在windows和Linux上均可以建立环境。

下面我主要说一下我仿真LDPC的心得。开始的时候，我自己做了一套环境，从编码到解码，完全是根据自己的实际需要出发，然而仿真的速度却不高。后来总结了一下，问题出在编码器的实现上。在给定一个原始输入向量（信源）x，编码后的码字为y，我的算法是y=x*G，其中G是生成矩阵。我采用了矩阵相乘的做法，这消耗掉了大量的计算资源。其实如果G是一个稀疏矩阵，可以使用很小的代价来解出y：yG’=x，这是一个矩阵方程。现在问题就归结在G的产生上了。这里就不得不提到一种利用系数矩阵实现LDPC的快速编码方法，这就是利用LU分解。

假设校验矩阵H是N*M的矩阵，可以分解为两部分：左部M*M的A矩阵和右部M*N的B矩阵，在列置换后将A变换为非奇异矩阵，即可逆矩阵。编码后码字x也可以分为M比特的校验比特c和N-M比特的信源比特s。从校验矩阵的性质H*x=0可以得出：[A|B]*[c/s]=0，这里的/符号其实表示上下关系。这样可以得到A*c+B*s=0，进而c=A’*B*s，这里的’符号表示矩阵的逆。通过这种方法可以算出校验比特c了。然而A’不是稀疏矩阵，运算量还是比较大的，为O(M*(N-M))。因为H是稀疏矩阵，那么B就是稀疏矩阵，那么用下面的方法就可以快一点算出c来了：

(1)    计算z=B*s，运算量是O(M)。

(2)    计算c=A’*z，运算量是O(M*M)

其实对于A*c=z，我们可以通过行变换把A化为上三角矩阵U，变为U*c=y。其中行变换过程我们可以将其表示为一个下三角矩阵L，这样y就可以通过求解L*y=z得出。这相当于对A做了一次LU分解。我们可以通过下面的步骤来进行分解（包括A的产生）：

(1)    通过行列变换将A变换为可逆矩阵；

(2)    将U和L初始化为零矩阵；

(3)    设立中间矩阵F，初值为H；

(4)    For i="1" to M做以下事情：

a)         在F中找到i行i列的非零元素，如果找不到就在以后的行或列中寻找；

b)        对F和H进行行列变换，把找到的那个元素放到第i行第i列；

c)        把F中第i列的第1行到第i行的内容复制到U的第i列；

d)        把F中第i列的第i行到最后一行的内容复制到L的第i列；

e)         把F中的第i行加到后面的行

(5)    行列变换后的H中的最后N-M列就是B矩阵

于是我们可以用L，U和B矩阵来解出c：

(1)    计算z=B*s

(2)    解L*y=z得出y

(3)    解U*c=y得出c

值得说明的是行列变换的过程中没有最有效的算法，但是有heuristic的算法，成为最小列或最小积原则。

该功能已经在Toronto大学的Radford M. Neal教授开发的程序包内实现。为了加速仿真，我使用了该程序包，该程序包用链表结构实现了稀疏矩阵，对于求解方程、进行矩阵乘法非常有效。但是该程序包有一些不足：

(1)    校验矩阵有着自己的格式，拿到一个校验矩阵后需要转换成这种格式；

(2)    解码算法只有最原始的算法，且以浮点实现，这对于解码器的开发是远远不够的。

于是我针对这两点作了自己的扩展，并在程序统计BER和FER，使之适应于我的研究。该环境我放在附件中，仅供参考。感谢Radford M. Neal教授的研究和贡献。