close all
clear all
%% H矩阵生成
load Matrix(2016,1008)Block56.mat
z=56;
m=1008;
n=2016;
k=n-m;
mb=m/z;
nb=n/z;
H=[];
for i=1:mb
H_row=[];
for j=1:nb
Hij=zeros(z,z);
if (H_block(i,j)~=0)
Hij(1,H_block(i,j))=1;
for q=2:z
Hij(q,:)=circshift(Hij(1,:),q-1); %循环移位
end
end
H_row=[H_row Hij]; %矩阵行拼接
end
H=[H;H_row]; %拼接成完整矩阵
end
H(1,m)=0; %调整a的右上角的值
Hs=H(:,m+1:n);
%% 不同信噪比下性能仿真
SNR=-1:1:2;
BER=zeros(1,length(SNR));
for snr=1:length(SNR)
if(SNR(snr)==2) %仿真停止条件
num_stop=5;
else
num_stop=50;
end
num_eb=0; %初始化错误比特数
num_ec=0; %初始化错误码数
num_code=0; %初始化总码数
while(num_ec<num_stop)
s=randi([0 1],1,k);
x=s*Hs';
%% LDPC编码
p=zeros(1,m);
for i=1:z
if (i==1)
p(i)=mod(x(i),2);
elseif(i>1&&i<=z)
p(i)=mod(x(i)+p((mb-1)*z+i-1),2);
end
for j=1:mb-1
p(j*z+i)=mod(x(j*z+i)+p((j-1)*z+i),2);
end
end
s_out=[p s];
%% 调制和通过信道
s_bpsk=-2*s_out+1; %bpsk调制
s_noise=awgn(s_bpsk,SNR(snr),'measured'); %通过AWGN信道
u0=4*s_noise(1:n).*10^((SNR(snr)-3)/10); %初始化
uji=zeros(m,n);
vij=zeros(n,m);
itera=0; %迭代次数
while(itera<30)
%% 和积算法解码
for i=1:n
v2c=find(H(:,i)==1)';
for j=1:length(v2c)
tmp=v2c;
tmp(j)=[];
temp=sum(uji(tmp(1:length(tmp)),i));
vij(i,v2c(j))=u0(i)+temp;
end
end
for j=1:m
v2c=find(H(j,:)==1);
for i=1:length(v2c)
tmp=v2c;
tmp(i)=[];
temp1=prod(tanh(vij(tmp(1:length(tmp)),j)/2));
uji(j,v2c(i))=2*atanh(temp1);
end
end
v=sum(uji)+u0;
v(v>0)=0;
v(v<0)=1;
if mod(v*H',2)==0
break
else
itera=itera+1;
end
end
%% 最小和算法
% for i=1:n
% v2c=find(H(:,i)==1)';
% for j=1:length(v2c)
% tmp=v2c;
% tmp(j)=[];
% temp=sum(uji(tmp(1:length(tmp)),i));
% vij(i,v2c(j))=u0(i)+temp;
% end
% end
%
% for j=1:m
% v2c=find(H(j,:)==1);
% for i=1:length(v2c)
% tmp=v2c;
% tmp(i)=[];
% minv=min(abs(vij(tmp(1:length(tmp)),j)));
% sym=prod(sign(vij(tmp(1:length(tmp)),j)));
% uji(j,v2c(i))=sym*minv;
% end
% end
%
% v=sum(uji)+u0;
% v(v>0)=0;
% v(v<0)=1;
%
% if mod(v*H',2)==0
% break
% else
% itera=itera+1;
% end
%
% end
%% Normalized Min-sum
% a=0.7;
% for i=1:n
% v2c=find(H(:,i)==1)';
% for j=1:length(v2c)
% tmp=v2c;
% tmp(j)=[];
% temp=sum(uji(tmp(1:length(tmp)),i));
% vij(i,v2c(j))=u0(i)+temp;
% end
% end
%
% for j=1:m
% v2c=find(H(j,:)==1);
% for i=1:length(v2c)
% tmp=v2c;
% tmp(i)=[];
% minv=min(abs(vij(tmp(1:length(tmp)),j)));
% sym=prod(sign(vij(tmp(1:length(tmp)),j)));
% uji(j,v2c(i))=sym*minv*a;
% end
% end
%
% v=sum(uji)+u0;
% v(v>0)=0;
% v(v<0)=1;
%
% if mod(v*H',2)==0
% break
% else
% itera=itera+1;
% end
%
% end
%% Offset Min-sum
% b=0.5;
% for i=1:n
% v2c=find(H(:,i)==1)';
% for j=1:length(v2c)
% tmp=v2c;
% tmp(j)=[];
% temp=sum(uji(tmp(1:length(tmp)),i));
% vij(i,v2c(j))=u0(i)+temp;
% end
% end
%
% for j=1:m
% v2c=find(H(j,:)==1);
% for i=1:length(v2c)
% tmp=v2c;
% tmp(i)=[];
% minv=min(abs(vij(tmp(1:length(tmp)),j)));
% sym=prod(sign(vij(tmp(1:length(tmp)),j)));
% maxv=max((minv-b),0);
% uji(j,v2c(i))=sym*maxv;
% end
% end
%
% v=sum(uji)+u0;
% v(v>0)=0;
% v(v<0)=1;
%
% if mod(v*H',2)==0
% break
% else
% itera=itera+1;
% end
%
% end
%% 统计误码率
ber=sum(abs(v(m+1:n)-s));
num_eb=ber+num_eb;
if ber~=0
cer=1;
else
cer=0;
end
num_ec=cer+num_ec;
num_code=num_code+1;
end
BER(snr)=num_eb/k/num_code;
end
评论1