用matlab来仿真ldpc 有构造检验矩阵 编码 译码的程序但是不知道怎么连贯的仿真

用matlab来仿真ldpc 有构造检验矩阵 编码 译码的程序但是不知道怎么连贯的仿真

sigm=0.1:0.1:5;
ber=zeros(size(sigm));
for l=1:100
for w=1:length(sigm)
i=load('H.txt');
a=i(:,1);
b=i(:,2);
a=a+1;
b=b+1;
H=zeros(130,260);
idx=sub2ind(size(H),a',b');
H(idx)=1;%generate the H matrix
x=zeros(1,260);
y=2*x-1;
z=y+sigm(w)*randn(1,260);
p1=ones(1,260)./(1+exp(-2*z./sigm(w)));
p0=1-p1;
qij1=H.*repmat(p1,130,1);
qij0=H.*repmat(p0,130,1);%intial step
for i=1:100
prdct=ones(130,1);
rji0=zeros(size(qij1));
tqij1=ones(size(qij1))-2*qij1;
for i=1:130
for j=1:260
if tqij1(i,j)~=0
prdct(i)=prdct(i)*tqij1(i,j);
else
prdct(i)=prdct(i);
end
end
end
for i=1:130
for j=1:260
if tqij1(i,j)~=0
rji0(i,j)=prdct(i)/tqij1(i,j);
else
rij0(i,j)=0;
end
end
end
rji0=1/2+1/2*rji0;
rji1=1-rji0;%horizontal step
prdct1=ones(1,260);
for j=1:260
for i=1:130
if rji1(i,j)~=0
prdct1(j)=prdct1(j)*rji1(i,j);
else
prdct1(j)=prdct1(j);
end
end
end
QQ1=p1.*prdct1;
prdct0=ones(1,260);
for j=1:260
for i=1:130
if rji0(i,j)~=0
prdct0(j)=prdct0(j)*rji0(i,j);
else
prdct0(j)=prdct0(j);
end
end
end
QQ0=p0.*prdct0;
k=ones(size(QQ1))./(QQ1+QQ0);
Q0=k.*QQ0;
Q1=k.*QQ1;
prdct2=ones(1,260);
qqij0=zeros(size(qij0));
for j=1:260
for i=1:130
if rji0(i,j)~=0
prdct2(j)=prdct2(j)*rji0(i,j);
else
prdct2(j)=prdct2(j);
end
end
end
for j=1:260
for i=1:130
if rji0(i,j)~=0
qqij0(i,j)=prdct2(j)/rji0(i,j);
else
qqij0(i,j)=0;
end
end
end
qqqij0=repmat(p0,130,1).*qqij0;
prdct3=ones(1,260);
qqij1=zeros(size(qij1));
for j=1:260
for i=1:130
if rji1(i,j)~=0
prdct3(j)=prdct3(j)*rji1(i,j);
else
prdct3(j)=prdct3(j);
end
end
end
for j=1:260
for i=1:130
if rji1(i,j)~=0
qqij1(i,j)=prdct3(j)/rji1(i,j);
else
qqij1(i,j)=0;
end
end
end
qqqij1=repmat(p1,130,1).*qqij1;
kk=ones(size(qij0))./(qqqij0+qqqij1);
qij0=kk.*qqqij0;
qij1=kk.*qqqij1;%vertical step
zz=(sign(Q1-Q0)+1)/2;%decide the code should be either 1 or 0
if rem(zz*H',2)==0
break;
end %calculate the syndrome
end
ber(w)=length(find(zz))/(260+ber(w));
end
end
snr=ones(size(sigm))./sigm;
ber=ber./(ones(size(ber))*100);
plot(20*log(snr),ber,'*-');
xlabel('SNR');
ylabel('BER');
我自己编的,用SPA解码,循环一百次观察信噪比和误码率的关系,复杂度非常高,我是菜鸟.