two path diveristy looking gd in Octave, starting to port into C
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Fri, 17 Apr 2015 00:59:08 +0000 (00:59 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Fri, 17 Apr 2015 00:59:08 +0000 (00:59 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@2119 01035d8c-6547-0410-b346-abe4f91aad63

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

index 0a62469461742a1f2d16e0f685c07921214f86c4..7f76e876100db284d0d2b942d6d5fb2491e1205a 100644 (file)
@@ -53,7 +53,7 @@ function sim_in = symbol_rate_init(sim_in)
     nhfdelay         = sim_in.hf_delay_ms*Rs/1000;
     hf_mag_only      = sim_in.hf_mag_only;
 
-    Nchip            = sim_in.Nchip;  % spread spectrum factor
+    Nd               = sim_in.Nd;     % diveristy
     Ns               = sim_in.Ns;     % step size between pilots
     ldpc_code        = sim_in.ldpc_code;
     rate             = sim_in.ldpc_code_rate; 
@@ -71,11 +71,11 @@ function sim_in = symbol_rate_init(sim_in)
       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);
+    sim_in.prev_sym_tx = qpsk_mod([0 0])*ones(1,Nc*Nd);
+    sim_in.prev_sym_rx = qpsk_mod([0 0])*ones(1,Nc*Nd);
 
-    sim_in.rx_symb_buf  = zeros(3*Nsymbrow, Nc*Nchip);
-    sim_in.rx_pilot_buf = zeros(3*Npilotsframe,Nc*Nchip);
+    sim_in.rx_symb_buf  = zeros(3*Nsymbrow, Nc*Nd);
+    sim_in.rx_pilot_buf = zeros(3*Npilotsframe,Nc*Nd);
     sim_in.tx_bits_buf  = zeros(1,2*framesize);
 
     % pilot sequence is used for phase and amplitude estimation, and frame sync
@@ -90,11 +90,11 @@ function sim_in = symbol_rate_init(sim_in)
     % we use first 2 pilots of next frame to help with frame sync and fine freq
 
     sim_in.Nct_sym_buf = 2*Nsymbrowpilot + 2;
-    sim_in.ct_symb_buf = zeros(sim_in.Nct_sym_buf, Nc);
+    sim_in.ct_symb_buf = zeros(sim_in.Nct_sym_buf, Nc*Nd);
 
     sim_in.ff_phase = 1;
 
-    sim_in.ct_symb_ff_buf = zeros(Nsymbrowpilot + 2, Nc);
+    sim_in.ct_symb_ff_buf = zeros(Nsymbrowpilot + 2, Nc*Nd);
 
     % Init LDPC --------------------------------------------------------------------
 
@@ -179,9 +179,9 @@ function [tx_symb tx_bits] = bits_to_qpsk_symbols(sim_in, tx_bits, code_param)
     Nc            = sim_in.Nc;
     Npilotsframe  = sim_in.Npilotsframe;
     Ns            = sim_in.Ns;
-    Nchip         = sim_in.Nchip;
     modulation    = sim_in.modulation;
     pilot         = sim_in.pilot;
+    Nd            = sim_in.Nd;
 
     if ldpc_code
         [tx_bits, tmp] = ldpc_enc(tx_bits, code_param);
@@ -205,15 +205,17 @@ function [tx_symb tx_bits] = bits_to_qpsk_symbols(sim_in, tx_bits, code_param)
     
     tx_symb = [pilot(1,:); pilot(2,:); tx_symb;];
 
-    % Optionally copy to other carriers (spreading)
+    % Optionally copy to other carriers (diversity)
 
-    for c=Nc+1:Nc:Nc*Nchip
-      tx_symb(:,c:c+Nc-1) = tx_symb(:,1:Nc);
+    tmp = tx_symb;
+    for d=1:Nd-1
+      tmp = [tmp tx_symb];
     end
+    tx_symb = tmp;
 
-    % ensures energy/symbol is normalised when spreading
+    % ensures energy/symbol is normalised with diveristy
 
-    tx_symb = tx_symb/sqrt(Nchip);
+    tx_symb = tx_symb/sqrt(Nd);
 end
 
 
@@ -225,27 +227,30 @@ function [rx_symb rx_bits rx_symb_linear amp_ phi_ EsNo_ cohpsk] = qpsk_symbols_
     Nsymbrow      = cohpsk.Nsymbrow;
     Nsymbrowpilot = cohpsk.Nsymbrowpilot;
     Nc            = cohpsk.Nc;
+    Nd            = cohpsk.Nd;
     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
+    % Use pilots to get phase and amplitude estimates We assume there
+    % are two samples at the start of each frame and two at the end
+    % Note: correlation (averging) method was used initially, but was
+    % poor at tracking fast phase changes that we experience on fading
+    % channels.  Linear regression (fitting a straight line) works
+    % better on fading channels, but increases BER slighlty for AWGN
+    % channels.
 
     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);
     
