tic
close all
clear all
clc
%ミ LDPC H
onePerCol = 3;
onePerRow = 6;
coderate = (onePerRow-onePerCol)/onePerRow;
%gallagerLDPC痻皚
k = 100;
H1 = zeros(k,k*onePerRow);
for i = 1:k
H1(i,(i-1)*onePerRow+1:i*onePerRow) = ones(1,onePerRow);
end
H22 = randperm(length(H1));
H2 = H1(:,H22(1,:));
H33 = randperm(length(H1));
H3 = H1(:,H33(1,:));
H = [H1;H2;H3];
[M N] = size(H);
%GF4H痻皚
GF4H = H.*randint(M,N,[1,3]);
%安砞message
bit_num = 1200;
message = zeros(1,bit_num);
%安砞竒筁LDPCcodeword
cw = zeros(1,N);
%兜跑计﹍
WER = [];
Symbol_error_rate = [];
SNR = 0.5:0.1:1;
iteration = 30;
loop_num = 20;
for ebn0 = SNR
ebn0
noe = 0;
error_num = 0;
error_block = 0;
for loop = 1:loop_num
%BPSK mod
bpskmod = (message*2-1);
%AWGN
ebi = 1/coderate;
ebc = ebi*coderate;
n0 = ebi.*10.^(-ebn0./10);
sigma_2 = n0/2;
sigma = sqrt(sigma_2);
AWGN=sigma*(randn(1,length(bpskmod)));
%钡Μ
yi = bpskmod + AWGN;
%==========FFT decode ====================
%兜跑计﹍
R0 = zeros(M,N);
R1 = zeros(M,N);
R2 = zeros(M,N);
R3 = zeros(M,N);
F_qnm0 = zeros(M,N);
F_qnm1 = zeros(M,N);
F_qnm2 = zeros(M,N);
F_qnm3 = zeros(M,N);
%﹍
%0 & 1诀瞯
prob1 = ones(size(yi))./(1+exp(-2*yi./sigma_2));
prob0 = 1 - prob1;
prob1a = prob1(1:2:end);
prob1b = prob1(2:2:end);
prob0a = prob0(1:2:end);
prob0b = prob0(2:2:end);
%衡GF4 0~3诀瞯
GF4_0 = prob0a.*prob0b;
GF4_1 = prob0a.*prob1b;
GF4_2 = prob1a.*prob0b;
GF4_3 = prob1a.*prob1b;
GF4 = [GF4_0;GF4_1;GF4_2;GF4_3];
%FFT
L_qnm0 = H.*repmat(GF4_0,M,1);
L_qnm1 = H.*repmat(GF4_1,M,1);
L_qnm2 = H.*repmat(GF4_2,M,1);
L_qnm3 = H.*repmat(GF4_3,M,1);
[x1 y1] = find(GF4H == 1);
[x2 y2] = find(GF4H == 2);
[x3 y3] = find(GF4H == 3);
%=======================================
for g = 1:iteration
repGF4_0 = L_qnm0;
repGF4_1 = L_qnm1;
repGF4_2 = L_qnm2;
repGF4_3 = L_qnm3;
for i = 1:length(x2)
L_qnm1(x2(i),y2(i)) = repGF4_3(x2(i),y2(i));
L_qnm2(x2(i),y2(i)) = repGF4_1(x2(i),y2(i));
L_qnm3(x2(i),y2(i)) = repGF4_2(x2(i),y2(i));
end
for i = 1:length(x3)
L_qnm1(x3(i),y3(i)) = repGF4_2(x3(i),y3(i));
L_qnm2(x3(i),y3(i)) = repGF4_3(x3(i),y3(i));
L_qnm3(x3(i),y3(i)) = repGF4_1(x3(i),y3(i));
end
%C_node(Horizontal) 癟穝
for i = 1:M
c1 = find(H(i,:));
trmatrix = [1 1;1 -1];
FFTmatrix = [trmatrix trmatrix;trmatrix -trmatrix];
F_qnm = [];
for k = 1:length(c1)
F_qnm(:,k) = FFTmatrix*[L_qnm0(i,c1(k));L_qnm1(i,c1(k));L_qnm2(i,c1(k));L_qnm3(i,c1(k))];
end
for k = 1:length(c1)
r_mn = 1;
R_mn = 1;
for l = 1:length(c1)
if l ~= k
r_mn = r_mn.*F_qnm(:,l);
end
R_mn = 0.25*FFTmatrix*r_mn;
R0(i,c1(k)) = R_mn(1,1);
R1(i,c1(k)) = R_mn(2,1);
R2(i,c1(k)) = R_mn(3,1);
R3(i,c1(k)) = R_mn(4,1);
end
end
end
deperR0 = R0;
deperR1 = R1;
deperR2 = R2;
deperR3 = R3;
for i = 1:length(x2)
R1(x2(i),y2(i)) = deperR2(x2(i),y2(i));
R2(x2(i),y2(i)) = deperR3(x2(i),y2(i));
R3(x2(i),y2(i)) = deperR1(x2(i),y2(i));
end
for i = 1:length(x3)
R1(x3(i),y3(i)) = deperR3(x3(i),y3(i));
R2(x3(i),y3(i)) = deperR1(x3(i),y3(i));
R3(x3(i),y3(i)) = deperR2(x3(i),y3(i));
end
%V_node(Vertical)癟穝
for j = 1:N
r1 = find(H(:,j));
for k = 1:length(r1)
beta0 = 1;
beta1 = 1;
beta2 = 1;
beta3 = 1;
for l = 1:length(r1)
if l ~= k
beta0 = beta0*R0(r1(l),j);
beta1 = beta1*R1(r1(l),j);
beta2 = beta2*R2(r1(l),j);
beta3 = beta3*R3(r1(l),j);
end
end
q_nm0(r1(k),j) = GF4_0(j)*beta0;
q_nm1(r1(k),j) = GF4_1(j)*beta1;
q_nm2(r1(k),j) = GF4_2(j)*beta2;
q_nm3(r1(k),j) = GF4_3(j)*beta3;
L_qnm0(r1(k),j) = q_nm0(r1(k),j)./(q_nm0(r1(k),j) + q_nm1(r1(k),j) + q_nm2(r1(k),j) + q_nm3(r1(k),j));
L_qnm1(r1(k),j) = q_nm1(r1(k),j)./(q_nm0(r1(k),j) + q_nm1(r1(k),j) + q_nm2(r1(k),j) + q_nm3(r1(k),j));
L_qnm2(r1(k),j) = q_nm2(r1(k),j)./(q_nm0(r1(k),j) + q_nm1(r1(k),j) + q_nm2(r1(k),j) + q_nm3(r1(k),j));
L_qnm3(r1(k),j) = q_nm3(r1(k),j)./(q_nm0(r1(k),j) + q_nm1(r1(k),j) + q_nm2(r1(k),j) + q_nm3(r1(k),j));
end
%=========================codeword浪琩===============
q0 = GF4_0(j)*prod(R0(r1(k),j));
q1 = GF4_1(j)*prod(R1(r1(k),j));
q2 = GF4_2(j)*prod(R2(r1(k),j));
q3 = GF4_3(j)*prod(R3(r1(k),j));
Q0 = q0/(q0+q1+q2+q3);
Q1 = q1/(q0+q1+q2+q3);
Q2 = q2/(q0+q1+q2+q3);
Q3 = q3/(q0+q1+q2+q3);
Qmatrix = [Q0;Q1;Q2;Q3];
[maxvalue maxGF] = max(Qmatrix);
if maxGF == 4
Q(j) = 3;
elseif maxGF == 3
Q(j) = 2;
elseif maxGF == 2
Q(j) = 1;
elseif maxGF == 1
Q(j) = 0;
end
end
if Q == cw
break
end
end
%璸衡岿粇瞯
noe = sum(Q ~= cw);
error_num = error_num + noe;
%璸衡WER
if noe == 0
error_block = error_block;
else
error_block = error_block+1;
end
% if error_num > 24000
% break
% end
end
nod = length(cw);
c_e = error_num/nod/loop_num;
Symbol_error_rate = [Symbol_error_rate c_e];
ber_bk = error_block/loop_num;
WER = [WER ber_bk];
Symbol_error_rate
WER
end
semilogy(SNR,Symbol_error_rate,'-*');
hold on
ber_th = (1/2)*erfc(sqrt(10.^(SNR/10)));
semilogy(SNR,ber_th,'r*-');
hold on
semilogy(SNR,WER,'g*-');
legend('Symbol error rate-36LDPC','berth','WER-36LDPC');
xlabel('EbN0');
title('FFT');
grid on
toc