octave version of modem using lin reg for phase est, curves look OK

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

diff --git a/codec2-dev/octave/cohpsk.m b/codec2-dev/octave/cohpsk.m
index 96fbec69c6465fdd419e23bf3f24e628bb2c2370..0a62469461742a1f2d16e0f685c07921214f86c4 100644 (file)
--- a/codec2-dev/octave/cohpsk.m
+++ b/codec2-dev/octave/cohpsk.m
@@ -65,9 +65,11 @@ function sim_in = symbol_rate_init(sim_in)
     sim_in.Npilotsframe  = Npilotsframe     = 2;
     sim_in.Nsymbrowpilot = Nsymbrowpilot    = Nsymbrow + Npilotsframe;
 
-    printf("Each frame contains %d data bits or %d data symbols, transmitted as %d symbols by %d carriers.", framesize, Nsymb, Nsymbrow, Nc);
-    printf("  There are %d pilot symbols in each carrier together at the start of each frame, then %d data symbols.", Npilotsframe, Ns); 
-    printf("  Including pilots, the frame is %d symbols long by %d carriers.\n\n", Nsymbrowpilot, Nc);
+    if verbose == 2
+      printf("Each frame contains %d data bits or %d data symbols, transmitted as %d symbols by %d carriers.", framesize, Nsymb, Nsymbrow, Nc);
+      printf("  There are %d pilot symbols in each carrier together at the start of each frame, then %d data symbols.", Npilotsframe, Ns); 
+      printf("  Including pilots, the frame is %d symbols long by %d carriers.\n\n", Nsymbrowpilot, Nc);
+    end
 
     sim_in.prev_sym_tx = qpsk_mod([0 0])*ones(1,Nc*Nchip);
     sim_in.prev_sym_rx = qpsk_mod([0 0])*ones(1,Nc*Nchip);
@@ -226,6 +228,7 @@ function [rx_symb rx_bits rx_symb_linear amp_ phi_ EsNo_ cohpsk] = qpsk_symbols_
     Npilotsframe  = cohpsk.Npilotsframe;
     pilot         = cohpsk.pilot;
     verbose       = cohpsk.verbose;
+    coh_en        = cohpsk.coh_en;
 
     % 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
@@ -234,12 +237,20 @@ function [rx_symb rx_bits rx_symb_linear amp_ phi_ EsNo_ cohpsk] = qpsk_symbols_
     pilot2 = [cohpsk.pilot(1,:); cohpsk.pilot(2,:); cohpsk.pilot(1,:); cohpsk.pilot(2,:);];
     phi_ = zeros(Nsymbrow, Nc);
     amp_ = zeros(Nsymbrow, Nc);
-   
+    
+    % corr method was used initially, but was poor at tracking fast phase 
+    % changes.  Linear regression works better.
+
     for c=1:Nc
-      corr = pilot2(:,c)' * ct_symb_buf(sampling_points,c);
-      mag  = sum(abs(ct_symb_buf(sampling_points,c)));
+      %corr = pilot2(:,c)' * ct_symb_buf(sampling_points,c);      
+      %phi_(:, c) = angle(corr);
+     
+      y = ct_symb_buf(sampling_points,c) .* pilot2(:,c);
+      [m b] = linreg(sampling_points, y, length(sampling_points));
+      yfit = m*[3 4 5 6] + b;
+      phi_(:, c) = angle(yfit);
       
-      phi_(:, c) = angle(corr);
+      mag  = sum(abs(ct_symb_buf(sampling_points,c)));
       amp_(:, c) = mag/length(sampling_points);
     end
 
@@ -251,8 +262,11 @@ function [rx_symb rx_bits rx_symb_linear amp_ phi_ EsNo_ cohpsk] = qpsk_symbols_
     for c=1:Nc
       for r=1:Nsymbrow
         i = (c-1)*Nsymbrow + r;
-        rx_symb(r,c) = ct_symb_buf(2+r,c)*exp(-j*phi_(r,c));
-        %rx_symb(r,c) = ct_symb_buf(2+r,c);
+        if coh_en == 1
+          rx_symb(r,c) = ct_symb_buf(2+r,c)*exp(-j*phi_(r,c));
+        else
+          rx_symb(r,c) = ct_symb_buf(2+r,c);
+        end
         rx_symb_linear(i) = rx_symb(i);
         rx_bits((2*(i-1)+1):(2*i)) = qpsk_demod(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))));