-    % corr method was used initially, but was poor at tracking fast phase 
-    % changes.  Linear regression works better.
-
-    for c=1:Nc
+    for c=1:Nc*Nd
       %corr = pilot2(:,c)' * ct_symb_buf(sampling_points,c);      
       %phi_(:, c) = angle(corr);
      
-      y = ct_symb_buf(sampling_points,c) .* pilot2(:,c);
+      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*[3 4 5 6] + b;
       phi_(:, c) = angle(yfit);
@@ -254,25 +259,37 @@ function [rx_symb rx_bits rx_symb_linear amp_ phi_ EsNo_ cohpsk] = qpsk_symbols_
       amp_(:, c) = mag/length(sampling_points);
     end
 
-    % now correct phase of data symbols and make decn on bits
+    % now correct phase of data symbols
 
     rx_symb = zeros(Nsymbrow, Nc);
     rx_symb_linear = zeros(1, Nsymbrow*Nc);
     rx_bits = zeros(1, framesize);
-    for c=1:Nc
+    for c=1:Nc*Nd
       for r=1:Nsymbrow
-        i = (c-1)*Nsymbrow + r;
         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));
+        i = (c-1)*Nsymbrow + r;
         %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 bits
+
+    for c=1:Nc
+      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
+        rx_symb_linear(i) = div_symb;
+        rx_bits((2*(i-1)+1):(2*i)) = qpsk_demod(div_symb);
+      end
+    end
+
     % 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
@@ -684,7 +701,7 @@ function sim_out = ber_test(sim_in)
     Npilotsframe     = sim_in.Npilotsframe;
     Ns               = sim_in.Ns;
     Np               = sim_in.Np;
-    Nchip            = sim_in.Nchip;
+    Nd               = sim_in.Nd;
     modulation       = sim_in.modulation;
     pilot            = sim_in.pilot;
     prev_sym_tx      = sim_in.prev_sym_tx;
@@ -697,8 +714,9 @@ function sim_out = ber_test(sim_in)
     nhfdelay         = sim_in.hf_delay_ms*Rs/1000;
     hf_mag_only      = sim_in.hf_mag_only;
     f_off            = sim_in.f_off;
+    div_time_shift   = sim_in.div_timeshift;
 
-    [spread spread_2ms hf_gain] = init_hf_model(Fs, Rs, Nsymbrowpilot*Ntrials);
+    [spread spread_2ms hf_gain] = init_hf_model(Fs, Rs, Nsymbrowpilot*(Ntrials+2));
 
     if strcmp(modulation,'dqpsk')
       Nsymbrowpilot = Nsymbrow;
@@ -738,8 +756,8 @@ function sim_out = ber_test(sim_in)
         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);
+        ct_symb_buf = zeros(2*Nsymbrowpilot, Nc*Nd);
+        prev_tx_symb = prev_rx_symb = ones(1, Nc*Nd);
 
         % simulation starts here-----------------------------------
  
@@ -780,11 +798,15 @@ function sim_out = ber_test(sim_in)
 
                 wsep = 2*pi*(1+0.5);  % e.g. 75Hz spacing at Rs=50Hz, alpha=0.5 filters
 
-                hf_model(hf_n, :) = zeros(1,Nc*Nchip);
+                hf_model(hf_n, :) = zeros(1,Nc*Nd);
                 
                 for r=1:Nsymbrowpilot
-                  for c=1:Nchip*Nc
-                    time_shift = floor((c-1)*Nsymbrowpilot);
+                  for c=1:Nd*Nc
+                    if c > Nc 
+                      time_shift = sim_in.div_timeshift;
+                    else
+                      time_shift = 1;
+                    end
                     ahf_model = hf_gain*(spread(hf_n+time_shift) + exp(-j*c*wsep*nhfdelay)*spread_2ms(hf_n+time_shift));
                     
                     if hf_mag_only
@@ -801,7 +823,7 @@ function sim_out = ber_test(sim_in)
             % keep a record of each tx symbol so we can check average power
 
             for r=1:Nsymbrow
-              for c=1:Nchip*Nc
+              for c=1:Nd*Nc
                  s_ch_tx_log = [s_ch_tx_log s_ch(r,c)];
               end
             end
@@ -809,7 +831,7 @@ function sim_out = ber_test(sim_in)
             % AWGN noise and phase/freq offset channel simulation
             % 0.5 factor ensures var(noise) == variance , i.e. splits power between Re & Im
 
