discovered BER measurement using pilots can't work as pilots get rotated away from...

git-svn-id: https://svn.code.sf.net/p/freetel/code@3039 01035d8c-6547-0410-b346-abe4f91aad63

diff --git a/codec2-dev/octave/cohpsk.m b/codec2-dev/octave/cohpsk.m
index 4aeb73339cff1f48292b0336b66def178fe756b9..34abfa7d6e4d5ba7d200c87da95236a63401d7e5 100644 (file)
--- a/codec2-dev/octave/cohpsk.m
+++ b/codec2-dev/octave/cohpsk.m
@@ -175,6 +175,8 @@ endfunction
 
 
 function [tx_symb tx_bits] = bits_to_qpsk_symbols(sim_in, tx_bits, code_param)
+    ldpc_code     = sim_in.ldpc_code;
+    rate          = sim_in.ldpc_code_rate;
     framesize     = sim_in.framesize;
     Nsymbrow      = sim_in.Nsymbrow;
     Nsymbrowpilot = sim_in.Nsymbrowpilot;
@@ -185,6 +187,10 @@ function [tx_symb tx_bits] = bits_to_qpsk_symbols(sim_in, tx_bits, code_param)
     pilot         = sim_in.pilot;
     Nd            = sim_in.Nd;
 
+    if ldpc_code
+        [tx_bits, tmp] = ldpc_enc(tx_bits, code_param);
+    end
+
     % modulate --------------------------------------------
 
     % organise symbols into a Nsymbrow rows by Nc cols
@@ -241,7 +247,7 @@ function [rx_symb rx_bits rx_symb_linear amp_ phi_ sig_rms noise_rms cohpsk] = q
 
     sampling_points = [1 2 cohpsk.Nsymbrow+3 cohpsk.Nsymbrow+4];
     pilot2 = [cohpsk.pilot(1,:); cohpsk.pilot(2,:); cohpsk.pilot(1,:); cohpsk.pilot(2,:);];
-    phi_ = zeros(Nsymbrowpilot, Nc*Nd);
+    phi_ = zeros(Nsymbrow, Nc*Nd);
     amp_ = zeros(Nsymbrow, Nc*Nd);
     
     for c=1:Nc*Nd
@@ -250,7 +256,7 @@ function [rx_symb rx_bits rx_symb_linear amp_ phi_ sig_rms noise_rms cohpsk] = q
      
       y = ct_symb_buf(sampling_points,c) .* pilot2(:,c-Nc*floor((c-1)/Nc));
       [m b] = linreg(sampling_points, y, length(sampling_points));
-      yfit = m*[1 2 3 4 5 6] + b;
+      yfit = m*[3 4 5 6] + b;
       phi_(:, c) = angle(yfit);
       %for r=1:Nsymbrow
       %  printf("  %f", phi_(r,c));
@@ -266,49 +272,31 @@ function [rx_symb rx_bits rx_symb_linear amp_ phi_ sig_rms noise_rms cohpsk] = q
     rx_symb_linear = zeros(1, Nsymbrow*Nc*Nd);
     rx_bits = zeros(1, framesize);
     for c=1:Nc*Nd
-      for r=1:Nsymbrowpilot
+      for r=1:Nsymbrow
         if coh_en == 1
-          rx_symb(r,c) = ct_symb_buf(r,c)*exp(-j*phi_(r,c));
+          rx_symb(r,c) = ct_symb_buf(2+r,c)*exp(-j*phi_(r,c));
         else
-          rx_symb(r,c) = ct_symb_buf(r,c);
-        end
-        i = (c-1)*Nsymbrow + r - Npilotsframe;
-        if r > Npilotsframe
-          rx_symb_linear(i) = rx_symb(r,c);
+          rx_symb(r,c) = ct_symb_buf(2+r,c);
         end
+        i = (c-1)*Nsymbrow + r;
+        rx_symb_linear(i) = rx_symb(r,c);
         %printf("phi_ %d %d %f %f\n", r,c,real(exp(-j*phi_(r,c))), imag(exp(-j*phi_(r,c))));
       end
     end
 
-    % and finally optional diversity combination and make decn on data and pilot bits
+    % and finally optional diversity combination and make decn on bits
 
-    tx_pilot_bits = rx_pilot_bits = zeros(Npilotsframe*Nc,1);
     for c=1:Nc