@@ -682,6 +696,7 @@ function sim_out = ber_test(sim_in)
     hf_sim           = sim_in.hf_sim;
     nhfdelay         = sim_in.hf_delay_ms*Rs/1000;
     hf_mag_only      = sim_in.hf_mag_only;
+    f_off            = sim_in.f_off;
 
     [spread spread_2ms hf_gain] = init_hf_model(Fs, Rs, Nsymbrowpilot*Ntrials);
 
@@ -719,8 +734,9 @@ function sim_out = ber_test(sim_in)
 
         hf_n = 1;
 
-        phase_offset = 0;
-        w_offset     = pi/16;
+        phase_offset_rect = 1;
+        w_offset      = 2*pi*f_off/Rs;
+        w_offset_rect = exp(j*w_offset);
 
         ct_symb_buf = zeros(2*Nsymbrowpilot, Nc);
         prev_tx_symb = prev_rx_symb = ones(1,Nc);
@@ -750,7 +766,6 @@ function sim_out = ber_test(sim_in)
               tx_bits_buf(1:framesize) = tx_bits;
             end
 
-
             s_ch = tx_symb;
 
             % HF channel simulation  ------------------------------------
@@ -797,8 +812,12 @@ function sim_out = ber_test(sim_in)
             noise = sqrt(variance*0.5)*(randn(Nsymbrowpilot,Nc*Nchip) + j*randn(Nsymbrowpilot,Nc*Nchip));
             noise_log = [noise_log noise];
 
-            s_ch = s_ch + noise;
-            
+            for r=1:Nsymbrowpilot
+              s_ch(r,:) *= phase_offset_rect;
+              phase_offset_rect *= w_offset_rect;
+            end
+            s_ch += noise;
+
             ct_symb_buf(1:Nsymbrowpilot,:) = ct_symb_buf(Nsymbrowpilot+1:2*Nsymbrowpilot,:);
             ct_symb_buf(Nsymbrowpilot+1:2*Nsymbrowpilot,:) = s_ch;
 
@@ -887,7 +906,7 @@ function sim_out = ber_test(sim_in)
           EsNodBSurface(find(EsNodBSurface < -5)) = -5;
           mesh(x,y,EsNodBSurface);
           grid
-          axis([1 (Nc+1)*Nchip 1 Rs*5 -5 15])
+          axis([1 (Nc+1)*Nchip 1 Rs*5 -5 25])
           title('HF Channel Es/No');
 
           if verbose 
@@ -907,7 +926,7 @@ function sim_out = ber_test(sim_in)
          p = Ns;
          for r=1:m1
            if p == Ns
-             phi_x_counter++;
+             phi_x_counter += Npilotsframe;
              p = 0;
            end
            p++;
@@ -919,6 +938,7 @@ function sim_out = ber_test(sim_in)
          figure(5);
          clf
          subplot(211)
+         [m n] = size(phi_log);
          plot(phi_x, phi_log(:,2),'r+;Estimated HF channel phase;')
          if hf_sim
            hold on;
@@ -928,6 +948,7 @@ function sim_out = ber_test(sim_in)
          end
          ylabel('Phase (rads)');
          legend('boxoff');
+         axis([1 m -1.1*pi 1.1*pi])
 
          subplot(212)
          plot(phi_x, amp_log(:,2),'r+;Estimated HF channel amp;')
@@ -939,17 +960,13 @@ function sim_out = ber_test(sim_in)
          ylabel('Amplitude');
          xlabel('Time (symbols)');
          legend('boxoff');
+         axis([1 m 0 3])
        end
 
        figure(4)
        clf
-       subplot(211)
        stem(Nerrs_log)
-       subplot(212)
-       if ldpc_code
-         stem(ldpc_Nerrs_log)
-       end
-
+       axis([1 length(Nerrs_log) 0 max(Nerrs_log)+1])
    end
 
 endfunction
@@ -958,6 +975,7 @@ endfunction
 
 function sim_in = standard_init
   sim_in.verbose          = 1;
+  sim_in.do_write_pilot_file = 0;
   sim_in.plot_scatter     = 0;
 
   sim_in.Esvec            = 50; 