-            noise = sqrt(variance*0.5)*(randn(Nsymbrowpilot,Nc*Nchip) + j*randn(Nsymbrowpilot,Nc*Nchip));
+            noise = sqrt(variance*0.5)*(randn(Nsymbrowpilot,Nc*Nd) + j*randn(Nsymbrowpilot,Nc*Nd));
             noise_log = [noise_log noise];
 
             for r=1:Nsymbrowpilot
@@ -901,12 +923,13 @@ function sim_out = ber_test(sim_in)
           clf;
         
           y = 1:(hf_n-1);
-          x = 1:Nc*Nchip;
+          x = 1:Nc*Nd;
           EsNodBSurface = 20*log10(abs(hf_model(y,:))) - 10*log10(variance);
           EsNodBSurface(find(EsNodBSurface < -5)) = -5;
+          EsNodBSurface(find(EsNodBSurface > 25)) = 25;
           mesh(x,y,EsNodBSurface);
           grid
-          axis([1 (Nc+1)*Nchip 1 Rs*5 -5 25])
+          axis([1 Nc*Nd 1 Rs*5 -5 25])
           title('HF Channel Es/No');
 
           if verbose 
@@ -991,5 +1014,5 @@ function sim_in = standard_init
   sim_in.hf_sim           = 0;
   sim_in.hf_mag_only      = 0;
 
-  sim_in.Nchip            = 1;
+  sim_in.Nd            = 1;
 endfunction
index d3cb239094120ce611e8debe8aacbfe2ee6dbe5d..3b381abe98ac1086351697b5f49ec90218185c8b 100644 (file)
@@ -30,12 +30,13 @@ EsNo = 10^(EsNodB/10);
 
 Rs = 50;
 Nc = 4;
+Nd = 2;
 framesize = 32;
 
 % FDMDV init ---------------------------------------------------------------
 
 afdmdv.Fs = 8000;
-afdmdv.Nc = Nc-1;
+afdmdv.Nc = Nd*Nc-1;
 afdmdv.Rs = Rs;
 afdmdv.M  = afdmdv.Fs/afdmdv.Rs;
 afdmdv.tx_filter_memory = zeros(afdmdv.Nc+1, Nfilter);
@@ -43,7 +44,7 @@ afdmdv.Nfilter =  Nfilter;
 afdmdv.gt_alpha5_root = gt_alpha5_root;
 afdmdv.Fsep = 75;
 afdmdv.phase_tx = ones(afdmdv.Nc+1,1);
-freq_hz = Fsep*( -Nc/2 - 0.5 + (1:Nc) );
+freq_hz = Fsep*( -afdmdv.Nc/2 - 0.5 + (1:afdmdv.Nc+1) );
 afdmdv.freq_pol = 2*pi*freq_hz/Fs;
 afdmdv.freq = exp(j*afdmdv.freq_pol);
 afdmdv.Fcentre = 1500;
@@ -79,7 +80,7 @@ acohpsk.Nc               = Nc;
 acohpsk.Rs               = Rs;
 acohpsk.Ns               = 4;
 acohpsk.coh_en           = 1;
-acohpsk.Nchip            = 1;
+acohpsk.Nd               = Nd;
 acohpsk.modulation       = 'qpsk';
 acohpsk.do_write_pilot_file = 0;
 acohpsk = symbol_rate_init(acohpsk);
@@ -185,7 +186,7 @@ for i=1:frames
 
   % sample rate demod processing
 
-  ch_symb = zeros(acohpsk.Nsymbrowpilot, Nc);
+  ch_symb = zeros(acohpsk.Nsymbrowpilot, Nc*Nd);
   for r=1:acohpsk.Nsymbrowpilot
 
     % downconvert each FDM carrier to Nc separate baseband signals
@@ -233,6 +234,7 @@ for i=1:frames
 
 end
 
+if 0
 stem_sig_and_error(1, 111, tx_bits_log_c(1:n), tx_bits_log(1:n) - tx_bits_log_c(1:n), 'tx bits', [1 n -1.5 1.5])
 stem_sig_and_error(2, 211, real(tx_symb_log_c(1:n)), real(tx_symb_log(1:n) - tx_symb_log_c(1:n)), 'tx symb re', [1 n -1.5 1.5])
 stem_sig_and_error(2, 212, imag(tx_symb_log_c(1:n)), imag(tx_symb_log(1:n) - tx_symb_log_c(1:n)), 'tx symb im', [1 n -1.5 1.5])
@@ -275,7 +277,8 @@ Nerrs_c = sum(xor(tx_bits_prev_log, rx_bits_log_c));
 Tbits_c = length(tx_bits_prev_log);
 ber_c = Nerrs_c/Tbits_c;
 ber = Nerrs/Tbits;
