一、简介
低密度校验码（LDPC码）是一种前向纠错码，LDPC码最早在20世纪60年代由Gallager在他的博士论文中提出，但限于当时的技术条件，缺乏可行的译码算法，此后的35年间基本上被人们忽略，其间由Tanner在1981年推广了LDPC码并给出了LDPC码的图表示，即后来所称的Tanner图。1993年Berrou等人发现了Turbo码，在此基础上，1995年前后MacKay和Neal等人对LDPC码重新进行了研究，提出了可行的译码算法，从而进一步发现了LDPC码所具有的良好性能，迅速引起强烈反响和极大关注。经过十几年来的研究和发展，研究人员在各方面都取得了突破性的进展，LDPC码的相关技术也日趋成熟，甚至已经开始有了商业化的应用成果，并进入了无线通信等相关领域的标准。
1 LDPC码的特点
LDPC码是一种分组码，其校验矩阵只含有很少量非零元素。正是校验矩阵的这种稀疏性，保证了译码复杂度和最小码距都只随码长呈现线性增加。除了校验矩阵是稀疏矩阵外，码本身与任何其它的分组码并无二致。其实如果现有的分组码可以被稀疏矩阵所表达，那么用于码的迭代译码算法也可以成功的移植到它身上。然而，一般来说，为现有的分组码找到一个稀疏矩阵并不实际。不同的是，码的设计是以构造一个校验矩阵开始的，然后才通过它确定一个生成矩阵进行后续编码。而LDPC的编码就是本文所要讨论的主体内容。
译码方法是LDPC码与经典的分组码之间的最大区别。经典的分组码一般是用ML类的译码算法进行译码的，所以它们一般码长较小，并通过代数设计以减低译码工作的复杂度。但是LDPC码码长较长，并通过其校验矩阵H的图像表达而进行迭代译码，所以它的设计以校验矩阵的特性为核心考虑之一。
1.2 LDPC码的构造
构造二进制LDPC码实际上就是要找到一个稀疏矩阵H作为码的校验矩阵，基本方法是将一个全零矩阵的一小部分元素替换成1，使得替换后的矩阵各行和各列具有所要求的数目的非零元素。如果要使构造出的码可用，则必须满足几个条件，分别是无短环，无低码重码字，码间最小距离要尽可能大。
1.3 Tanner图
LDPC码常常通过图来表示，而Tanner图所表示的其实是LDPC码的校验矩阵。Tanner图包含两类顶点:n个码字比特顶点(称为比特节点)，分别与校验矩阵的各列相对应和m个校验方程顶点(称为校验节点)，分别与校验矩阵的各行对应。校验矩阵的每行代表一个校验方程，每列代表一个码字比特。所以，如果一个码字比特包含在相应的校验方程中，那么就用一条连线将所涉及的比特节点和校验节点连起来，所以Tanner图中的连线数与校验矩阵中的1的个数相同。以下图是矩阵
的Tanner图，其中比特节点用圆形节点表示，校验节点用方形节点表示，加黑线显示的是一个6循环：
Tanner图中的循环是由图中的一群相互连接在一起的顶点所组成的，循环以这群顶点中的一个同时作为起点和终点，且只经过每个顶点一次。循环的长度定义为它所包含的连线的数量，而图形的围长，也可叫做图形的尺寸，定义为图中最小的循环长度。如上图中，图形的尺寸，即围长为6，如加黑线所示。
1.4 LDPC编码
方法一：设H=[A | B]，对H进行高斯消元可得到H=[I| P]，设编码完成的码字为u=[c| s]，其中c为监督位，s为信息位。因为Hu’ = uH’ = 0，可以得到Ic’ + Ps’ = 0 即 Ic’ = Ps’ (在GF(2)上)，从而可求 c’ = P*s’。如果高斯消元过程中进行了列交换，则只需记录列交换，并以相反次序对编码后的码字同样进行列交换即可。解码时先求出u，再进行列交换得到u1=[c| s]，后面部分即是想要的信息。
1.5 LDPC码H矩阵的构造
规则的LDPC码和非规则的LDPC码
如果校验矩阵中各行非零元素的个数相同，各列中非零元素的个数也相同，这样的LDPC码称为规则码，与规则码对照，如果校验矩阵的各行中非零元素的个数不同或各列中非零元素个数不同，此时的LDPC码称为非规则的LDPC码。
1.5.1 非规则的LDPC码
构造方法：码率R=M/N=0.5，构造的H矩阵列重固定，而行重是随机的。下面通过一个大小为6*12，列重为2的H矩阵来说明。因为列重固定为2，所以其列坐标集合为[1,1,2,2,3,3,4,4…….12,12]，而行坐标可以随机构造。
1.5.2 规则的LDPC 码
对于码率R=0.5的规则的LDPC码其行重与列重有关，行重是列重的2倍。以612的H矩阵为例，假设其列重为2，每列1的个数为2，总的1的个数为212=24，H矩阵有6行，行重相同则每行1的个数为24/6=4，即行重为4。假设其列重为3，其行重为6。
构造方法：与构造行重类似，行重列重固定之后，其在矩阵中的位置坐标集合固定，如列重为2，行重为4，列坐标集合为[1，1，2，2，3，3，……]，行坐标集合为[1,1,1,1,2,2,2,2,3,3,3,3,……]，下面考虑的就是如何将行列坐标组合起来而不重复。这里采用的是保持行坐标集合不变，列坐标集合打乱重排。以构造大小5*10，列重为2，行重为4的H矩阵为例说明。
1.6 LDPC译码
Gallager 在描述 LDPC 码的时候，分别提出了硬判决译码算法和软判决译码算法两种。经过不断发展，如今的硬判决算法已在 Gallager 算法基础上进展很多，包含许多种加权比特翻转译码算法及其改进形式。硬判决和软判决各有优劣，可以适用于不同的应用场合。
1.6.1 比特翻转算法（BF）
硬判决译码算法最早是 Gallager 在提出 LDPC 码软判决算法时的一种补充。硬判决译码的基本假设是当校验方程不成立时，说明此时必定有比特位发生了错误，而所有可能发生错误的比特中不满足校验方程个数最多的比特发生错误的概率最大。在每次迭代时均翻转发生错误概率最大的比特并用更新之后的码字重新进行译码。具体步骤如下：
二、源代码
close all;
clc;
clear all
warning off
%-----------------读入"隐藏的图片"---------------------
I=imread('W.bmp');
%-----------------------读入"载体图像"-------------------------
cover_image=imread('lena.bmp');
%------------------------------------------------------------------
I0=rgb2gray(I);%灰度化
cover_image=rgb2gray(cover_image);%灰度化
[wm0,watermarked_image,wm]=ldpc_dct(I0,cover_image);%ldpc_dct嵌入提取
e=wm0-wm;
[m,n]=size(e);
mse=sum((e(:).^2))/(m*n);
psnr=10*log10(255^2/mse);%原始水印与提取水印的峰值信噪比
% disp(['ldpc改进dct提取水印的峰值信噪比psnr=',num2str(psnr)])
figure(1)
subplot(221)
imshow(cover_image);
title('原图');
subplot(222);
imshow(I0);
title('水印图');
title('水印图');
%显示嵌入水印后的图象
subplot(223);
uint8_watermarked_image=uint8(watermarked_image);
imshow(uint8_watermarked_image)
title('ldpc编码译码改进后嵌入水印图')
subplot(224);
imshow(double(wm));
title('ldpc编码译码改进后提取水印图')
%% 剪切攻击
I_jianqie=I0;%剪切图
I_jianqie(20:30,20:40)=256;
[wm_jianqie0,watermarked_image_jianqie,wm_jianqie]=ldpc_dct(I_jianqie,cover_image);%ldpc_dct嵌入提取
e_jianqie=wm_jianqie0-wm_jianqie;
[m,n]=size(e_jianqie);
mse_jianqie=sum((e_jianqie(:).^2))/(m*n);
psnr_jianqie=10*log10(255^2/mse_jianqie);%原始水印与提取水印的峰值信噪比
% disp(['攻击后峰值信噪比psnr=',num2str(psnr_jianqie)])
figure(2)
subplot(221)
imshow(cover_image);
title('原图');
subplot(222);
imshow(I0);
title('水印图');
%显示嵌入水印后的图象
subplot(223);
uint8_watermarked_image_jianqie=uint8(watermarked_image_jianqie);
imshow(uint8_watermarked_image_jianqie)
title('剪切攻击后嵌入水印图')
subplot(224);
imshow(double(wm_jianqie));
title('剪切攻击后提取水印图')
%% 高斯噪声
I_gaosi=imnoise(I0,'gaussian',0,0.01);%高斯加噪
[wm_gaosi0,watermarked_image_gaosi,wm_gaosi]=ldpc_dct(I_gaosi,cover_image);%ldpc_dct嵌入提取
e_gaosi=wm_gaosi0-wm_gaosi;
[m,n]=size(e_gaosi);
mse_gaosi=sum((e_gaosi(:).^2))/(m*n);
psnr_gaosi=10*log10(255^2/mse_gaosi);%原始水印与提取水印的峰值信噪比
% disp(['高斯噪声攻击后峰值信噪比psnr=',num2str(psnr_gaosi)])
figure(3)
subplot(221)
imshow(cover_image);
title('原图');
subplot(222);
imshow(I0);
title('水印图');
%显示嵌入水印后的图象
subplot(223);
uint8_watermarked_image_gaosi=uint8(watermarked_image_gaosi);
imshow(uint8_watermarked_image_gaosi)
title('高斯噪声攻击后嵌入水印图')
subplot(224);
imshow(double(wm_gaosi));
title('高斯噪声攻击后提取水印图')
%% 旋转攻击
%%9.rotate 45 旋转
I_xuanzhuan=imrotate(I0,45,'bilinear','crop');%旋转45度
[wm_xuanzhuan0,watermarked_image_xuanzhuan,wm_xuanzhuan]=ldpc_dct(I_gaosi,cover_image);%ldpc_dct嵌入提取
e_xuanzhuan=wm_xuanzhuan0-wm_xuanzhuan;
[m,n]=size(e_xuanzhuan);
mse_xuanzhuan=sum((e_xuanzhuan(:).^2))/(m*n);
psnr_xuanzhuan=10*log10(255^2/mse_xuanzhuan);%原始水印与提取水印的峰值信噪比
% disp(['旋转攻击后峰值信噪比psnr=',num2str(psnr_xuanzhuan)])
figure(4)
subplot(221)
imshow(cover_image);
title('原图');
subplot(222);
imshow(I0);
title('水印图');
function H = makeLdpc(M, N, method, noCycle, onePerCol)
% Create R = 1/2 low density parity check matrix
%
% M : Number of row
% N : Number of column
% method : Method for distributing non-zero element
% {0} Evencol : For each column, place 1s uniformly at random
% {1} Evenboth: For each column and row, place 1s uniformly at random
% noCyle : Length-4 cycle
% {0} Ignore (do nothing)
% {1} Eliminate
% onePerCol: Number of ones per column
%
% H : Low density parity check matrix
%
%
% Copyright Bagawan S. Nugroho, 2007
% http://bsnugroho.googlepages.com
% Number of ones per row (N/M ratio must be 2)
if N/M ~= 2
fprintf('Code rate must be 1/2\n');
end
onePerRow = (N/M)*onePerCol;
% fprintf('Creating LDPC matrix...\n');
switch method
% Evencol
case {0}
% Distribute 1s uniformly at random within column
for i = 1:N
onesInCol(:, i) = randperm(M)';
end
% Create non zero elements (1s) index
r = reshape(onesInCol(1:onePerCol, :), N*onePerCol, 1);
tmp = repmat([1:N], onePerCol, 1);
c = reshape(tmp, N*onePerCol, 1);
% Create sparse matrix H
H = full(sparse(r, c, 1, M, N));
% Evenboth
case {1}
% Distribute 1s uniformly at random within column
for i = 1:N
onesInCol(:, i) = randperm(M)';
end
% Create non zero elements (1s) index
r = reshape(onesInCol(1:onePerCol, :), N*onePerCol, 1);
tmp = repmat([1:N], onePerCol, 1);
c = reshape(tmp, N*onePerCol, 1);
% Make the number of 1s between rows as uniform as possible
% Order row index
[r, ix] = sort(r);
% Order column index based on row index
for i = 1:N*onePerCol
cSort(i, :) = c(ix(i));
end
% Create new row index with uniform weight
tmp = repmat([1:M], onePerRow, 1);
r = reshape(tmp, N*onePerCol, 1);
% Create sparse matrix H
% Remove any duplicate non zero elements index using logical AND
S = and(sparse(r, cSort, 1, M, N), ones(M, N));
H = full(S);
end % switch
% Check rows that have no 1 or only have one 1
for i = 1:M
n = randperm(N);
% Add two 1s if row has no 1
if length(find(r == i)) == 0
H(i, n(1)) = 1;
H(i, n(2)) = 1;
% Add one 1 if row has only one 1
elseif length(find(r == i)) == 1
H(i, n(1)) = 1;
end
end % for i
% If desired, eliminate any length-4 cycle
if noCycle == 1
for i = 1:M
% Look for pair of row - column
for j = (i + 1):M
w = and(H(i, :), H(j, :));
c1 = find(w);
lc = length(c1);
if lc > 1
% If found, flip one 1 to 0 in the row with less number of 1s
if length(find(H(i, :))) < length(find(H(j, :)))
% Repeat the process until only one column left
for cc = 1:lc - 1
H(j, c1(cc)) = 0;
end
三、运行结果
四、备注
版本：2014a
