clear all
clc
s=[1 0 1 0 0 1 0 0 1]'
rows=length(s); %计算输入序列s的长度记做校验矩阵的行,即96
cols=2*length(s); %行乘2为校验矩阵的列
row_flag(1:rows)=0; %定义一个一行96个0的数组 方便后面使用
parity_check=zeros(rows,cols); %校验矩阵初始化
%%%%%%%%使每列随机产生3个1即列重为3%%%%%%%%%%%%%
bits_per_col=3; %每列3个1
for i=1:cols
a=randperm(rows); %把1~96随机打乱
for j=1:bits_per_col
parity_check(a(j),i)=1; %每列3个1
row_flag(a(j))=row_flag(a(j))+1;
end
end
%计算每行1的最多个数
max_ones_per_row=ceil(cols*bits_per_col/rows); %ceil是向离它最近的大整数圆整192*3/96=6 每行1的最多数是6
%add 1's to rows having no 1(a redundant row) or only one 1(that bit in the codeword becomes
%zero irrespective of the input)
for i=1:rows
if row_flag(i)==0 %如果该行没有1,则随机添加两个1
for k=1:2
j=unidrnd(cols); %产生从1到192所指定的最大数数之间的离散均匀随机数
while parity_check(i,j)==1
j=unidrnd(cols);
end
parity_check(i,j)=1; %在找到的新位置上置1
row_flag(i)=row_flag(i)+1; %行重加1
end
end
if row_flag(i)==1 %如果该行只有1个1,则随机再添加1个1
j=unidrnd(cols);
while parity_check(i,j)==1
j=unidrnd(cols);
end
parity_check(i,j)=1;
row_flag(i)=row_flag(i)+1;
end
end
%尝试在列上分散1的位置,使得每行1的个数均衡(相近或相一致)
for i=1:rows
j=1;
a=randperm(cols);
while row_flag(i)>max_ones_per_row; %如果该行行重大于允许的最大行重,则进行处理
if parity_check(i,a(j))==1 %随机选择某一该行上为1的列来处理,将该列该行上的1分散到其他的行
%随机查找该列上适合放置1(行重小于允许的最大行重,且该位置上为0)的行
newrow=unidrnd(rows);
k=0;
while (row_flag(newrow)>=max_ones_per_row || parity_check(newrow,a(j))==1) && k<rows
newrow=unidrnd(rows);
k=k+1;
end
if parity_check(newrow,a(j))==0
%将待处理行上的1转放到找到的行上
parity_check(newrow,a(j))=1;
row_flag(newrow)=row_flag(newrow)+1;
parity_check(i,a(j))=0;
row_flag(i)=row_flag(i)-1;
end
end
j=j+1;
end
end
%尝试删除4环
for loop=1:10
chkfinish=1;
for r=1:rows
ones_position=find(parity_check(r,:)==1);
ones_count=length(ones_position);
for i=[1:r-1 r+1:rows]
common=0;
for j=1:ones_count
if parity_check(i,ones_position(j))==1
common=common+1 ;
if common==1
thecol=ones_position(j);
end
end
if common==2
chkfinish=0; %如果还存在4环,则不结束循环,还进入下一次循环
common=common-1;
if (round(rand)==0) % 随机决定是保留前面的列还是后面的列
coltoberearranged=thecol; %保留后面的列,交换前面的列
thecol=ones_position(j);
else
coltoberearranged=ones_position(j); %保留前面的列,交换后面的列
end
parity_check(i,coltoberearranged)=3; %make this entry 3 so that we dont use
%of this entry again while getting rid
%of other cylces
newrow=unidrnd(rows);
iteration=0; %尝试5次在待交换的列中随机查找0
while parity_check(newrow,coltoberearranged)~=0 & iteration<5
newrow=unidrnd(rows);
iteration=iteration+1;
end
if iteration>=5 %超过5次后则扩大范围随机查找非1的0或3,直到找到为止
while parity_check(newrow,coltoberearranged)==1
newrow=unidrnd(rows);
end
end
%把该列中找到的0或3置为1
parity_check(newrow,coltoberearranged)=1;
end%if common==2
end%for j=1:ones_count
end%for i=[1:r-1 r+1:rows]
end%for r=1:rows
%如果本次循环已不存在4环,则结束循环,不进入下一次循环
if chkfinish
break
end
end%for loop=1:10
%replace the 3's with 0's
parity_check=parity_check==1;
%Get the Parity Checks
H=parity_check;
%第二部分::求p
dim=size(H); %H是96行192列的大小
rows=dim(1);
cols=dim(2);
rearranged_cols=zeros(1, rows);
%逐行进行高斯消元,前面的rows行 x rows列形成单位矩阵
for k=1:rows
vec = [k:cols];
%查找可交换的列
x = k;
while (x<=cols & H(k,x)==0)
ind = find(H(k+1:rows, x) ~= 0);
if ~isempty(ind)
break
end
x = x + 1;
end
%如果找不到可交换的列则视为非法的H矩阵并退出
if x>cols
error('Invalid H matrix.');
end
%如果不是当前列则进行列交换,同时保存交换记录
if (x~=k)
rearranged_cols(k)=x;
temp=H(:,k);
H(:,k)=H(:,x);
H(:,x)=temp;
end
%高斯消元,使G(k,k)==1
if (H(k,k)==0)
ind = find(H(k+1:rows, k) ~= 0);
ind_major = ind(1);
x = k + ind_major;
H(k, vec) = rem(H(x, vec) + H(k, vec), 2);
end
%高斯消元,使得第k列除G(k,k)==1外其他位置为0
ind = find(H(:, k) ~= 0)';
for x = ind
if x ~= k
H(x, vec) = rem(H(x, vec) + H(k, vec), 2);
end
end
end
P=H(:,rows+1:cols);
assignin('base','H',H);
assignin('base','rearranged_cols',rearranged_cols);
%第三部分::高斯消元法编码
% 高斯消元的步骤
%设H=[A | B] ==========> [I | P]
% u=[c | s] %c是校验位 s是信息位 u是输出码字
%∵ H*u' = u*H' = 0
%代入得:
% _ _
% | c' |
% [I | P]| | = 0
% | s' |
% - -
%∴I*c' + P*s' = 0
%∴I*c' = P*s' (在GF(2)上)
%∴ c' = P*s'
%再由u=[c | s]即可得到编码后的码字。
%如果高斯消元过程中进行了列交换,
%则只需记录列交换,并以相反次序对编码后的码字同样进行列交换即可。
%解码时先求出u,再进行列交换得到uu=[c | s],后面部分即是想要的信息。
c=P*s; % c=P*s';
c=mod(c, 2); % c=mod(c, 2);
u1=[c s]; % u1=[c' s];
rows=length(rearranged_cols);
for i=rows:-1:1
if rearranged_cols(i)~=0
temp=u1(i);
u1(i)=u1(rearranged_cols(i));
u1(rearranged_cols(i))=temp;
end
end
u=u1;
u=u(:)
assignin('base','u',u);
amp=1;
tx_waveform=bpsk(u,amp);
rx_waveform=awgn(tx_waveform,10);
%需要将之前的H和rearranged_cols和u参数传递进来
dim=size(H);
rows=dim(1);
cols=dim(2);
vhat(1,1:cols)=0;
zero(1,1:rows)=0;
s=struct('qmn0',0,'qmn1',0,'dqmn',0,'rmn0',0,'rmn1',0,'qn0',0,'qn1',0,'alphamn',1);
%associate this structure with all non zero elements of H
amp=1;
scale(1:length(u))=1;
SNR=10;
%Prior Probabilities
%先验概率
pl1=1./(1+exp(-2*amp*rx_waveform.*scale/(SNR/2)));
pl0=1-pl1;
%initialization
%初始化:对特定的信道预设信息比特的先验概率。
newh(1:rows, 1:cols)=s;
[h1i h1j]=find(H==1);
h1num=length(h1i);
for i=1:h1num
newh(h1i(i),h1j(i)).qmn0=pl0(h1j(i));
newh(h1i(i),h1j(i)).qmn1=pl1(h1j(i));
end
%迭代次数100
for iteration=1:100
%horizontal step
%横向步骤:由信息节点的先验概率按置信传播算法得出各校验节点的后验概率。
for i=1:h1num %计算概率差
newh(h1i(i),h1j(i)).dqmn=newh(h1i(i),h1j(i)).qmn0-newh(h1i(i),h1j(i)).qmn1;
end
for i=1:rows
colind=find(h1i==i);
colnum=length(colind);
for j=1:colnum
drmn=1;
for k=1:colnum
if k~=j
drmn=drmn*newh(i,h1j(colind(k))).dqmn;
end
end
newh(i,h1j(colind(j))).rmn0=(1+drmn)/2;
newh(i,h1j(colind(j))).rmn1=(1-drmn)/2;
end
end
%vertical step
%纵向步骤:由校验节点的后验概率推算出信息节点的后验概率。
评论0