-printf("EsNodB: %4.1f ber..: %3.2f Nerrs..: %d Tbits..: %d\n", EsNodB, ber, Nerrs, Tbits);
+end
+
 printf("EsNodB: %4.1f ber_c: %3.2f Nerrs_c: %d Tbits_c: %d\n", EsNodB, ber_c, Nerrs_c, Tbits_c);
 
 % some other useful plots
index cb83e688cd7a202536063c130c93f34e4d5b11eb..c680f267ba70b4afd73206c00b349ed6fcf63469 100644 (file)
@@ -38,12 +38,12 @@ function test_curves
 
   sim_in.Esvec            = 20; 
   sim_in.framesize        = 32;
-  sim_in.Ntrials          = 100;
+  sim_in.Ntrials          = 400;
   sim_in.Rs               = 50;
   sim_in.Nc               = 4;
   sim_in.Np               = 2;
   sim_in.Ns               = 4;
-  sim_in.Nchip            = 1;
+  sim_in.Nd               = 1;
   sim_in.modulation       = 'qpsk';
   sim_in.ldpc_code_rate   = 1;
   sim_in.ldpc_code        = 0;
@@ -52,6 +52,7 @@ function test_curves
   sim_in.hf_sim           = 1;
   sim_in.hf_mag_only      = 0;
   sim_in.f_off            = 0;
+  sim_in.div_timeshift    = 1;
 
   sim_qpsk                = ber_test(sim_in);
 
@@ -89,6 +90,19 @@ function test_curves
   sim_in.hf_mag_only      = 0;
   sim_qpsk_pilot_hf       = ber_test(sim_in, 'qpsk');
 
+  sim_in.Nd               = 2;
+  sim_in.div_timeshift    = 1;
+  sim_qpsk_pilot_hf_div   = ber_test(sim_in, 'qpsk');
+
+  sim_in.div_timeshift    = sim_in.Rs;
+  sim_qpsk_pilot_hf_div2  = ber_test(sim_in, 'qpsk');
+
+  sim_in.modulation       = 'qpsk';
+  sim_in.coh_en           = 0;
+  sim_in.hf_mag_only      = 1;
+  sim_in.div_timeshift    = 1;
+  sim_qpsk_hf_div         = ber_test(sim_in, 'qpsk');
+
   % plot results ---------------------------------------------------
 
   figure(1); 
@@ -100,7 +114,10 @@ function test_curves
   semilogy(sim_qpsk_pilot.Ebvec, sim_qpsk_pilot.BERvec,'b;QPSK pilot AWGN;')
 
   semilogy(sim_qpsk_hf.Ebvec, sim_qpsk_hf.BERvec,'r;QPSK HF ideal;')
+  semilogy(sim_qpsk_hf_div.Ebvec, sim_qpsk_hf_div.BERvec,'r;QPSK HF ideal div;')
   semilogy(sim_qpsk_pilot_hf.Ebvec, sim_qpsk_pilot_hf.BERvec,'b;QPSK pilot HF;')
+  semilogy(sim_qpsk_pilot_hf_div.Ebvec, sim_qpsk_pilot_hf_div.BERvec,'g;QPSK pilot Nd=2 HF;')
+  semilogy(sim_qpsk_pilot_hf_div2.Ebvec, sim_qpsk_pilot_hf_div2.BERvec,'g;QPSK pilot Nd=2 Rs HF;')
   semilogy(sim_dqpsk_hf.Ebvec, sim_dqpsk_hf.BERvec,'c;DQPSK HF;')
 
   hold off;
@@ -125,18 +142,20 @@ function test_single
   sim_in.Rs               = 50;
   sim_in.Ns               = 4;
   sim_in.Np               = 2;
-  sim_in.Nchip            = 1;
+  sim_in.Nd               = 2;
   sim_in.ldpc_code_rate   = 1;
   sim_in.ldpc_code        = 0;
 
-  sim_in.Ntrials          = 100;
-  sim_in.Esvec            = 10
+  sim_in.Ntrials          = 400;
+  sim_in.Esvec            = 13
   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.div_timeshift    = 1;
+
   %sim_in.modulation      = 'dqpsk';
 
   sim_qpsk_hf             = ber_test(sim_in);
@@ -596,8 +615,8 @@ endfunction
 
 more off;
 %close all;
-test_curves();
-%test_single();
+%test_curves();
+test_single();
 %rate_Fs_tx("tx_zero.raw");
 %rate_Fs_tx("tx.raw");
 %rate_Fs_rx("tx_-4dB.wav")