-      for r=1:Nsymbrowpilot
+      for r=1:Nsymbrow
+        i = (c-1)*Nsymbrow + r;
         div_symb = rx_symb(r,c);
         for d=1:Nd-1
           div_symb += rx_symb(r,c + Nc*d);
         end
-        if r > Npilotsframe
-          % data symbols, ouput of demod
-          i = (c-1)*Nsymbrow + r - Npilotsframe;
-          rx_bits((2*(i-1)+1):(2*i)) = qpsk_demod(div_symb);
-        else 
-          % demodulated pilot symbols, which we can use to estimate BER
-          %printf("%d %d (%f %f) (%f %f)\n", r, c, real(div_symb), imag(div_symb), real(cohpsk.pilot(r,c)), imag(cohpsk.pilot(r,c)));
-          i = (c-1)*Npilotsframe + r;
-          tx_pilot_bits((2*(i-1)+1):(2*i)) = qpsk_demod(cohpsk.pilot(r,c));
-          rx_pilot_bits((2*(i-1)+1):(2*i)) = qpsk_demod(div_symb);
-        end
+        rx_bits((2*(i-1)+1):(2*i)) = qpsk_demod(div_symb);
       end
     end
 
-    % estimate BER from pilot bits
-
-    cohpsk.npilotbits += Npilotsframe*Nc;
-    nerr = sum(xor(tx_pilot_bits, rx_pilot_bits));
-    cohpsk.npilotbiterrors += nerr;
-
     % Estimate noise power from demodulated symbols.  One method is to
     % calculate the distance of each symbol from the average symbol
     % position. However this is complicated by fading, which means the
@@ -744,18 +732,17 @@ function sim_out = ber_test(sim_in)
 
         ct_symb_buf = zeros(2*Nsymbrowpilot, Nc*Nd);
         prev_tx_symb = prev_rx_symb = ones(1, Nc*Nd);
-      
+
         % simulation starts here-----------------------------------
  
         for nn = 1:Ntrials+2
-            printf("%d %d\n", nn, sim_in.framesize);
-                 
+                  
             if ldpc_code
               tx_bits = round(rand(1,framesize*rate));                       
             else
               tx_bits = round(rand(1,framesize));                       
             end
-            
+
             if strcmp(modulation,'qpsk')
 
               [tx_symb tx_bits] = bits_to_qpsk_symbols(sim_in, tx_bits, code_param);
@@ -771,7 +758,6 @@ function sim_out = ber_test(sim_in)
               tx_bits_buf(1:framesize) = tx_bits;
             end
 
-            printf("b\n");
             s_ch = tx_symb;
 
             % HF channel simulation  ------------------------------------
@@ -808,8 +794,6 @@ function sim_out = ber_test(sim_in)
                 end
             end
            
-            printf("c\n");
-
             % keep a record of each tx symbol so we can check average power
 
             for r=1:Nsymbrow
@@ -834,7 +818,7 @@ function sim_out = ber_test(sim_in)
             ct_symb_buf(Nsymbrowpilot+1:2*Nsymbrowpilot,:) = s_ch;
 
             if strcmp(modulation,'qpsk')
-              [rx_symb rx_bits rx_symb_linear amp_ phi_] = qpsk_symbols_to_bits(sim_in, ct_symb_buf(1:Nsymbrowpilot+Npilotsframe,:));                                 
+              [rx_symb rx_bits rx_symb_linear amp_ phi_ EsNo_ sim_in] = qpsk_symbols_to_bits(sim_in, ct_symb_buf(1:Nsymbrowpilot+Npilotsframe,:));                                 
               phi_log = [phi_log; phi_];
               amp_log = [amp_log; amp_];
             end
@@ -842,7 +826,6 @@ function sim_out = ber_test(sim_in)
               [rx_symb rx_bits rx_symb_linear prev_rx_symb] = dqpsk_symbols_to_bits(sim_in, s_ch, prev_rx_symb);                                 
             end
                         
-            printf("d\n");
             % Wait until we have enough frames to do pilot assisted phase estimation
 
             if nn > 1
@@ -874,15 +857,12 @@ function sim_out = ber_test(sim_in)
                 Tbitsldpc += framesize*rate;
               end
             end
-
-            printf("e\n");
-          end  % 1:Ntrials
- 
-          printf("f\n");
+          end
+           
           TERvec(ne) = Terrs;
           BERvec(ne) = Terrs/Tbits;
 
