added analog signal simulation and interleaving of (224,112) code, and a basic sync...
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Tue, 16 May 2017 03:03:30 +0000 (03:03 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Tue, 16 May 2017 03:03:30 +0000 (03:03 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@3137 01035d8c-6547-0410-b346-abe4f91aad63

codec2-dev/octave/ofdm_rx.m
codec2-dev/octave/ofdm_tx.m

index ba3b1425c460b8d85e10942e0eb715cb38bd9bc8..c99b8b0cfe1a5800fdc1efba463a855f9f86f54e 100644 (file)
@@ -11,7 +11,7 @@
     [ ] way to fall out of sync
 #}
 
-function ofdm_rx(filename)
+function ofdm_rx(filename, interleave_frames = 1)
   ofdm_lib;
   ldpc;
   more off;
@@ -35,18 +35,6 @@ function ofdm_rx(filename)
   [code_param framesize rate] = ldpc_init_user(HRA_112_112, modulation, mod_order, mapping);
   assert(Nbitsperframe == code_param.code_bits_per_frame);
 
-  % fixed test frame of tx bits
-
-  rand('seed', 100);
-  tx_bits = round(rand(1,code_param.data_bits_per_frame));
-  tx_codeword = LdpcEncode(tx_bits, code_param.H_rows, code_param.P_matrix);
-
-  % init logs and BER stats
-
-  rx_bits = []; rx_np = []; timing_est_log = []; delta_t_log = []; foff_est_hz_log = [];
-  phase_est_pilot_log = [];
-  Terrs = Tbits = 0;
-
   % load real samples from file
 
   Ascale= 2E5;
@@ -54,6 +42,46 @@ function ofdm_rx(filename)
   Nsam = length(rx); Nframes = floor(Nsam/Nsamperframe);
   prx = 1;
 
+  % buffers for interleaved frames
+
+  Nsymbolsperframe = code_param.code_bits_per_frame/bps
+  Nsymbolsperinterleavedframe = interleave_frames*Nsymbolsperframe
+  rx_np = zeros(1, Nsymbolsperinterleavedframe);
+  rx_amp = zeros(1, Nsymbolsperinterleavedframe);
+
+  % OK generate tx frame for BER calcs
+
+  rand('seed', 100);
+  tx_bits = []; tx_codewords = [];
+  for f=1:interleave_frames
+    atx_bits = round(rand(1,code_param.data_bits_per_frame));
+    tx_bits = [tx_bits atx_bits];
+    acodeword = LdpcEncode(atx_bits, code_param.H_rows, code_param.P_matrix);
+    tx_codewords = [tx_codewords acodeword];
+  end
+
+  % used for rx frame sync on interleaved symbols - we demod the
+  % entire interleaved frame to raw bits
+
+  tx_symbols = [];
+  for s=1:Nsymbolsperinterleavedframe
+    tx_symbols = [tx_symbols qpsk_mod( tx_codewords(2*(s-1)+1:2*s) )];
+  end
+
+  tx_symbols = gp_interleave(tx_symbols);
+
+  tx_bits_raw = [];
+  for s=1:Nsymbolsperinterleavedframe
+    tx_bits_raw = [tx_bits_raw qpsk_demod(tx_symbols(s))];
+  end
+
+  % init logs and BER stats
+
+  rx_bits = []; rx_np_log = []; timing_est_log = []; delta_t_log = []; foff_est_hz_log = [];
+  phase_est_pilot_log = [];
+  Terrs = Tbits = Terrs_coded = Tbits_coded = Tpackets = Tpacketerrs = 0;
+  Nbitspervocframe = 28;
+
   % 'prime' rx buf to get correct coarse timing (for now)
 
   prx = 1;
@@ -61,7 +89,7 @@ function ofdm_rx(filename)
   states.rxbuf(Nrxbuf-nin+1:Nrxbuf) = rx(prx:nin);
   prx += nin;
 
-  state = 'searching';
+  state = 'searching'; frame_count = 0;
 
   % main loop ----------------------------------------------------------------
 
@@ -80,30 +108,53 @@ function ofdm_rx(filename)
     prx += states.nin;
 
     [rx_bits_raw states aphase_est_pilot_log arx_np arx_amp] = ofdm_demod(states, rxbuf_in);
-    rx_codeword = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, arx_np, min(EsNo,30), arx_amp);
-    rx_bits = rx_codeword(1:code_param.data_bits_per_frame);
+    frame_count++;
 
-    % measure errors and iterate start machine
+    % update sliding windows of rx-ed symbols and symbol amplitudes
 
-    errors = xor(tx_bits, rx_bits);
-    Nerrs = sum(errors);
-    next_state = state;
+    rx_np(1:Nsymbolsperinterleavedframe-Nsymbolsperframe) = rx_np(Nsymbolsperframe+1:Nsymbolsperinterleavedframe);
+    rx_np(Nsymbolsperinterleavedframe-Nsymbolsperframe+1:Nsymbolsperinterleavedframe) = arx_np;
+    rx_amp(1:Nsymbolsperinterleavedframe-Nsymbolsperframe) = rx_amp(Nsymbolsperframe+1:Nsymbolsperinterleavedframe);
+    rx_amp(Nsymbolsperinterleavedframe-Nsymbolsperframe+1:Nsymbolsperinterleavedframe) = arx_amp;
+
+    printf("f: %d state: %s frame_count: %d\n", f, state, frame_count);
+
+    % If looking for sync: check raw BER on frame just received
+    % against all possible positions in the interleaver frame.
 
+    % iterate state machine ------------------------------------
+
+    next_state = state;
     if strcmp(state,'searching')  
-      if Nerrs == 0
-        next_state = 'synced';
+
+      % If looking for sync: check raw BER on frame just received
+      % against all possible positions in the interleaver frame.
+
+      for ff=1:interleave_frames
+        st = (ff-1)*Nbitsperframe+1; en = st + Nbitsperframe - 1;
+        errors = xor(tx_bits_raw(st:en), rx_bits_raw);
+        Nerrs = sum(errors); 
+        printf("  ff: %d Nerrs: %d\n", ff, Nerrs);
+        if Nerrs/Nbitsperframe < 0.1
+          next_state = 'synced';
+          frame_count = ff;
+          if f < interleave_frames
+            % point trying a LDPC decode if we don't have a full frame!
+            frame_count -= interleave_frames;
+          end
+        end
       end
     end
     state = next_state;
 
     if strcmp(state,'searching') 
 
-      % attempt coarse timing estimate (i.e. detect start of frame)
+      % still searching? Attempt coarse timing estimate (i.e. detect start of frame)
 
       st = M+Ncp + Nsamperframe + 1; en = st + 2*Nsamperframe; 
       [ct_est foff_est] = coarse_sync(states, states.rxbuf(st:en), states.rate_fs_pilot_samples);
       if states.verbose
-        printf("Nerrs: %d ct_est: %4d foff_est: %3.1f\n", Nerrs, ct_est, foff_est);
+        printf("  Nerrs: %d ct_est: %4d foff_est: %3.1f\n", Nerrs, ct_est, foff_est);
       end
 
       % calculate number of samples we need on next buffer to get into sync
@@ -114,29 +165,75 @@ function ofdm_rx(filename)
 
       states.sample_point = states.timing_est = 1;
       states.foff_est_hz = foff_est;
+    end
+    
+    if strcmp(state,'synced') && (frame_count == interleave_frames)
+
+      % de-interleave QPSK symbols
 
-    else
+      arx_np = gp_deinterleave(rx_np);
+      arx_amp = gp_deinterleave(rx_amp);
 
-      % we are in sync so log states and bit errors
+      % LDPC decode
 
-      rx_np = [rx_np arx_np];
+      rx_bits = [];
+      for ff=1:interleave_frames
+        st = (ff-1)*Nbitsperframe/bps+1; en = st + Nbitsperframe/bps - 1;
+        rx_codeword = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, arx_np(st:en), min(EsNo,30), arx_amp(st:en));
+        rx_bits = [rx_bits rx_codeword(1:code_param.data_bits_per_frame)];
+      end
+
+      errors_coded = xor(tx_bits, rx_bits);
+      Nerrs_coded = sum(errors_coded);
+      Terrs_coded += Nerrs_coded;
+      Tbits_coded += code_param.data_bits_per_frame*interleave_frames;
+      Nerrs_coded_log(f) = Nerrs_coded;
+
+      printf("  Nerrs_coded: %d\n", Nerrs_coded);
+
+      % we are in sync so log states and bit/packet error stats
+
+      rx_np_log = [rx_np_log arx_np];
       timing_est_log = [timing_est_log states.timing_est];
       delta_t_log = [delta_t_log states.delta_t];
       foff_est_hz_log = [foff_est_hz_log states.foff_est_hz];
       phase_est_pilot_log = [phase_est_pilot_log; aphase_est_pilot_log];
 