diff --git a/codec2-dev/octave/linreg.m b/codec2-dev/octave/linreg.m
new file mode 100644 (file)
index 0000000..7a44b8b
--- /dev/null
+++ b/codec2-dev/octave/linreg.m
@@ -0,0 +1,33 @@
+% linreg.m
+% David Rowe April 2015
+%
+% Based on:
+%    http://stackoverflow.com/questions/5083465/fast-efficient-least-squares-fit-algorithm-in-c
+
+function [m b] = linreg(x,y,n)
+  sumx = 0.0;   % sum of x
+  sumx2 = 0.0;  % sum of x^2
+  sumxy = 0.0;  % sum of x * y
+  sumy = 0.0;   % sum of y
+  sumy2 = 0.0;  % sum of y**2
+
+  for i=1:n   
+    sumx  += x(i);       
+    sumx2 += x(i)^2;  
+    sumxy += x(i) * y(i);
+    sumy  += y(i);      
+    sumy2 += y(i)^2; 
+  end 
+
+  denom = (n * sumx2 - sumx*sumx);
+
+  if denom == 0
+    % singular matrix. can't solve the problem.
+    m = 0;
+    b = 0;
+  else
+    m = (n * sumxy  -  sumx * sumy) / denom;
+    b = (sumy * sumx2  -  sumx * sumxy) / denom;
+  end
+
+endfunction

diff --git a/codec2-dev/octave/tcohpsk.m b/codec2-dev/octave/tcohpsk.m
index 109f125d1d5fd7c051c2d142443e076a40934834..d3cb239094120ce611e8debe8aacbfe2ee6dbe5d 100644 (file)
--- a/codec2-dev/octave/tcohpsk.m
+++ b/codec2-dev/octave/tcohpsk.m
@@ -78,6 +78,7 @@ acohpsk.ldpc_code_rate   = 1;
 acohpsk.Nc               = Nc;
 acohpsk.Rs               = Rs;
 acohpsk.Ns               = 4;
+acohpsk.coh_en           = 1;
 acohpsk.Nchip            = 1;
 acohpsk.modulation       = 'qpsk';
 acohpsk.do_write_pilot_file = 0;

diff --git a/codec2-dev/octave/test_cohpsk.m b/codec2-dev/octave/test_cohpsk.m
index 50832bceefa6f9ebdf1cfb454511265e325830ae..cb83e688cd7a202536063c130c93f34e4d5b11eb 100644 (file)
--- a/codec2-dev/octave/test_cohpsk.m
+++ b/codec2-dev/octave/test_cohpsk.m
@@ -36,9 +36,8 @@ function test_curves
   sim_in.verbose          = 1;
   sim_in.plot_scatter     = 1;
 
-  sim_in.Esvec            = 10; 
+  sim_in.Esvec            = 20; 
   sim_in.framesize        = 32;
-  sim_in.hf_sim           = 1;
   sim_in.Ntrials          = 100;
   sim_in.Rs               = 50;
   sim_in.Nc               = 4;
@@ -48,12 +47,17 @@ function test_curves
   sim_in.modulation       = 'qpsk';
   sim_in.ldpc_code_rate   = 1;
   sim_in.ldpc_code        = 0;
+  sim_in.coh_en           = 1;
+
+  sim_in.hf_sim           = 1;
+  sim_in.hf_mag_only      = 0;
+  sim_in.f_off            = 0;
 
   sim_qpsk                = ber_test(sim_in);
 
   % AWGN curves ----------------------------------------------------
 
-  sim_in.Ntrials          = 500;
+  sim_in.Ntrials          = 400;
   sim_in.hf_sim           = 0;
   sim_in.plot_scatter     = 0;
   sim_in.Esvec            = 5:10; 
@@ -65,26 +69,24 @@ function test_curves
   sim_dqpsk               = ber_test(sim_in, 'dqpsk');
 
   sim_in.modulation       = 'qpsk';
-  sim_in.Ns               = 4;
-  sim_in.Np               = 2;
   sim_qpsk_pilot          = ber_test(sim_in, 'qpsk');
 
   % HF curves ----------------------------------------------------
 
-  sim_in.Ntrials          = 200;
+  sim_in.Ntrials          = 400;
   sim_in.hf_sim           = 1;
   sim_in.plot_scatter     = 0;
   sim_in.Esvec            = 5:20; 
-
-  Ebvec = sim_in.Esvec - 10*log10(2);
-  BER_theory = 0.5*erfc(sqrt(10.^(Ebvec/10)));
-  
   sim_in.modulation       = 'dqpsk';
   sim_dqpsk_hf            = ber_test(sim_in, 'dqpsk');
 
   sim_in.modulation       = 'qpsk';
