1% ofdm_ldpc_rx.m
2% David Rowe April 2017
3%
4% OFDM file based rx, with LDPC and interleaver, Octave version of src/ofdm_demod.c
5
6#{
7    1. Streaming mode operation:
8
9       ofdm_ldpc_rx("test_700d.raw","700D")
10
11    2. Burst mode, tell state machine there is one packet in each burst:
12
13       ofdm_ldpc_rx("test_datac0.raw","datac0","packetsperburst",1)
14
15#}
16
17function ofdm_ldpc_rx(filename, mode="700D", varargin)
18  ofdm_lib;
19  ldpc;
20  gp_interleaver;
21  more off;
22  pkg load signal;
23
24  % init modem
25
26  config = ofdm_init_mode(mode);
27  states = ofdm_init(config);
28  ofdm_load_const;
29  states.verbose = 1;
30  pass_packet_count = 0;
31
32  i=1;
33  while i <= length(varargin)
34    if strcmp(varargin{i},"packetsperburst")
35      states.data_mode = "burst"; % use pre/post amble based sync
36      states.packetsperburst = varargin{i+1}; i++;
37      states.postambledetectoren = 1;
38    elseif strcmp(varargin{i},"passpacketcount")
39      pass_packet_count = varargin{i+1}; i++;
40    else
41      printf("\nERROR unknown argument: [%d] %s \n", i ,varargin{i});
42      return;
43    end
44    i++;
45  end
46
47  % some constants used for assembling modem frames
48
49  [code_param Nbitspercodecframe Ncodecframespermodemframe] = codec_to_frame_packing(states, mode);
50
51  % load real samples from file
52
53  Ascale= states.amp_scale/2.0;  % /2 as real signal has half amplitude
54  frx=fopen(filename,"rb"); rx = fread(frx, Inf, "short")/Ascale; fclose(frx);
55  Nsam = length(rx);
56  prx = 1;
57
58  % Generate tx frame for BER calcs
59
60  payload_bits = round(ofdm_rand(code_param.data_bits_per_frame)/32767);
61  tx_bits = fec_encode(states, code_param, mode, payload_bits, Ncodecframespermodemframe, Nbitspercodecframe);
62
63  % Some handy constants
64
65  Nsymsperframe = Nbitsperframe/bps;
66  Nsymsperpacket = Nbitsperpacket/bps;
67  Ncodedbitsperpacket = code_param.coded_bits_per_frame;
68  Ncodedsymsperpacket = code_param.coded_syms_per_frame;
69
70  % init logs and BER stats
71
72  rx_bits = []; rx_np_log = []; timing_est_log = []; delta_t_log = []; foff_est_hz_log = [];
73  channel_est_pilot_log = []; snr_log = []; mean_amp_log = [];
74  Terrs = Tbits = Terrs_coded = Tbits_coded = Perrs_coded = 0;
75  Nerrs_coded_log = Nerrs_log = [];
76  error_positions = [];
77  Nerrs_coded = Nerrs_raw = 0;
78  paritychecks = [0];
79  EsNo = 1;
80  rx_uw = zeros(1,states.Nuwbits);
81
82  rx_syms = zeros(1,Nsymsperpacket); rx_amps = zeros(1,Nsymsperpacket);
83  packet_count = frame_count = 0;
84
85  % main loop ----------------------------------------------------------------
86
87  f = 1;
88  while(prx < Nsam)
89
90    % insert samples at end of buffer, set to zero if no samples
91    % available to disable phase estimation on future pilots on last
92    % frame of simulation
93
94    lnew = min(Nsam-prx,states.nin);
95    rxbuf_in = zeros(1,states.nin);
96
97    if lnew
98      rxbuf_in(1:lnew) = rx(prx:prx+lnew-1);
99    end
100    prx += states.nin;
101
102    if states.verbose
103      printf("f: %3d nin: %4d st: %-6s ", f, states.nin, states.sync_state);
104    end
105
106    if strcmp(states.sync_state,'search')
107      [timing_valid states] = ofdm_sync_search(states, rxbuf_in);
108    else
109      % accumulate a buffer of data symbols for this packet
110      rx_syms(1:end-Nsymsperframe) = rx_syms(Nsymsperframe+1:end);
111      rx_amps(1:end-Nsymsperframe) = rx_amps(Nsymsperframe+1:end);
112      [states rx_bits achannel_est_pilot_log arx_np arx_amp] = ofdm_demod(states, rxbuf_in);
113      rx_syms(end-Nsymsperframe+1:end) = arx_np;
114      rx_amps(end-Nsymsperframe+1:end) = arx_amp;
115
116      rx_uw = extract_uw(states, rx_syms(end-Nuwframes*Nsymsperframe+1:end), rx_amps(end-Nuwframes*Nsymsperframe+1:end));
117
118      % We need the full packet of symbols before disassembling and checking for bit errors
119      if (states.modem_frame == (states.Np-1))
120        packet_count++;
121
122        % unpack, de-interleave PSK symbols and symbol amplitudes
123        [rx_uw_unused payload_syms payload_amps txt_bits] = disassemble_modem_packet(states, rx_syms, rx_amps);
124        payload_syms_de = gp_deinterleave(payload_syms);
125        payload_amps_de = gp_deinterleave(payload_amps);
126
127        % Count uncoded (raw) errors
128        rx_bits = zeros(1,Ncodedbitsperpacket);
129        for s=1:Ncodedsymsperpacket
130          if bps == 2 rx_bits(2*s-1:2*s) = qpsk_demod(payload_syms_de(s)); end
131          if bps == 4 rx_bits(bps*(s-1)+1:bps*s) = qam16_demod(states.qam16,payload_syms_de(s), payload_amps_de(s)); end
132        end
133        errors = xor(tx_bits, rx_bits);
134        Nerrs = sum(errors);
135        Nerrs_log = [Nerrs_log Nerrs]; Nerrs_raw = Nerrs;
136        Terrs += Nerrs;
137        Tbits += Nbitsperpacket;
138
139        % LDPC decode
140
141        % keep earlier mean amplitude estimator for compatability with 700D
142        if states.amp_est_mode == 0
143          mean_amp = states.mean_amp;
144        else
145          mean_amp = mean(payload_amps_de)+1E-12;
146        end
147        mean_amp_log = [mean_amp_log mean_amp];
148
149        % used fixed EsNo est, as EsNo estimator for QAM not working very well at this stage
150        EsNo = 10^(states.EsNodB/10);
151
152        % TODO 2020 support for padding with known data bits
153
154        [rx_codeword paritychecks] = ldpc_dec(code_param, mx_iter=100, demod=0, dec=0, ...
155                                              payload_syms_de/mean_amp, EsNo, payload_amps_de/mean_amp);
156        rx_bits = rx_codeword(1:code_param.data_bits_per_frame);
157        errors = xor(payload_bits, rx_bits);
158        Nerrs_coded  = sum(errors);
159
160        if Nerrs_coded Perrs_coded++; end
161        Terrs_coded += Nerrs_coded;
162        Tbits_coded += code_param.data_bits_per_frame;
163        Nerrs_coded_log = [Nerrs_coded_log Nerrs_coded];
164
165        % per-packet SNR estimate
166        EsNo_estdB = esno_est_calc(rx_syms);
167        SNR_estdB = snr_from_esno(states, EsNo_estdB);
168        snr_log = [snr_log SNR_estdB];
169      end
170
171      % we are in sync so log modem states
172
173      rx_np_log = [rx_np_log arx_np];
174      timing_est_log = [timing_est_log states.timing_est];
175      delta_t_log = [delta_t_log states.delta_t];
176      foff_est_hz_log = [foff_est_hz_log states.foff_est_hz];
177      channel_est_pilot_log = [channel_est_pilot_log; achannel_est_pilot_log];
178      frame_count++;
179    end
180
181    states = sync_state_machine(states, rx_uw);
182
183    if states.verbose
184      if strcmp(states.last_sync_state,'search') == 0
185        pcc = max(paritychecks);
186        iter = 0;
187        for i=1:length(paritychecks)
188          if paritychecks(i) iter=i; end
189        end
190        % complete logging line
191        if (states.modem_frame == 0) && (strcmp(states.last_sync_state, "trial") == 0)
192            printf("euw: %3d %d mf: %2d pbw: %s foff: %4.1f eraw: %3d ecod: %3d iter: %3d pcc: %3d snr: %5.2f",
193                    states.uw_errors, states.sync_counter, states.modem_frame, states.phase_est_bandwidth(1), states.foff_est_hz,
194                    Nerrs_raw, Nerrs_coded, iter, pcc, SNR_estdB);
195        else
196            printf("euw: %3d %d mf: %2d pbw: %s foff: %4.1f",
197                    states.uw_errors, states.sync_counter, states.modem_frame, states.phase_est_bandwidth(1), states.foff_est_hz);
198        end
199      end
200      printf("\n");
201    end
202
203    % reset stats if in streaming mode, don't reset if in burst mode
204    if strcmp(states.data_mode, "streaming") && states.sync_start
205      Nerrs_raw = Nerrs_coded = 0;
206      Nerrs_log = [];
207      Terrs = Tbits = 0;
208      Tpacketerrs = Tpackets = 0;
209      Terrs_coded = Tbits_coded = 0;
210      error_positions = Nerrs_coded_log = [];
211    end
212    f++;
213  end
214  Nframes = f;
215
216  printf("Raw BER..: %5.4f Tbits: %5d Terrs: %5d SNR3k: %5.2f\n", Terrs/(Tbits+1E-12), Tbits, Terrs, mean(snr_log));
217  printf("Coded BER: %5.4f Tbits: %5d Terrs: %5d\n", Terrs_coded/(Tbits_coded+1E-12), Tbits_coded, Terrs_coded);
218  printf("Coded PER: %5.4f Pckts: %5d Perrs: %5d Npre: %d Npost: %d\n",
219         Perrs_coded/(packet_count+1E-12), packet_count, Perrs_coded,  states.npre, states.npost);
220
221  if length(rx_np_log)
222      figure(1); clf;
223      plot(exp(j*pi/4)*rx_np_log(floor(end/4):floor(end-end/8)),'+');
224      mx = 2*mean(abs(channel_est_pilot_log(:)));
225      axis([-mx mx -mx mx]);
226      title('Scatter');
227
228      figure(2); clf;
229      plot(angle(channel_est_pilot_log),'g+', 'markersize', 5);
230      title('Phase est');
231      axis([1 length(channel_est_pilot_log) -pi pi]);
232
233      figure(3); clf;
234      amp_est = abs(channel_est_pilot_log);
235      plot(amp_est,'g+', 'markersize', 5);
236      title('Amp est');
237      axis([1 length(channel_est_pilot_log) min(amp_est(:)) max(amp_est(:))]);
238
239      figure(4); clf;
240      subplot(211); plot(snr_log); ylabel('SNR3kdB');
241      subplot(212); plot(mean_amp_log); ylabel('mean amp');
242
243      figure(5); clf;
244      subplot(211)
245      stem(delta_t_log)
246      title('delta t');
247      subplot(212)
248      plot(timing_est_log);
249      title('timing est');
250
251      figure(6); clf;
252      plot(foff_est_hz_log)
253      mx = max(max(abs(foff_est_hz_log)),1);
254      axis([1 max(Nframes,2) -mx mx]);
255      title('Fine Freq');
256      ylabel('Hz')
257  end
258
259  if length(Nerrs_log) > 1
260    figure(7); clf;
261    subplot(211)
262    stem(Nerrs_log);
263    title('Uncoded errrors/modem frame')
264    axis([1 length(Nerrs_log) 0 Nbitsperpacket*0.2]);
265    if length(Nerrs_coded_log)
266      subplot(212)
267      stem(Nerrs_coded_log);
268      title('Coded errors/mode frame')
269      axis([1 length(Nerrs_coded_log) 0 Nbitsperpacket*0.2]);
270    end
271  end
272
273  figure(9); clf; plot_specgram(rx);
274
275  if pass_packet_count > 0
276    if packet_count >= pass_packet_count printf("Pass!\n"); else printf("Fail!\n"); end;
277  end
278endfunction
279