-      % measure bit errors
+      % measure uncoded bit errors
 
+      rx_bits_raw = [];
+      for s=1:Nsymbolsperinterleavedframe
+        rx_bits_raw = [rx_bits_raw qpsk_demod(rx_np(s))];
+      end
+      errors = xor(tx_bits_raw, rx_bits_raw);
+      Nerrs = sum(errors);
       Terrs += Nerrs;
       Nerrs_log(f) = Nerrs;
-      Tbits += code_param.data_bits_per_frame;
+      Tbits += code_param.code_bits_per_frame*interleave_frames;
+
+      % measure packet errors based on Codec 2 packet size
+
+      Nvocframes = floor(code_param.data_bits_per_frame*interleave_frames/Nbitspervocframe);
+      for fv=1:Nvocframes
+        st = (fv-1)*Nbitspervocframe + 1;
+        en = st + Nbitspervocframe - 1;
+        Nvocpacketerrs = sum(xor(tx_bits(st:en), rx_bits(st:en)));
+        if Nvocpacketerrs
+          Tpacketerrs++;
+        end
+        Tpackets++;
+      end
+
+      frame_count = 0;
     end
   end
 
-  printf("Coded BER: %5.4f Tbits: %d Terrs: %d\n", Terrs/Tbits, Tbits, Terrs);
+  printf("Coded BER: %5.4f Tbits: %5d Terrs: %5d PER: %5.4f Tpacketerrs: %5d Tpackets: %5d\n", 
+         Terrs_coded/Tbits_coded, Tbits_coded, Terrs_coded, Tpacketerrs/Tpackets, Tpacketerrs, Tpackets);
+  printf("Raw BER..: %5.4f Tbits: %5d Terrs: %5d\n", Terrs/Tbits, Tbits, Terrs);
 
   figure(1); clf; 