-  sim_in.Ns               = 4;
-  sim_in.Np               = 2;
+  sim_in.coh_en           = 0;
+  sim_in.hf_mag_only      = 1;
+  sim_qpsk_hf             = ber_test(sim_in, 'qpsk');
+
+  sim_in.coh_en           = 1;
+  sim_in.hf_mag_only      = 0;
   sim_qpsk_pilot_hf       = ber_test(sim_in, 'qpsk');
 
   % plot results ---------------------------------------------------
@@ -97,16 +99,16 @@ function test_curves
   semilogy(sim_dqpsk.Ebvec, sim_dqpsk.BERvec,'c;DQPSK AWGN;')
   semilogy(sim_qpsk_pilot.Ebvec, sim_qpsk_pilot.BERvec,'b;QPSK pilot AWGN;')
 
-  %semilogy(sim_qpsk_hf_ideal.Ebvec, sim_qpsk_hf_ideal.BERvec,'b;QPSK HF ideal;')
-  semilogy(sim_dqpsk_hf.Ebvec, sim_dqpsk_hf.BERvec,'c;DQPSK HF;')
+  semilogy(sim_qpsk_hf.Ebvec, sim_qpsk_hf.BERvec,'r;QPSK HF ideal;')
   semilogy(sim_qpsk_pilot_hf.Ebvec, sim_qpsk_pilot_hf.BERvec,'b;QPSK pilot HF;')
+  semilogy(sim_dqpsk_hf.Ebvec, sim_dqpsk_hf.BERvec,'c;DQPSK HF;')
 
   hold off;
 
   xlabel('Eb/N0')
   ylabel('BER')
   grid("minor")
-  axis([min(Ebvec) max(Ebvec) 1E-3 1])
+  axis([min(sim_qpsk_hf.Ebvec) max(sim_qpsk_hf.Ebvec) 1E-3 1])
   legend("boxoff");
 endfunction
 
@@ -129,13 +131,13 @@ function test_single
 
   sim_in.Ntrials          = 100;
   sim_in.Esvec            = 10; 
-  sim_in.hf_sim           = 0;
+  sim_in.hf_sim           = 1;
   sim_in.hf_mag_only      = 0;
   sim_in.modulation       = 'qpsk';
+  sim_in.coh_en           = 1;
+  sim_in.f_off            = 0;
 
-  sim_in.modulation       = 'dqpsk';
-
-  sim_in.do_write_pilot_file = 0;
+  %sim_in.modulation      = 'dqpsk';
 
   sim_qpsk_hf             = ber_test(sim_in);
 
@@ -593,6 +595,7 @@ endfunction
 % Start simulations ---------------------------------------
 
 more off;
+%close all;
 test_curves();
 %test_single();
 %rate_Fs_tx("tx_zero.raw");

diff --git a/codec2-dev/octave/tlinreg.m b/codec2-dev/octave/tlinreg.m
index 81656e687151b0108b216385ca0a09e76241ada2..1cde350d0dc4dee456348802144799c7c7e830a7 100644 (file)
--- a/codec2-dev/octave/tlinreg.m
+++ b/codec2-dev/octave/tlinreg.m
@@ -1,33 +1,7 @@
-
-1;
-
-function [m b] = linreg(x,y,n)
-  sumx = 0.0;   % sum of x
-  sumx2 = 0.0;  % sum of x^2
-  sumxy = 0.0;  % sum of x * y
-  sumy = 0.0;   % sum of y
-  sumy2 = 0.0;  % sum of y**2
-
-  for i=1:n   
-    sumx  += x(i);       
-    sumx2 += x(i)^2;  
-    sumxy += x(i) * y(i);
-    sumy  += y(i);      
-    sumy2 += y(i)^2; 
-  end 
-
-  denom = (n * sumx2 - sumx*sumx);
-
-  if denom == 0
-    % singular matrix. can't solve the problem.
-    m = 0;
-    b = 0;
-  else
-    m = (n * sumxy  -  sumx * sumy) / denom;
-    b = (sumy * sumx2  -  sumx * sumxy) / denom;
-  end
-
-endfunction
+% tlinreg
+% David Rowe April 2015
+%
+% Unit test for linear regression
 
 x = [1 2 7 8];
 y = [ -0.70702 + 0.70708i 0.77314 - 0.63442i -0.98083 + 0.19511i 0.99508 - 0.09799i] .* [1 -1 1 -1];