N=input('N=')%输入信源符号的个数
s=0l=0H=0
for i=1:N
fprintf('第%d个',i)
p(i)=input('p=')%输入信源符号概率分布矢量,p(i)<1
if p(i)<=0
error('不符合概率分布')
end
s=s+p(i)
H=H+(- p(i)*log2(p(i)))%计算信源信息熵
end
if (s<=0.999999||s>=1.000001)
error('不符合概率分布')
end
tic
for i=1:N-1 %按概率分布大小对信源排序
for j=i+1:N
if p(i)<p(j)
m=p(j)p(j)=p(i)p(i)=m
end
end
end
x=f1(1,N,p,1)
for i=1:N %计算平均码长
L(i)=length(find(x(i,:)))
l=l+p(i)*L(i)
end
n=H/l%计算编码效率
fprintf('按概率降序排列的码字:\n')
disp(x) %显示按概率降序排列的码字
fprintf('平均码长:\n')
disp(l)% 显示平均码长
fprintf('信源信息熵:\n')
disp(H)%显示信源信息熵
fprintf('编码效率:\n')
disp(n) %显示编码效率
fprintf('计算耗时time= %f\n',toc)
再建立两个M文件：%函数f1存放于f1.m
function x=f1(i,j,p,r)
global x
x=char(x)
if(j<=i)
return
else
q=0
for t=i:j %对于区间[i,j]自上而下求累加概率值
q=p(t)+qy(t)=q
end
for t=i:j%把所有自上而下的累加概率值与该区间总概率值减该累加概率值之差取绝对值存在一数组
v(t)=abs(y(t)-(q-y(t)))
end
for t=i:j
if(v(t)==min(v)) %求该数组中最小的一个值来确定分界点位置
for k=i:t%赋值码字
x(k,r)='0'
end
for k=(t+1):j
x(k,r)='1'
end
d=t
f1(i,d,p,r+1)%递归调用及相互调用
f2(d+1,j,p,r+1)
f1(d+1,j,p,r+1)
f2(i,d,p,r+1)
else
end
end
end
return第二个：%函数f2存放于f2.m
function x=f2(i,j,p,r)
global x
x=char(x)
if(j<=i)
return
else
q=0
for t=i:j %对于区间[i,j]自上而下求累加概率值
q=p(t)+qy(t-i+1)=q
end
for t=1:j-(i-1)%把所有自上而下的累加概率值与该区间总概率值减该累加概率值之差取绝对值存在一数组
v(t)=abs(y(t)-(q-y(t)))
end
for t=1:j-(i-1)
if(v(t)==min(v)) %求该数组中最小的一个值来确定分界点位置
d=t+i-1
for k=i:d %赋值码字
x(k,r)='0'
end
for k=(d+1):j
x(k,r)='1'
end
f2(d+1,j,p,r+1)%递归调用及相互调用
f1(i,d,p,r+1)
f2(i,d,p,r+1)
f1(d+1,j,p,r+1)
else
end
end
end
return
LDPC码编码是在通信系统的发送端进行的，在接收端进行相应的译码，这样才能实现编码的纠错。LDPC 码由于其奇偶校验矩阵的稀疏性，使其存在高效的译码算法，其复杂度与码长成线性关系，克服了分组码在码长很大时，所面临的巨大译码算法复杂度问题，使长码分组的应用成为可能。而且由于校验矩阵稀疏，使得在长码时，相距很远的信息比特参与统一校验，这使得连续的突发差错对译码的影响不大，编码本身就具有抗突发错误的特性。
LDPC码的译码算法种类很多，其中大部分可以被归结到信息传递〔Mesaseg Prpagation，MP)算法集中。这一类译码算法由于具有良好的性能和严格的数学结构，使得译码性能的定量分析成为可能，因此特别受到关注。MP算法集中的置信传播(BP)算法是Gallager提出的一种软输入迭代译码算法，具有最好的性能。如果我们首先理解并掌握了一些很简单的硬判决算法后，对BP算法的理解会更加容易。同时，通过一些常用的数学手段，我们可以对BP译码算法作一些简化，从而在一定的性能损失内获得对运算量和存储量需求的降低。
sigm=0.1:0.1:5ber=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解码，循环一百次观察信噪比和误码率的关系，复杂度非常高，我是菜鸟。
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）