-  plot(rx_np,'+');
+  plot(rx_np_log,'+');
   axis([-2 2 -2 2]);
   title('Scatter');
 
@@ -155,11 +252,14 @@ function ofdm_rx(filename)
 
   figure(4); clf;
   plot(foff_est_hz_log)
-  mx = max(abs(foff_est_hz_log));
-  axis([1 max(Nframes,2) -mx mx]);
+  axis([1 max(Nframes,2) -2 2]);
   title('Fine Freq');
+  ylabel('Hz')
 
   figure(5); clf;
-  %plot(Nerrs_log);
-  
+  subplot(211)
+  title('Nerrs Log')
+  stem(Nerrs_log);
+  subplot(212)
+  stem(Nerrs_coded_log);
 endfunction
index a414fc4af81ceeac5a5a8e3add61f5f340055b1f..d706a491a745594e798eec8a200c78abb6df7b94 100644 (file)
  
 #}
 
-function ofdm_tx(filename, Nsec, EbNodB=100, channel='awgn', freq_offset_Hz=0)
+function ofdm_tx(filename, Nsec, interleave_frames = 1, EbNodB=100, channel='awgn', freq_offset_Hz=0)
   ofdm_lib;
   ldpc;
+  gp_interleaver;
 
   % init modem
 
@@ -38,20 +39,56 @@ function ofdm_tx(filename, Nsec, EbNodB=100, channel='awgn', freq_offset_Hz=0)
   [code_param framesize rate] = ldpc_init_user(HRA_112_112, modulation, mod_order, mapping);
   assert(Nbitsperframe == code_param.code_bits_per_frame);
 
