1% fsk_lib_ldpc_rx.m
2%
3% LDPC coded 4FSK modem rx, reads 8 kHz 16 bit short raw file of real samples and demodulates
4
5function fsk_lib_ldpc_rx(filename, Rs=100, coderate=0.5)
6  fsk_lib_ldpc;
7
8  % set up LDPC code
9  init_cml();
10  if coderate == 0.5
11    load H_256_512_4.mat;
12  elseif coderate == 0.75
13    load HRAa_1536_512.mat; H=HRA;
14  else
15    disp("unknown code rate");
16  end
17  [states code_param] = fsk_lib_ldpc_init (H, Rs, Fs=8000);
18  n = code_param.coded_bits_per_frame; k = code_param.data_bits_per_frame;
19
20  % known transmitted bits for BER estimation
21  rand('seed',1);
22  data_bits = round(rand(1,code_param.data_bits_per_frame));
23  codeword_bits = LdpcEncode(data_bits, code_param.H_rows, code_param.P_matrix);
24
25  frx=fopen(filename,"rb"); rx = fread(frx, Inf, "short"); fclose(frx);
26
27  % freq estimator and demod
28  run_frames = floor(length(rx)/states.N)-1;
29  st = 1; f_log = []; rx_bits = []; rx_filt = []; SNRest_log = []; rx_timing_log = [];
30  for f=1:run_frames
31
32    % extract nin samples from input stream
33    nin = states.nin;
34    en = st + states.nin - 1;
35
36    % due to nin variations it's possible to overrun buffer
37    if en < length(rx)
38      sf = rx(st:en);
39      states = est_freq(states, sf, states.M);  states.f = states.f2;
40      [arx_bits states] = fsk_demod(states, sf);
41      rx_bits = [rx_bits arx_bits];
42      rx_filt = [rx_filt abs(states.f_int_resample)];
43      f_log = [f_log; states.f];
44      SNRest_log = [SNRest_log states.SNRest];
45      rx_timing_log = [rx_timing_log states.norm_rx_timing];
46      st += nin;
47    end
48  end
49
50  % count bit errors in test frames
51
52  num_frames=floor(length(rx_bits)/code_param.coded_bits_per_frame);
53  log_nerrs = []; num_frames_rx = 0; Tbits = Terrs = Tperr = Tpackets = 0;
54  uber = cber = 0.5; cper = 1;
55  for f=1:num_frames-1
56    st = (f-1)*code_param.coded_bits_per_frame + 1; en = (f+1)*code_param.coded_bits_per_frame;
57    states = ber_counter(states, codeword_bits, rx_bits(st:en));
58    log_nerrs = [log_nerrs states.nerr];
59    if states.ber_state num_frames_rx++; end
60
61    % Using sync provided by ber_counter() state machine for LDPC frame alignment
62    if states.ber_state
63      st_bit = (f-1)*code_param.coded_bits_per_frame + states.coarse_offset;
64      st_symbol = (st_bit-1)/states.bitspersymbol + 1;
65      en_symbol = st_symbol +  code_param.coded_bits_per_frame/states.bitspersymbol - 1;
66
67      % map FSK filter ouputs to LLRs, then LDPC decode (see also fsk_cml_sam.m)
68      symL = DemodFSK(1/states.v_est*rx_filt(:,st_symbol:en_symbol), states.SNRest, 1);
69      llr = -Somap(symL);
70      [x_hat, PCcnt] = MpDecode(llr, code_param.H_rows, code_param.H_cols, max_iterations=100, decoder_type=0, 1, 1);
71      Niters = sum(PCcnt~=0);
72      detected_data = x_hat(Niters,:);
73      Nerrs = sum(xor(data_bits, detected_data(1:code_param.data_bits_per_frame)));
74      Terrs += Nerrs;
75      Tbits += code_param.data_bits_per_frame;
76      if Nerrs Tperr++; end
77      Tpackets++;
78    end
79  end
80
81  SNRestdB_log = 10*log10(SNRest_log);
82  if states.Tbits
83    printf("Fs: %d Rs: %d rate %4.2f (%d,%d) frames received: %3d SNRav: %4.2f\n",
84    Fs, Rs, coderate, n, k, num_frames_rx, mean(SNRestdB_log));
85    uber = states.Terrs/states.Tbits; cber = Terrs/Tbits; cper = Tperr/Tpackets;
86    printf("  Uncoded: nbits: %6d nerrs: %6d ber: %4.3f\n", states.Tbits, states.Terrs, uber);
87    printf("  Coded..: nbits: %6d nerrs: %6d ber: %4.3f\n", Tbits, Terrs, cber);
88    printf("  Coded..: npckt: %6d perrs: %6d per: %4.3f\n", Tpackets, Tperr, cper);
89  else
90    printf("No frames detected....\n");
91  end
92
93  figure(1); clf; subplot(211); plot(rx); axis([1 length(rx) -32767 32767]); subplot(212); plot_specgram(rx);
94  figure(2); clf;
95  subplot(211); plot(f_log); axis([1 length(f_log) states.fest_fmin states.fest_fmax]); ylabel('Tone Freq (Hz)');
96  subplot(212); plot(rx_timing_log); axis([1 length(rx_timing_log) -0.5 0.5]); ylabel('Timing');
97  figure(3); clf;
98  mx_SNRestdB = 5*ceil(max(SNRestdB_log)/5);
99  subplot(211); plot(SNRestdB_log); axis([1 length(SNRestdB_log) 0 mx_SNRestdB]); ylabel('SNRest (dB)');
100  subplot(212); stem(log_nerrs); ylabel('Uncoded errors');
101end
102
103