FFT+LDPC

tic

close all

clear all

clc

%ミ LDPC H

onePerCol = 3;

onePerRow = 6;

coderate = (onePerRow-onePerCol)/onePerRow;

%gallagerLDPC痻皚

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);

%GF4H痻皚

GF4H = H.*randint(M,N,[1,3]);

%安砞message

bit_num = 1200;

message = zeros(1,bit_num);

%安砞竒筁LDPCcodeword

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