-           if verbose 
+            if verbose 
               av_tx_pwr = (s_ch_tx_log * s_ch_tx_log')/length(s_ch_tx_log);
 
               printf("EsNo (dB): %3.1f Terrs: %d Tbits: %d BER %5.3f QPSK BER theory %5.3f av_tx_pwr: %3.2f",
@@ -894,7 +874,6 @@ function sim_out = ber_test(sim_in)
               end
               printf("\n");
             end
-
     end
     
     Ebvec = Esvec - 10*log10(bps);

diff --git a/codec2-dev/octave/tcohpsk.m b/codec2-dev/octave/tcohpsk.m
index 0674af2f953e8cd70ad36b831efbaf165fbc624a..59ddae01d42157734b700da1e741f84db403a62b 100644 (file)
--- a/codec2-dev/octave/tcohpsk.m
+++ b/codec2-dev/octave/tcohpsk.m
@@ -1,7 +1,7 @@
 % tcohpsk.m
 % David Rowe Oct 2014
 %
-% Octave coherent PSK modem script that has two modes:
+% Octave coherent PSK modem script that hs two modes:
 %
 % i) tests the C port of the coherent PSK modem.  This script loads
 %    the output of unittest/tcohpsk.c and compares it to the output of
@@ -64,8 +64,8 @@ randn('state',1);
 
 % select which test  ----------------------------------------------------------
 
-%test = 'compare to c';
-test = 'awgn';
+test = 'compare to c';
+%test = 'awgn';
 %test = 'fading';
 
 % some parameters that can be over ridden, e.g. to disable parts of modem
@@ -153,7 +153,7 @@ Fcentre = afdmdv.Fcentre = 1500;
 afdmdv.Fsep = afdmdv.Rs*(1+excess_bw);
 afdmdv.phase_tx = ones(afdmdv.Nc+1,1);
 % non linear carrier spacing, combined with clip, helps PAPR a lot!
-%freq_hz = afdmdv.Fsep*( -Nc*Nd/2 - 0.5 + (1:Nc*Nd).^0.98 )
+freq_hz = afdmdv.Fsep*( -Nc*Nd/2 - 0.5 + (1:Nc*Nd).^0.98 )
 afdmdv.freq_pol = 2*pi*freq_hz/Fs;
 afdmdv.freq = exp(j*afdmdv.freq_pol);
 afdmdv.Fcentre = 1500;
@@ -199,8 +199,6 @@ acohpsk.do_write_pilot_file = 1;      % enable this to dump pilot symbols to C .
 acohpsk = symbol_rate_init(acohpsk);
 acohpsk.Ndft = 1024;
 acohpsk.f_est = afdmdv.Fcentre;
-acohpsk.npilotbits       = 0;
-acohpsk.npilotbiterrors  = 0;
 
 ch_fdm_frame_buf = zeros(1, Nsw*acohpsk.Nsymbrowpilot*afdmdv.M);
 
@@ -492,7 +490,7 @@ for f=1:frames;
   % if we are in sync complete demodulation with symbol rate processing
 
   if (next_sync == 1) || (sync == 1)
-    [rx_symb rx_bits rx_symb_linear amp_ phi_ sig_rms noise_rms acohpsk] = qpsk_symbols_to_bits(acohpsk, acohpsk.ct_symb_ff_buf);
+    [rx_symb rx_bits rx_symb_linear amp_ phi_ sig_rms noise_rms] = qpsk_symbols_to_bits(acohpsk, acohpsk.ct_symb_ff_buf);
     rx_symb_log = [rx_symb_log; rx_symb];
     rx_amp_log = [rx_amp_log; amp_];
     rx_phi_log = [rx_phi_log; phi_];
@@ -547,9 +545,7 @@ for f=1:frames;
 end
 
 ber = Nerrs/Tbits;
-printf("\nOctave EsNodB: %4.1f BER......: %4.3f Nerrs..: %d Tbits..: %d\n", EsNodB, ber, Nerrs, Tbits);
-pilot_ber = acohpsk.npilotbiterrors/acohpsk.npilotbits;
-printf("                    Pilot BER: %4.3f Nerrs..: %d Tbits..: %d\n", pilot_ber, acohpsk.npilotbiterrors, acohpsk.npilotbits);
+printf("\nOctave EsNodB: %4.1f ber..: %4.3f Nerrs..: %d Tbits..: %d\n", EsNodB, ber, Nerrs, Tbits);
 
 if compare_with_c