1% Test MFSK at symbol level in CML using RA LDPC codes,  Bill, June 2020
2% Simulate in AWGM and plot BERs.    Assumes that CML is installed!
3%
4% If required setup numb of codewords (Ncw) and channel bits/sym (bps, 1 to 4)
5% may also select FEC code with Ctype (1 to 3); use plt to for debug plots of LLR values
6
7% July 1 version allows M=8 and pads required vectors if required to makeup 3rd bit
8
9ldpc;
10
11%rand('seed',1);
12%randn('seed',1);
13format short g
14more off
15init_cml();
16
17if exist('Ncw')==0,    Ncw=100, end     %setup defaults
18if exist('plt')==0,    plt=0;   end
19if exist('bps')==0,    bps=4;   end
20if exist('Ctype')==0,  Ctype=1, end
21
22if Ctype==1
23   load H_256_768_22.txt;   HRA = H_256_768_22; % rate 1/3
24elseif Ctype==2
25   load H_256_512_4.mat;   HRA=H;  % K=256,  rate 1/2 code
26   % above code might be improved -- but still works better than rate 1/3
27elseif Ctype==3
28   load HRAa_1536_512.mat; % rate 3/4,  N=2k code
29else
30   error('bad Ctype');
31end
32
33M=2^bps;   nos =0;    clear res
34modulation = 'FSK'; mod_order=M; mapping = 'gray';
35code_param = ldpc_init_user(HRA, modulation, mod_order, mapping);
36
37[Nr Nc] = size(HRA);
38Nbits = Nc - Nr;
39Krate = (Nc-Nr)/Nc
40framesize = Nc;
41%{
42[H_rows, H_cols] = Mat2Hrows(HRA);
43code_param.H_rows = H_rows;
44code_param.H_cols = H_cols;
45code_param.P_matrix = [];
46%}
47
48
49S =CreateConstellation('FSK', M);
50if M==2,   Ebvec=[7:0.2: 8.7],   end
51if M==4,   Ebvec=[4:0.25: 7],  end
52if M==8,   Ebvec =[3.5: 0.25: 5.5];  end
53if M==16,  Ebvec=[2.5:0.25:4.8];   end
54
55disp(['Symbol-based ' num2str(M) 'FSK sim with K=' ...
56       num2str(Nbits) ', code rate=' num2str(Krate) ', #CWs=' num2str(Ncw)])
57
58
59
60% if M=8, for 3 bits/symbol, may need to pad codeword with bits ...
61if floor(Nc/bps)*bps ~= Nc
62   Npad = ceil(Nc/bps) *bps-Nc
63   disp('padding codeword')
64else
65   Npad=0;
66end
67
68for Eb = Ebvec
69
70    Ec =  Eb + 10*log10(Krate);
71    Es =  Ec + 10*log10(bps);
72    Eslin = 10^(Es/10);       %Es/N0 = 1/2k_n^2
73
74    Terrs =0;
75    for nn = 1:Ncw
76
77        txbits = randi(2,1,Nbits) -1;
78
79        codeword = LdpcEncode( txbits, code_param.H_rows, code_param.P_matrix );
80        code_param.code_bits_per_frame = length( codeword );
81        code_param.data_bits_per_frame = length(txbits);
82        Nsymb = (code_param.code_bits_per_frame+Npad)/bps;
83
84        if Npad; codeword = [codeword   zeros(1,Npad)]; end
85        Tx = Modulate(codeword, S);
86
87        kn = sqrt(1/(2*Eslin));
88        Rx = Tx + kn * (randn(size(Tx)) + j*randn(size(Tx)));
89
90        SNRlin = Eslin;  % Valenti calls this SNR, but seems to be Es/N0
91        symL = DemodFSK(Rx, SNRlin, 2);    %demod type is nonCOH, without estimate amplitudes
92        bitL = Somap(symL);
93
94	if Npad, bitL(end-Npad+1:end)=[]; end
95        if plt>0, figure(110);   hist(bitL);   title('bit LLRs')
96            figure(111);   hist(bitL);   title('Sym Ls'),  pause,
97        end
98        max_it =100;   decoder_type =0;
99
100        [x_hat, PCcnt] = MpDecode( -bitL, code_param.H_rows, code_param.H_cols, ...
101            max_it, decoder_type, 1, 1);
102        Niters = sum(PCcnt~=0);
103        detected_data = x_hat(Niters,:);
104        error_positions = xor( detected_data(1:code_param.data_bits_per_frame), txbits );
105
106        Nerrs = sum( error_positions);
107        if plt>1, figure(121);  plot(error_positions);  Nerrs,  end
108        Terrs = Terrs + Nerrs;
109    end
110
111    BER = Terrs/ (Ncw*Nbits);
112
113    %HDs = (sign(bitL)+1)/2;
114    %NerrsHD  = sum(codeword~=HDs);
115    %BER_HD = Nerrs/Nbits;
116
117    nos = nos+1;
118    disp('Eb    Nerrs     BER')
119    res(nos, :) = [Eb, Terrs,  BER]
120
121end
122figure(91)
123semilogy(res(:,1), res(:,3), '-x'); grid on; hold on;
124%semilogy(res(:,1), berawgn(res(:,1), 'fsk', M, 'noncoherent'), 'g');
125title([num2str(M) 'FSK BER with LDPC FEC'])
126xlabel('Eb/N0'); ylabel('BER')
127
128
129
130
131
132
133