-  % generate fixed test frame of tx bits and run OFDM modulator
-  % todo: maybe extend this to 4 or 8 frames, one is a bit short
+  % Generate fixed test frame of tx bits and run OFDM modulator
 
   Nrows = Nsec*Rs;
   Nframes = floor((Nrows-1)/Ns);
 
-  rand('seed', 100);
-  tx_bits = round(rand(1,code_param.data_bits_per_frame));
-  tx_codeword = LdpcEncode(tx_bits, code_param.H_rows, code_param.P_matrix);
+  % Adjust Nframes so we have an integer number of interleaver frames
+  % in simulation
+
+  Nframes = interleave_frames*round(Nframes/interleave_frames);
 
+  % OK generate an interleaver frame of codewords from random data bits
+
+  rand('seed', 100);
+  tx_bits = tx_symbols = [];
+  for f=1:interleave_frames
+    atx_bits = round(rand(1,code_param.data_bits_per_frame));
+    tx_bits = [tx_bits atx_bits];
+    codeword = LdpcEncode(atx_bits, code_param.H_rows, code_param.P_matrix);
+    for b=1:2:Nbitsperframe
+      tx_symbols = [tx_symbols qpsk_mod(codeword(b:b+1))];
+    end
+  end
+  tx_symbols = gp_interleave(tx_symbols);
+  
+  atx = [];
+  for f=1:interleave_frames
+    st = (f-1)*Nbitsperframe/bps+1; en = st + Nbitsperframe/bps-1;
+    atx = [atx ofdm_txframe(states, tx_symbols(st:en))];
+  end
   tx = [];
-  for f=1:Nframes
-    tx = [tx ofdm_mod(states, tx_codeword)];
+  for f=1:Nframes/interleave_frames
+    tx = [tx atx];
   end
+
+  % not very pretty way to process analog signals with exactly the same channel
+
+  analog_hack = 0;
+  if analog_hack
+    s = load_raw('../raw/ve9qrp_10s.raw')';
+    tx = hilbert(s);
+
+    % normalise power to same as ofdm tx
+
+    nom_tx_pwr = 2/(Ns*(M*M)) + Nc/(M*M);
+    tx_pwr = var(tx);
+    tx *= sqrt(nom_tx_pwr/tx_pwr);
+  end
+
   Nsam = length(tx);
 
   % channel simulation
@@ -61,7 +98,7 @@ function ofdm_tx(filename, Nsec, EbNodB=100, channel='awgn', freq_offset_Hz=0)
   woffset = 2*pi*freq_offset_Hz/Fs;
 
   SNRdB = EbNodB + 10*log10(700/3000);
-  printf("EbNo: %3.1f dB  SNR(3k) est: %3.1f dB  foff: %3.1fHz\n", EbNodB, SNRdB, freq_offset_Hz);
+  printf("EbNo: %3.1f dB  SNR(3k) est: %3.1f dB  foff: %3.1fHz ", EbNodB, SNRdB, freq_offset_Hz);
 
   % set up HF model ---------------------------------------------------------------
 
@@ -108,7 +145,7 @@ function ofdm_tx(filename, Nsec, EbNodB=100, channel='awgn', freq_offset_Hz=0)
 
   noise = sqrt(variance/2)*0.5*randn(1,Nsam);
   rx = real(rx) + noise;
-  printf("measured SNR: %3.2f dB\n", 10*log10(var(real(tx))/var(noise)));
+  printf("measured SNR: %3.2f dB\n", 10*log10(var(real(tx))/var(noise))+10*log10(4000) - 10*log10(3000));
 
   Ascale = 2E5;
   frx=fopen(filename,"wb"); fwrite(frx, Ascale*rx, "short"); fclose(frx);