reasonable results with HF channel model, about 1dB IL on HF at 6dB Eb/No. Next...
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Thu, 30 Mar 2017 04:40:35 +0000 (04:40 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Thu, 30 Mar 2017 04:40:35 +0000 (04:40 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@3085 01035d8c-6547-0410-b346-abe4f91aad63

codec2-dev/octave/bpsk_hf.m

index 7bf0f5ed192ad8304e88fb8eb1687b740d492c63..fa5ffcd60e12e9939b2badd1ca1a5039785017f8 100644 (file)
@@ -6,8 +6,16 @@
 
 #{
  TODO:
-   [ ] sim pilot based phase est using known symbols
-   [ ] test AWGN BER with averaging pilots from adj carriers
+   [X] sim pilot based phase est using known symbols
+   [X] test AWGN BER with averaging pilots from adj carriers
+   [X] refactor to insert pilot rows
+   [X] add border cols, not used for data
+   [X] centre est on current carrier, extend to > 3
+   [ ] test single points
+   [ ] first pass curves taking into account pilot losses
+   [ ] consider interpolation
+   [ ] consider combining mod stripping phase est techniques inside frame
+   [ ] remove border carriers, interpolate edge carrier
 #}
 
 1;
@@ -19,13 +27,44 @@ function sim_out = run_sim(sim_in)
   EbNodB = sim_in.EbNodB;
   verbose = sim_in.verbose;
   hf_en = sim_in.hf_en;
+  hf_phase = sim_in.hf_phase;
   Ns = sim_in.Ns;          % step size for pilots
   Nc = sim_in.Nc;          % Number of cols, aka number of carriers
   Nr = Nbits/Nc;           % Number of rows to get Nbits total
   phase_offset = sim_in.phase_offset;
 
+  if verbose
+    printf("Nbits: %d\n", Nbits);
+    printf("Nc: %d\n", Nc);
+    printf("Nr: %d\n", Nr);
+    printf("Ns: %d (step size for pilots, Ns-1 data symbols between pilots)\n", Ns);
+  end
+
+  % check if Nbits fit neatly into carriers
+
   assert(Nbits/Nc == floor(Nbits/Nc), "Nbits/Nc must be an integer");
+
+  % check if bits fit neatly into frames with rows of pilots above and below
+  % PPP
+  % DDD
+  % DDD
+  % PPP
   
+  Nbitsperframe = (Ns-1)*Nc;
+  printf("Nbitsperframe: %d\n", Nbitsperframe);
+  Nframes = Nbits/Nbitsperframe;
+  Nrowsperframe = Nbitsperframe/Nc;
+  printf("Nrowsperframe: %d\n", Nrowsperframe);
+
+  % check if Nbits fit neatly into frames delineated by pilots
+
+  assert(Nframes == floor(Nframes), "Nbits/Nbits/frame must be an integer");
+  printf("Nframes: %d\n", Nframes);
+
+  Nrp = Nr + Nframes + 1;  % number of rows once pilots inserted
+                           % extra row of pilots at end
+  printf("Nrp: %d (number of rows including pilots)\n", Nrp);
+
   % set up HF model
 
   if hf_en
@@ -34,13 +73,31 @@ function sim_out = run_sim(sim_in)
 
     dopplerSpreadHz = 1.0; path_delay = 1E-3*Rs;
 
-    spread1 = doppler_spread(dopplerSpreadHz, Rs, Nr*2);
-    spread2 = doppler_spread(dopplerSpreadHz, Rs, Nr*2);
+    spread1 = doppler_spread(dopplerSpreadHz, Rs, Nrp+10);
+    spread2 = doppler_spread(dopplerSpreadHz, Rs, Nrp+10);
+
+    % sometimes doppler_spread() doesn't return exactly the number of samples we need
+    assert(length(spread1) >= Nrp, "not enough doppler spreading samples");
+    assert(length(spread2) >= Nrp, "not enough doppler spreading samples");
+
     hf_gain = 1.0/sqrt(var(spread1)+var(spread2));
     % printf("nsymb: %d lspread1: %d\n", nsymb, length(spread1));
   end
   
-  % simulate for each Eb/No point
+  % construct an artificial phase countour for testing, linear across freq and time
+
+  if sim_in.phase_test
+    phase_test = ones(Nrp, Nc+2);
+    for r=1:Nrp
+      for c=1:Nc+2
+        phase_test(r,c) = -pi/2 + c*pi/(Nc+2) + r*0.01*2*pi;
+        phase_test(r,c) = phase_test(r,c) - 2*pi*floor((phase_test(r,c)+pi)/(2*pi));
+      end
+    end
+  end
+
+  % simulate for each Eb/No point ------------------------------------
 
   for nn=1:length(EbNodB)
     rand('seed',1);
@@ -48,62 +105,56 @@ function sim_out = run_sim(sim_in)
 
     EbNo = 10 .^ (EbNodB(nn)/10);
     variance = 1/(EbNo/2);
-    noise = sqrt(variance)*(0.5*randn(Nr,Nc) + j*0.5*randn(Nr,Nc));
+    noise = sqrt(variance)*(0.5*randn(Nrp,Nc+2) + j*0.5*randn(Nrp,Nc+2));
 
-    tx_bits = rand(Nr,Nc) > 0.5;
-    tx = 1 - 2*tx_bits;
+    % generate tx bits and insert pilot rows and border cols
+
+    tx_bits = []; tx_bits_np = [];
+    for f=1:Nframes
+      tx_bits = [tx_bits; ones(1,Nc+2)];
+      for r=1:Nrowsperframe
+        arowofbits = rand(1,Nc) > 0.5;
+        tx_bits = [tx_bits; [1 arowofbits 1]];
+        tx_bits_np = [tx_bits_np; arowofbits];
+      end
+    end      
+    tx_bits = [tx_bits; [1 ones(1,Nc) 1]]; % final row of pilots
+    [nr nc] = size(tx_bits);
+    assert(nr == Nrp);
+    %tx_bits
+
+    tx = 2*tx_bits - 1;
 
     rx = tx * exp(j*phase_offset);
 
+    if sim_in.phase_test
+      rx = rx .* exp(j*phase_test);
+    end
+
     if hf_en
 
       % simplified rate Rs simulation model that doesn't include
-      % ISI, just freq filtering.  We assume perfect phase estimation
-      % so it's just amplitude distortion.
+      % ISI, just freq filtering.
       
       % Note Rs carrier spacing, sample rate is Rs
 
-      hf_model = zeros(Nr,Nc); phase_est = zeros(Nr,Nc);
-      for i=1:Nr
-        for c=1:Nc
+      hf_model = zeros(Nr,Nc+2); phase_est = zeros(Nr,Nc);
+      for r=1:Nrp
+        for c=1:Nc+2
           w = 2*pi*c*Rs/Rs; 
-          hf_model(i,c) = hf_gain*(spread1(i) + exp(-j*w*path_delay)*spread2(i));
+          hf_model(r,c) = hf_gain*(spread1(r) + exp(-j*w*path_delay)*spread2(r));
         end
-
-        if sim_in.av_phase
-
-          % Simulate obtaining phase est by averaging hf channel model
-          % samples at carrier frequencies adjacent to each carrier.
-          % In case of edge carriers we sample HF channel model to one
-          % side of carrier for now, for a real world estimate we'd
-          % need to interpolate phase to edge carriers somehow
-
-          for c=1:Nc
-            phase_sum = 0;
-            for cc = c-1:c+1
-              w = 2*pi*cc*Rs/Rs; 
-              ahf_model = hf_gain*(spread1(i) + exp(-j*w*path_delay)*spread2(i));
-              phase_sum += ahf_model;
-            end
-            phase_est(i,c) = angle(phase_sum); 
-          end
-          rx(i,:) = rx(i,:) .* hf_model(i,:) .* exp(-j*phase_est(i,:));
-
+        
+        if hf_phase
+          rx(r,:) = rx(r,:) .* hf_model(r,:);
         else
-
-          % just apply amplitude fading, assume we have ideal phase est
-
-          rx(i,:) = rx(i,:) .* abs(hf_model(i,:));
+          rx(r,:) = rx(r,:) .* abs(hf_model(r,:));
         end
       end
     end
 
     rx += noise;
 
-    Nbits_np = Nbits;
-    rx_np = rx;
-    tx_bits_np = tx_bits;
-
     % pilot based phase est, we use known tx symbols as pilots ----------
 
     if sim_in.pilot_phase_est
@@ -112,14 +163,19 @@ function sim_out = run_sim(sim_in)
       % adjust phase of data symbol
       % demodulate and count errors of just data
  
-      phase_est_log = zeros(Nr,Nc);
-      for c=1:Nc
-        for r=1:Ns:Nr-Ns
-
-          % estimate phase using average of two pilots at either end of frame
-
-          %aphase_est_rect = rx(r,c)*tx(r,c)' +  rx(r+Ns,c)*tx(r+Ns,c)';
-          aphase_est_rect = sum(rx(r,:)*tx(r,:)') +  sum(rx(r+Ns,:)*tx(r+Ns,:)');
+      phase_est_log = 10*ones(Nrp,Nc+2);
+      for c=2:Nc+1
+        for r=1:Ns:Nrp-Ns
+
+          % estimate phase using average of 6 pilots in a rect 2D window centred
+          % on this carrier
+          % PPP
+          % DDD
+          % DDD
+          % PPP
+          
+          cr = c-1:c+1;
+          aphase_est_rect = sum(rx(r,cr)*tx(r,cr)') +  sum(rx(r+Ns,cr)*tx(r+Ns,cr)');
           aphase_est = angle(aphase_est_rect);
 
           % apply this phase estimate to symbols in frame, and check for error
@@ -130,41 +186,57 @@ function sim_out = run_sim(sim_in)
           end
         end
       end
+    end
 
-      % build up tx_bits and rx symbols matrices without pilots
+    % strip pilots to give us just data symbols
       
-      tx_bits_np = []; rx_np = []; Nbits_np = 0;
-      for r=1:Nr
-        if mod(r-1,Ns) != 0
-          tx_bits_np = [tx_bits_np; tx_bits(r,:)];
-          rx_np = [rx_np; rx(r,:)];
-          Nbits_np += Nc;
-        end
+    rx_np = [];
+    for r=1:Nrp
+      if mod(r-1,Ns) != 0
+        rx_np = [rx_np; rx(r,2:Nc+1)];
       end
     end
 
+   %phase_test
+   %phase_est_log
+
     % calculate BER stats as a block, after pilots extracted
 
-    rx_bits_np = real(rx_np) < 0;
+    rx_bits_np = real(rx_np) > 0;
     errors = xor(tx_bits_np, rx_bits_np);
     Nerrs = sum(sum(errors));
 
-    printf("EbNodB: %3.2f BER: %4.3f Nbits: %d Nerrs: %d\n", EbNodB(nn), Nerrs/Nbits_np, Nbits_np, Nerrs);
+    printf("EbNodB: %3.2f BER: %4.3f Nbits: %d Nerrs: %d\n", EbNodB(nn), Nerrs/Nbits, Nbits, Nerrs);
 
     if verbose
       figure(1); clf; 
       plot(rx_np,'+');
       axis([-2 2 -2 2]);
-      if sim_in.pilot_phase_est
-        figure(2); clf;
-        plot(phase_est_log,'+');
-      end
       if hf_en
         figure(2); clf; 
-        plot(abs(hf_model));
-        figure(3); clf; 
-        plot(angle(hf_model));
-        hold on; plot(phase_est,'+'); hold off;
+        plot(abs(hf_model(:,2:Nc+1)));
+        if sim_in.hf_phase
+          figure(3); clf; 
+          plot(angle(hf_model(:,2:Nc+1)));
+          if sim_in.pilot_phase_est
+            hold on; plot(phase_est_log(:,2:Nc+1),'+', 'markersize', 10); hold off;
+          end
+          axis([1 Nrp -pi pi]);
+        end
+      else
+        if sim_in.pilot_phase_est
+          figure(3); clf;
+          plot(phase_est_log,'+');
+        end
+      end
+
+      if sim_in.phase_test
+          figure(3); clf;
+          plot(phase_test(:,2:Nc+1));
+          if sim_in.pilot_phase_est
+            hold on; plot(phase_est_log(:,2:Nc+1),'+', 'markersize', 10); hold off;
+          end
+          axis([1 Nrp -pi pi]);
       end
     end
 
@@ -199,16 +271,16 @@ end
 
 
 function run_single
-  data_rows = 1000;
   sim_in.Nc = 3;
-  sim_in.Ns = 6;
-  sim_in.Nbits = (data_rows*sim_in.Ns+1)*sim_in.Nc;
-  sim_in.EbNodB = 2;
+  sim_in.Ns = 8;
+  sim_in.Nbits = 5000*sim_in.Nc*(sim_in.Ns-1);
+  sim_in.EbNodB = 6;
   sim_in.verbose = 1;
   sim_in.pilot_phase_est = 1;
-  sim_in.phase_offset = -pi;
-  sim_in.hf_en = 0;
-  sim_in.av_phase = 0;
+  sim_in.phase_offset = 0;
+  sim_in.phase_test = 0;
+  sim_in.hf_en = 1;
+  sim_in.hf_phase = 1;
 
   run_sim(sim_in);
 end