octave modem with freq offset est working OK, good BER, time to re-port C code
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Sat, 28 Mar 2015 22:43:19 +0000 (22:43 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Sat, 28 Mar 2015 22:43:19 +0000 (22:43 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@2094 01035d8c-6547-0410-b346-abe4f91aad63

codec2-dev/octave/cohpsk.m
codec2-dev/octave/tcohpsk.m

index a8e990c9908976c49dd4fd12b93072451e2a005a..f7978dd50947688b38f2d5d80624ed1f1a8a4062 100644 (file)
@@ -93,6 +93,8 @@ function sim_in = symbol_rate_init(sim_in)
     sim_in.Nct_sym_buf = 2*Nsymbrowpilot + 2;
     sim_in.ct_symb_buf = zeros(sim_in.Nct_sym_buf, Nc);
 
+    sim_in.ff_phase = 1;
+
     % Init LDPC --------------------------------------------------------------------
 
     if ldpc_code
@@ -189,121 +191,42 @@ end
 
 % Symbol rate processing for rx side (demodulator) -------------------------------------------------------
 
-function [rx_symb rx_bits rx_symb_linear amp_linear amp_ phi_ EsNo_ prev_sym_rx sim_in] = qpsk_symbols_to_bits(sim_in, s_ch, prev_sym_rx)
-    framesize     = sim_in.framesize;
-    Nsymb         = sim_in.Nsymb;
-    Nsymbrow      = sim_in.Nsymbrow;
-    Nsymbrowpilot = sim_in.Nsymbrowpilot;
-    Nc            = sim_in.Nc;
-    Npilotsframe  = sim_in.Npilotsframe;
-    Ns            = sim_in.Ns;
-    Np            = sim_in.Np;
-    Nchip         = sim_in.Nchip;
-    modulation    = sim_in.modulation;
-    pilot         = sim_in.pilot;
-    rx_symb_buf   = sim_in.rx_symb_buf;
-    rx_pilot_buf  = sim_in.rx_pilot_buf;
-    tx_pilot_buf  = sim_in.tx_pilot_buf;
-    verbose       = sim_in.verbose;
-
-    % demodulate stage 1
-
-    for r=1:Nsymbrowpilot
-      for c=1:Nc*Nchip
-        rx_symb(r,c) = s_ch(r, c);
-        if strcmp(modulation,'dqpsk')
-          tmp = rx_symb(r,c);
-          rx_symb(r,c) *= conj(prev_sym_rx(c)/abs(prev_sym_rx(c)));
-          prev_sym_rx(c) = tmp;
-        end
-      end
-    end
-           
-    % strip out pilots
-
-    rx_symb_pilot = rx_symb;
-    rx_symb = zeros(Nsymbrow, Nc*Nchip);
-    rx_pilot = zeros(Npilotsframe, Nc*Nchip);
+function [rx_symb rx_bits amp_ phi_ EsNo_ cohpsk] = qpsk_symbols_to_bits(cohpsk, ct_symb_buf)
+    framesize     = cohpsk.framesize;
+    Nsymb         = cohpsk.Nsymb;
+    Nsymbrow      = cohpsk.Nsymbrow;
+    Nsymbrowpilot = cohpsk.Nsymbrowpilot;
+    Nc            = cohpsk.Nc;
+    Npilotsframe  = cohpsk.Npilotsframe;
+    pilot         = cohpsk.pilot;
+    verbose       = cohpsk.verbose;
+
+    % average pilots to get phase and amplitude estimates
+    % we assume there are two samples at the start of each frame and two at the end
+
+    sampling_points = [1 2 7 8];
+    pilot2 = [cohpsk.pilot(1,:); cohpsk.pilot(2,:); cohpsk.pilot(1,:); cohpsk.pilot(2,:);];
+    phi_ = zeros(Nsymbrow, Nc);
+    amp_ = zeros(Nsymbrow, Nc);
+   
+    for c=1:Nc
+      corr = pilot2(:,c)' * ct_symb_buf(sampling_points,c);
+      mag  = sum(abs(ct_symb_buf(sampling_points,c)));
 
-    for p=1:Npilotsframe
-      % printf("%d %d %d %d %d\n", (p-1)*Ns+1, p*Ns, (p-1)*(Ns+1)+2, p*(Ns+1), (p-1)*(Ns+1)+1);
-      rx_symb((p-1)*Ns+1:p*Ns,:) = rx_symb_pilot((p-1)*(Ns+1)+2:p*(Ns+1),:);
-      rx_pilot(p,:) = rx_symb_pilot((p-1)*(Ns+1)+1,:);
+      phi_(:, c) = angle(corr);
+      amp_(:, c) = mag/length(sampling_points);
     end
 
-    % buffer three frames of symbols (and pilots) for phase recovery
-
-    rx_symb_buf(1:2*Nsymbrow,:) = rx_symb_buf(Nsymbrow+1:3*Nsymbrow,:);
-    rx_symb_buf(2*Nsymbrow+1:3*Nsymbrow,:) = rx_symb;
-    rx_pilot_buf(1:2*Npilotsframe,:) = rx_pilot_buf(Npilotsframe+1:3*Npilotsframe,:);
-    rx_pilot_buf(2*Npilotsframe+1:3*Npilotsframe,:) = rx_pilot;
-    sim_in.rx_symb_buf = rx_symb_buf;
-    sim_in.rx_pilot_buf = rx_pilot_buf;
-
-    % pilot assisted phase estimation and correction of middle frame in rx symb buffer
+    % now correct phase of data symbols and make decn on bits
 
-    rx_symb = rx_symb_buf(Nsymbrow+1:2*Nsymbrow,:);
-            
-    phi_ = zeros(Nsymbrow, Nc*Nchip);
-    amp_ = ones(Nsymbrow, Nc*Nchip);
-
-    for c=1:Nc*Nchip
-
-      if verbose > 2
-        printf("phi_   : ");
-      end
-
-      for r=1:Nsymbrow
-        st = Npilotsframe+1+floor((r-1)/Ns) - floor(Np/2) + 1;
-        en = st + Np - 1;
-        ch_est = tx_pilot_buf(st:en,c)'*rx_pilot_buf(st:en,c)/Np;
-        phi_(r,c) = angle(ch_est);
-        amp_(r,c) = abs(ch_est);
-        %amp_(r,c) = abs(rx_symb(r,c));
-        if verbose > 2
-          printf("% 4.3f ", phi_(r,c))
-        end
-        rx_symb(r,c) *= exp(-j*phi_(r,c));
-      end
-
-      if verbose > 2
-        printf("\nrx_symb: ");
-        for r=1:Nsymbrow
-          printf("% 4.3f ", angle(rx_symb(r,c)))
-        end
-        printf("\nindexes: ");
-        for r=1:Nsymbrow
-          st = Npilotsframe+1+floor((r-1)/Ns) - floor(Np/2) + 1;
-          en = st + Np - 1;
-          printf("%2d,%2d  ", st,en)
-        end
-        printf("\npilots : ");
-        for p=1:3*Npilotsframe
-          printf("% 4.3f ", angle(rx_pilot_buf(p,c)));
-        end 
-        printf("\n\n");
-      end
-    end 
-    
-    % de-spread
-            
-    for r=1:Nsymbrow
-      for c=Nc+1:Nc:Nchip*Nc
-        rx_symb(r,1:Nc) = rx_symb(r,1:Nc) + rx_symb(r,c:c+Nc-1);
-        amp_(r,1:Nc)    = amp_(r,1:Nc) + amp_(r,c:c+Nc-1);
-      end
-    end
-           
-    % demodulate stage 2
-
-    rx_symb_linear = zeros(1,Nsymb);
-    amp_linear = zeros(1,Nsymb);
+    rx_symb = zeros(Nsymbrow, Nc);
+    rx_symb_linear = zeros(1, Nsymbrow*Nc);
     rx_bits = zeros(1, framesize);
     for c=1:Nc
       for r=1:Nsymbrow
         i = (c-1)*Nsymbrow + r;
-        rx_symb_linear(i) = rx_symb(r,c);
-        amp_linear(i) = amp_(r,c);
+        rx_symb(i) = ct_symb_buf(2+r,c)*exp(-j*phi_(c));
+        rx_symb_linear(i) = rx_symb(i);
         rx_bits((2*(i-1)+1):(2*i)) = qpsk_demod(rx_symb(r,c));
       end
     end
@@ -337,18 +260,11 @@ function [rx_symb rx_bits rx_symb_linear amp_linear amp_ phi_ EsNo_ prev_sym_rx
 
     % Estimate signal power
     
-    Es_ = mean(amp_linear .^ 2);
+    Es_ = mean(amp_ .^ 2);
  
     EsNo_ = Es_/No_;
     %printf("Es_: %f No_: %f  Es/No: %f  Es/No dB: %f\n", Es_, No_, Es_/No_, 10*log10(EsNo_));
-  
-    % LDPC decoder requires some amplitude normalisation
-    % (AGC), was found to break ow.  So we adjust the symbol
-    % amplitudes so that they are an averge of 1
-
-    rx_symb_linear /= mean(amp_linear);
-    amp_linear /= mean(amp_linear);
-    
+      
 endfunction
 
 
@@ -625,6 +541,7 @@ function [next_sync cohpsk] = frame_sync_fine_timing_est(cohpsk, ch_symb, sync,
           max_mag = mag;
           cohpsk.ct = t;
           cohpsk.f_fine_est = f_fine;
+          cohpsk.ff_rect = exp(-j*f_fine*2*pi/Rs);
         end
       end
     end
index 5f1e26c6d90baa8e92a41686f2deb9632cc65ad2..593dbfb120a79a8db8f76fe873927343bb366ae3 100644 (file)
@@ -25,7 +25,7 @@ randn('state',1);
 n = 2000;
 frames = 35*4;
 framesize = 32;
-foff = -80;
+foff = 0;
 
 EsNodB = 8;
 EsNo = 10^(EsNodB/10);
@@ -233,36 +233,71 @@ for i=1:frames
   % coarse timing (frame sync) and initial fine freq est ---------------------------------------------
   
   [next_sync acohpsk] = frame_sync_fine_timing_est(acohpsk, ch_symb, sync, next_sync);
+  %acohpsk.ff_rect = exp(j*2*pi*(-1.0)/Rs);
 
   if sync == 1
     next_sync = 2;
   end
 
-  printf("i: %d sync: %d next_sync: %d\n", i, sync, next_sync);
+  %printf("i: %d sync: %d next_sync: %d\n", i, sync, next_sync);
   sync = next_sync;
 
-  if (i==10)
+  if (i==1000)
     xx
   end
 
-  if 0
-  [rx_symb rx_bits rx_symb_linear amp_linear amp_ phi_ EsNo_ prev_sym_rx acohpsk] = qpsk_symbols_to_bits(acohpsk, ct_symb_buf(ct+1:ct+acohpsk.Nsymbrowpilot,:), []);
-  rx_symb_log = [rx_symb_log; rx_symb];
-  rx_amp_log = [rx_amp_log; amp_];
-  rx_phi_log = [rx_phi_log; phi_];
-  rx_bits_log = [rx_bits_log rx_bits];
+  % We can decode first that we get sync on.  Need to fine freq correct all of it's symbols, 
+  % including pilots.  From then on, just correct new symbols into frame.
+  % make copy, so if we lose sync we havent fine freq corrected ct_symb_buf
+  % if next_sync == 4 correct all 8
+  % if sync == 2 correct latest 6
 
-  % BER stats
+  if (next_sync == 4) || (sync == 4)
 
-  if i > 3
-    error_positions = xor(prev_tx_bits2, rx_bits);
-    Nerrs  += sum(error_positions);
-    nerr_log = [nerr_log sum(error_positions)];
-    Tbits += length(error_positions);
+    if next_sync == 4
+
+      % first frame, we've just gotten sync so fine freq correct all Nsymbrowpilot+2 samples
+
+      ct_symb_ff_buf = acohpsk.ct_symb_buf(acohpsk.ct+1:acohpsk.ct+acohpsk.Nsymbrowpilot+2,:);
+      for r=1:acohpsk.Nsymbrowpilot+2
+        acohpsk.ff_phase *= acohpsk.ff_rect';
+        ct_symb_ff_buf(r,:) *= acohpsk.ff_phase;
+      end
+
+    else
+      
+      % second and subsequent frames, just fine freq correct the latest Nsymbrowpilot
+
+      ct_symb_ff_buf(1:2,:) = ct_symb_ff_buf(Nsymbrowpilot+1:Nsymbrowpilot+2,:);
+      ct_symb_ff_buf(3:Nsymbrowpilot+2,:) = acohpsk.ct_symb_buf(acohpsk.ct+3:acohpsk.ct+acohpsk.Nsymbrowpilot+2,:);
+      for r=3:acohpsk.Nsymbrowpilot+2
+        afdmdv.ff_phase *= afdmdv.ff_rect';
+        ct_symb_ff_buf(r,:) *= afdmdv.ff_phase;
+      end
+
+    end
+
+    mag = abs(acohpsk.ff_phase);
+    acohpsk.ff_phase /= mag;
+      
+    [rx_symb rx_bits amp_ phi_ EsNo_] = qpsk_symbols_to_bits(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_];
+    rx_bits_log = [rx_bits_log rx_bits];
+
+    % BER stats
+
+    if i > 2
+      error_positions = xor(prev_tx_bits2, rx_bits);
+      Nerrs  += sum(error_positions);
+      nerr_log = [nerr_log sum(error_positions)];
+     Tbits += length(error_positions);
+    end 
   end
+
   prev_tx_bits2 = prev_tx_bits;
   prev_tx_bits = tx_bits;
-  end
 
 end
 
@@ -281,6 +316,7 @@ stem_sig_and_error(6, 211, real(rx_symb_log_c(1:n)), real(rx_symb_log(1:n) - rx_
 stem_sig_and_error(6, 212, imag(rx_symb_log_c(1:n)), imag(rx_symb_log(1:n) - rx_symb_log_c(1:n)), 'rx symb im', [1 n -1.5 1.5])
 stem_sig_and_error(7, 111, rx_bits_log_c(1:n), rx_bits_log(1:n) - rx_bits_log_c(1:n), 'rx bits', [1 n -1.5 1.5])
 
+if 0
 check(tx_bits_log, tx_bits_log_c, 'tx_bits');
 check(tx_symb_log, tx_symb_log_c, 'tx_symb');
 check(tx_fdm_log, tx_fdm_log_c, 'tx_fdm');
@@ -292,6 +328,7 @@ check(rx_amp_log, rx_amp_log_c, 'rx_amp_log',0.01);
 check(rx_phi_log, rx_phi_log_c, 'rx_phi_log');
 check(rx_symb_log, rx_symb_log_c, 'rx_symb',0.01);
 check(rx_bits_log, rx_bits_log_c, 'rx_bits');
+end
 
 % Determine bit error rate