1% lpdc_fsk_lib.m
2% April 2015
3%
4% Library version of ldpc4.m written by vk5dsp.  Application is high bit rate
5% balloon telemtry
6%
7% LDPC demo
8% Call the CML routines and simulate one set of SNRs.
9% This fucntion is an updated version of ldpc3() which uses less
10% of the CML functions
11%
12% sim_in the input parameter structure
13% sim_out contains BERs and other stats for each value of SNR
14% resfile is the result file
15%
16
171;
18
19function sim_out = ldpc5(sim_in, resfile, testmode, genie_Es, logging=0);
20
21    if nargin<4, testmode = 0; end
22    estEsN0 = 0;
23
24    HRA       = sim_in.HRA;
25    framesize = sim_in.framesize;
26    rate      = sim_in.rate;
27    mod_order = sim_in.mod_order;
28
29    Lim_Ferrs = sim_in.Lim_Ferrs;
30    Ntrials   = sim_in.Ntrials;
31    Esvec     = sim_in.Esvec;
32
33    demod_type = 0;
34    decoder_type = 0;
35    max_iterations = 100;
36    code_param = ldpc_init(HRA, mod_order);
37    bps = code_param.bits_per_symbol;
38
39
40    if (logging)
41       fod = fopen('decode.log', 'w');
42       fwrite(fod, 'Es estEs Its  secs \n');
43    end
44
45
46    for ne = 1:length(Esvec)
47        Es = Esvec(ne);
48        EsNo = 10^(Es/10);
49
50
51        Terrs = 0;  Tbits =0;  Ferrs =0;
52        for nn = 1: Ntrials
53
54          data = round( rand( 1, code_param.data_bits_per_frame ) );
55          codeword = ldpc_encode(code_param, data);
56
57          code_param.code_bits_per_frame = length( codeword );
58          Nsymb = code_param.code_bits_per_frame/bps;
59
60            if testmode==1
61               f1 = fopen("dat_in2064.txt", "w");
62               for k=1:length(data);  fprintf(f1, "%u\n", data(k)); end
63               fclose(f1);
64               system("./ra_enc");
65
66               load("dat_op2064.txt");
67               pbits = codeword(length(data)+1:end);   %   print these to compare with C code
68               dat_op2064(1:16)', pbits(1:16)
69               differences_in_parity =  sum(abs(pbits - dat_op2064'))
70               pause;
71            end
72
73
74            % modulate
75            % s = Modulate( codeword, code_param.S_matrix );
76            s= 1 - 2 * codeword;
77            code_param.symbols_per_frame = length( s );
78
79            variance = 1/(2*EsNo);
80            noise = sqrt(variance)* randn(1,code_param.symbols_per_frame);
81            %  +  j*randn(1,code_param.symbols_per_frame) );
82            r = s + noise;
83            Nr = length(r);
84
85            [detected_data Niters] = ldpc_decode(r, code_param, max_iterations, decoder_type);
86
87            error_positions = xor( detected_data(1:code_param.data_bits_per_frame), data );
88            Nerrs = sum( error_positions);
89
90            t = clock;   t =  fix(t(5)*60+t(6));
91            if (logging)
92              fprintf(fod, ' %3d  %4d\n', Niters, t);
93              end
94
95            if Nerrs>0, fprintf(1,'x'),  else fprintf(1,'.'),  end
96            if (rem(nn, 50)==0),  fprintf(1,'\n'),  end
97
98            if Nerrs>0,  Ferrs = Ferrs +1;  end
99            Terrs = Terrs + Nerrs;
100            Tbits = Tbits + code_param.data_bits_per_frame;
101
102            if Ferrs > Lim_Ferrs, disp(['exit loop with #cw errors = ' ...
103                    num2str(Ferrs)]);  break,  end
104        end
105
106        TERvec(ne) = Terrs;
107        FERvec(ne) = Ferrs;
108        BERvec(ne) = Terrs/ Tbits;
109        Ebvec = Esvec - 10*log10(code_param.bits_per_symbol * rate);
110
111        cparams= [code_param.data_bits_per_frame  code_param.symbols_per_frame ...
112            code_param.code_bits_per_frame];
113
114        sim_out.BERvec = BERvec;
115        sim_out.Ebvec = Ebvec;
116        sim_out.FERvec = FERvec;
117        sim_out.TERvec  = TERvec;
118        sim_out.cpumins = cputime/60;
119
120        if nargin > 2
121            save(resfile,  'sim_in',  'sim_out',  'cparams');
122            disp(['Saved results to ' resfile '  at Es =' num2str(Es) 'dB']);
123        end
124      end
125end
126
127
128function code_param = ldpc_init(HRA, mod_order)
129  code_param.bits_per_symbol = log2(mod_order);
130  [H_rows, H_cols] = Mat2Hrows(HRA);
131  code_param.H_rows = H_rows;
132  code_param.H_cols = H_cols;
133  code_param.P_matrix = [];
134  code_param.data_bits_per_frame = length(code_param.H_cols) - length( code_param.P_matrix );
135  code_param.symbols_per_frame = length(HRA);
136end
137
138
139function codeword = ldpc_encode(code_param, data)
140      codeword = LdpcEncode( data, code_param.H_rows, code_param.P_matrix );
141endfunction
142
143
144% Takes soft decision symbols (e.g. output of 2fsk demod) and converts
145% them to LLRs.  Note we calculate mean and var manually instead of
146% using internal functions.  This was required to get a bit exact
147% results against the C code.
148
149function llr = sd_to_llr(sd)
150    sd = sd / mean(abs(sd));
151    x = sd - sign(sd);
152    sumsq = sum(x.^2);
153    summ = sum(x);
154    mn = summ/length(sd);
155    estvar = sumsq/length(sd) - mn*mn;
156    estEsN0 = 1/(2* estvar + 1E-3);
157    llr = 4 * estEsN0 * sd;
158endfunction
159
160
161% LDPC decoder - note it estimates EsNo from received symbols
162
163function [detected_data Niters] = ldpc_decode(r, code_param, max_iterations, decoder_type)
164  % in the binary case the LLRs are just a scaled version of the rx samples ..
165
166 #{
167  r = r / mean(abs(r));       % scale for signal unity signal
168  estvar = var(r-sign(r));
169  estEsN0 = 1/(2* estvar + 1E-3);
170  input_decoder_c = 4 * estEsN0 * r;
171 #}
172  llr = sd_to_llr(r);
173
174  [x_hat, PCcnt] = MpDecode(llr, code_param.H_rows, code_param.H_cols, ...
175                            max_iterations, decoder_type, 1, 1);
176  Niters = sum(PCcnt!=0);
177  detected_data = x_hat(Niters,:);
178
179  if isfield(code_param, "c_include_file")
180    ldpc_gen_h_file(code_param, max_iterations, decoder_type, llr, x_hat, detected_data);
181  end
182end
183
184
185% One application of FSK LDPC work is SSTV.  This function generates a
186% simulated frame for testing
187
188function frame_rs232 = gen_sstv_frame
189  load('H2064_516_sparse.mat');
190  HRA = full(HRA);
191  mod_order = 2;
192  code_param = ldpc_init(HRA, mod_order);
193
194  % generate payload data bytes and checksum
195
196  data = floor(rand(1,256)*256);
197  %data = zeros(1,256);
198  checksum = crc16(data);
199  data = [data hex2dec(checksum(3:4)) hex2dec(checksum(1:2))];
200
201  % unpack bytes to bits and LPDC encode
202
203  mask = 2.^(7:-1:0);       % MSB to LSB unpacking to match python tx code.
204  unpacked_data = [];
205  for b=1:length(data)
206    unpacked_data = [unpacked_data bitand(data(b), mask) > 0];
207  end
208  codeword = [ldpc_encode(code_param, unpacked_data) 0 0 0 0];  % pad with 0s to get integer number of bytes
209
210  % pack back into bytes to match python code
211
212  lpacked_codeword = length(codeword)/8;
213  packed_codeword = zeros(1,lpacked_codeword);
214  for b=1:lpacked_codeword
215    st = (b-1)*8 + 1;
216    packed_codeword(b) = sum(codeword(st:st+7) .* mask);
217  end
218
219  % generate header bits
220
221  header = [hex2dec('55')*ones(1,16) hex2dec('ab') hex2dec('cd') hex2dec('ef') hex2dec('01')];
222
223  % now construct entire unpacked frame
224
225  packed_frame = [header packed_codeword];
226  mask = 2.^(0:7);        % LSB to MSB packing for header
227  lpacked_frame = length(packed_frame);
228  frame = [];
229  for b=1:lpacked_frame
230    frame = [frame bitand(packed_frame(b), mask) > 0];
231  end
232
233  % insert rs232 framing bits
234
235  frame_rs232 = [];
236  for b=1:8:length(frame)
237    frame_rs232 = [frame_rs232 0 frame(b:b+7) 1];
238  end
239
240  %printf("codeword: %d unpacked_header: %d frame: %d frame_rs232: %d \n", length(codeword), length(unpacked_header), length(frame), length(frame_rs232));
241endfunction
242
243
244% calculates and compares the checksum of a SSTV frame, that has RS232
245% start and stop bits
246
247function checksum_ok = sstv_checksum(frame_rs232)
248  l = length(frame_rs232);
249  expected_l = (256+2)*10;
250  assert(l == expected_l);
251
252  % extract rx bytes
253
254  rx_data = zeros(1,256);
255  mask = 2.^(0:7);          % LSB to MSB
256  k = 1;
257  for i=1:10:expected_l
258    rx_bits = frame_rs232(i+1:i+8);
259    rx_data(k) = sum(rx_bits .* mask);
260    k++;
261  end
262
263  % calc rx checksum and extract tx checksum
264
265  rx_checksum = crc16(rx_data(1:256));
266  tx_checksum = sprintf("%02X%02X", rx_data(258), rx_data(257));
267  %printf("tx_checksum: %s rx_checksum: %s\n", tx_checksum, rx_checksum);
268  checksum_ok = strcmp(tx_checksum, rx_checksum);
269endfunction
270