added code to output error pattern files for initial tests on speech quality through...

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

diff --git a/codec2-dev/octave/test_qpsk3.m b/codec2-dev/octave/test_qpsk3.m
index 60d355bfd835914b898dccdd3d0949b00f10fd12..37efad87ec197c4691d838e0cfec9c6e1de0d4c3 100644 (file)
--- a/codec2-dev/octave/test_qpsk3.m
+++ b/codec2-dev/octave/test_qpsk3.m
@@ -1,12 +1,15 @@
-% test_qps2k.m\r
-% David Rowe Feb 2014\r
+% test_qps3k.m\r
+% David Rowe March 2014\r
 %\r
 % QPSK modem simulation, version 2.  Simplifed version of\r
-% test_qpsk. initially based on code by Bill Cowley Generates curves\r
+% test_qpsk. Initially based on code by Bill Cowley Generates curves\r
 % BER versus E/No curves for different modems.  Design to test\r
 % coherent demodulation ideas on HF channels without building a full\r
 % blown modem.  Uses 'genie provided' estimates for timing estimation,\r
 % frame sync.\r
+%\r
+% Compared to test_qsk2.m this version supports phase estimation\r
+% (coherent demod)\r
 \r
 1;\r
 \r
@@ -30,11 +33,14 @@ function sim_out = ber_test(sim_in, modulation)
     hf_phase_only    = sim_in.hf_phase_only;\r
     hf_mag_only      = sim_in.hf_mag_only;\r
     Nc               = sim_in.Nc;\r
+    sim_coh_dpsk     = sim_in.sim_coh_dpsk;\r
 \r
     bps              = 2;\r
     Nsymb            = framesize/bps;\r
-    prev_sym_tx      = qpsk_mod([0 0]);\r
-    prev_sym_rx      = qpsk_mod([0 0]);\r
+    for k=1:Nc\r
+      prev_sym_tx(k) = qpsk_mod([0 0]);\r
+      prev_sym_rx(k) = qpsk_mod([0 0]);\r
+    end\r
 \r
     phase_est_method = sim_in.phase_est_method;\r
     if phase_est_method == 2\r
@@ -78,14 +84,14 @@ function sim_out = ber_test(sim_in, modulation)
 \r
         rate = sim_in.ldpc_code_rate; \r
         mod_order = 4; \r
-        modulation = 'QPSK';\r
+        modulation2 = 'QPSK';\r
         mapping = 'gray';\r
 \r
         demod_type = 0;\r
         decoder_type = 0;\r
         max_iterations = 100;\r
 \r
-        code_param = ldpc_init(rate, framesize, modulation, mod_order, mapping);\r
+        code_param = ldpc_init(rate, framesize, modulation2, mod_order, mapping);\r
         code_param.code_bits_per_frame = framesize;\r
         code_param.symbols_per_frame = framesize/bps;\r
     else\r
@@ -150,7 +156,7 @@ function sim_out = ber_test(sim_in, modulation)
         rx_symb_log      = [];\r
         noise_log        = [];\r
         mod_strip_log    = [];\r
-\r
+        \r
         % init HF channel\r
 \r
         hf_n = 1;\r
@@ -158,6 +164,9 @@ function sim_out = ber_test(sim_in, modulation)
         hf_fading = ones(1,Nsymb);             % default input for ldpc dec\r
         hf_model = ones(Ntrials*Nsymb/Nc, Nc); % defaults for plotting surface\r
 \r
+        sim_out.errors_log = [];\r
+        sim_out.ldpc_errors_log = [];\r
+\r
         for nn = 1: Ntrials\r
                   \r
             tx_bits = round( rand( 1, framesize*rate ) );\r
@@ -166,16 +175,17 @@ function sim_out = ber_test(sim_in, modulation)
 \r
             if ldpc_code\r
                 [tx_bits, s] = ldpc_enc(tx_bits, code_param);\r
-            else\r
-                s = zeros(1, Nsymb);\r
-                for i=1:Nsymb\r
-                    tx_symb = qpsk_mod(tx_bits(2*(i-1)+1:2*i));\r
-                    if strcmp(modulation,'dqpsk')\r
-                        tx_symb *= prev_sym_tx;\r
-                        prev_sym_tx = tx_symb;\r
-                    end \r
-                    s(i) = tx_symb;\r
-                end\r
+            end\r
+            s = zeros(1, Nsymb);\r
+            for i=1:Nc:Nsymb\r
+              for k=1:Nc\r
+                tx_symb = qpsk_mod(tx_bits(2*(i-1+k-1)+1:2*(i+k-1)));\r
+                if strcmp(modulation,'dqpsk')\r
+                    tx_symb *= prev_sym_tx(k);\r
+                    prev_sym_tx(k) = tx_symb;\r
+                end \r
+                s(i+k-1) = tx_symb;\r
+              end\r
             end\r
             tx_bits_buf(1:framesize) = tx_bits_buf(framesize+1:2*framesize);\r
             tx_bits_buf(framesize+1:2*framesize) = tx_bits;\r
@@ -270,7 +280,7 @@ function sim_out = ber_test(sim_in, modulation)
                         % s_ch(i+k-1) = r_delay_line(k,centre);\r
                      end\r
 \r
-                      if phase_est_method == 3\r
+                     if phase_est_method == 3\r
                         % QPSK modulation strip and phase est with original symbol mags\r
 \r
                         mod_strip_pol  = angle(r_delay_line(k,:)) * 4;\r
@@ -308,15 +318,22 @@ function sim_out = ber_test(sim_in, modulation)
             % de-modulate\r
 \r
             rx_bits = zeros(1, framesize);\r
-            for i=1:Nsymb\r
-                rx_symb = s_ch(i);\r
+            for i=1:Nc:Nsymb\r
+              for k=1:Nc\r
+                rx_symb = s_ch(i+k-1);\r
                 if strcmp(modulation,'dqpsk')\r
                     tmp = rx_symb;\r
-                    rx_symb *= conj(prev_sym_rx/abs(prev_sym_rx));\r
-                    prev_sym_rx = tmp;\r
+                    rx_symb *= conj(prev_sym_rx(k)/abs(prev_sym_rx(k)));\r
+                    if sim_coh_dpsk\r
+                      prev_sym_rx(k) = qpsk_mod(qpsk_demod(tmp));\r
+                    else\r
+                      prev_sym_rx(k) = tmp;\r
+                    end\r
+                    s_ch(i+k-1) = rx_symb;\r
                 end\r
-                rx_bits((2*(i-1)+1):(2*i)) = qpsk_demod(rx_symb);\r
+                rx_bits((2*(i-1+k-1)+1):(2*(i+k-1))) = qpsk_demod(rx_symb);\r
                 rx_symb_log = [rx_symb_log rx_symb];\r
+              end\r
             end\r
 \r
 if newldpc\r
@@ -347,6 +364,7 @@ if newldpc
               %rx_bits_buf(st_rx_bits:st_rx_bits+20-1)\r
               error_positions = xor(rx_bits_buf(st_rx_bits:en_rx_bits), tx_bits_buf(1:framesize));\r
               Nerrs = sum(error_positions);\r
+              sim_out.errors_log = [sim_out.errors_log error_positions];\r
               Terrs += Nerrs;\r
               Tbits += length(tx_bits);\r
 \r
@@ -359,6 +377,7 @@ if newldpc
                 %detected_data = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, s_ch, min(100,EsNo), hf_fading);\r
                 %error_positions = xor( detected_data(1:framesize*rate), tx_bits(1:framesize*rate) );\r
                 Nerrs = sum(error_positions);\r
+                sim_out.ldpc_errors_log = [sim_out.ldpc_errors_log error_positions];\r
                 if Nerrs\r
                     Ferrsldpc++;\r
                 end\r
@@ -384,6 +403,7 @@ else
                 end\r
                 Terrsldpc += Nerrs;\r
                 Tbitsldpc += framesize*rate;\r
+\r
             end\r
         end\r
 end\r
@@ -413,6 +433,13 @@ end
     end\r
     \r
     Ebvec = Esvec - 10*log10(bps);\r
+\r
+    % account for extra power rqd for pilot symbols\r
+\r
+    if (phase_est_method == 2) && (phase_est)\r
+      Ebvec += 10*log10(Ns/(Ns-1));\r
+    end\r
+\r
     sim_out.BERvec          = BERvec;\r
     sim_out.Ebvec           = Ebvec;\r
     sim_out.TERvec          = TERvec;\r
@@ -486,6 +513,8 @@ if 0
 end\r
      end\r
 \r
+size(sim_out.errors_log)\r
+\r
 endfunction\r
 \r
 % Gray coded QPSK modulation function\r
@@ -648,20 +677,21 @@ endfunction
 function phase_est_hf\r
   sim_in = standard_init();\r
 \r
-  sim_in.Rs               = 200;\r
-  sim_in.Nc               = 4;\r
+  sim_in.Rs               = 100;\r
+  sim_in.Nc               = 8;\r
 \r
   sim_in.verbose          = 1;\r
   sim_in.plot_scatter     = 0;\r
 \r
-  sim_in.Esvec            = 2:8; \r
-  sim_in.Ntrials          = 30;\r
+  sim_in.Esvec            = 5:15; \r
+  sim_in.Ntrials          = 100;\r
 \r
   sim_in.newldpc          = 1;\r
   sim_in.ldpc_code_rate   = 1/2;\r
   sim_in.ldpc_code        = 1;\r
 \r
   sim_in.phase_est        = 0;\r
+  sim_in.sim_coh_dpsk     = 0;\r
   sim_in.phase_est_method = 2;\r
   sim_in.Np               = 3;\r
   sim_in.phase_offset     = 0;\r
@@ -674,47 +704,40 @@ function phase_est_hf
 \r
   baseline                = ber_test(sim_in, 'qpsk');\r
 \r
-  sim_in.hf_mag_only      = 1;\r
+  sim_in.hf_mag_only      = 0;\r
   sim_in.phase_est_method = 2;\r
   sim_in.phase_est        = 1;\r
+  sim_in.Np               = 3;\r
+  pilot_3                 = ber_test(sim_in, 'qpsk');\r
+  sim_in.Np               = 5;\r
+  pilot_5                 = ber_test(sim_in, 'qpsk');\r
   sim_in.Np               = 7;\r
   pilot_7                 = ber_test(sim_in, 'qpsk');\r
 \r
-  sim_in.hf_mag_only      = 0;\r
-  pilot_7b                = ber_test(sim_in, 'qpsk');\r
-\r
-if 0\r
+if 1\r
   sim_in.phase_est        = 0;\r
   dqpsk                   = ber_test(sim_in, 'dqpsk');\r
 \r
-  sim_in.phase_est        = 1;\r
-  sim_in.phase_est_method = 3;\r
-  sim_in.Np               = 41;\r
-  dqpsk_strip_41          = ber_test(sim_in, 'dqpsk');\r
-end\r
-  \r
   figure(1); \r
   clf;\r
   semilogy(baseline.Ebvec, baseline.BERvec,'r;QPSK CCIR poor;')\r
   hold on;\r
   semilogy(baseline.Ebvec, baseline.BERldpcvec,'r;QPSK CCIR poor ldpc;')\r
-  semilogy(pilot_7.Ebvec, pilot_7.BERvec,'b;QPSK CCIR poor ldpc pilot 7;')\r
-  semilogy(pilot_7.Ebvec, pilot_7.BERldpcvec,'b;QPSK CCIR poor ldpc pilot 7;')\r
-  semilogy(pilot_7b.Ebvec, pilot_7b.BERvec,'g;QPSK CCIR poor ldpc pilot 7b;')\r
-  semilogy(pilot_7b.Ebvec, pilot_7b.BERldpcvec,'g;QPSK CCIR poor ldpc pilot 7b;')\r
-if 0\r
-  semilogy(dqpsk.Ebvec, dqpsk.BERvec,'c;DQPSK CCIR poor ldpc;')\r
-  semilogy(dqpsk.Ebvec, dqpsk.BERldpcvec,'c;DQPSK CCIR poor ldpc;')\r
-  semilogy(dqpsk_strip_41.Ebvec, dqpsk_strip_41.BERvec,'m;DQPSK CCIR poor ldpc strip 41;')\r
-  semilogy(dqpsk_strip_41.Ebvec, dqpsk_strip_41.BERldpcvec,'m;DQPSK CCIR poor ldpc strip 41;')\r
-end\r
+  semilogy(pilot_3.Ebvec, pilot_3.BERvec,'b;QPSK CCIR poor ldpc pilot 3;')\r
+  semilogy(pilot_3.Ebvec, pilot_3.BERldpcvec,'b;QPSK CCIR poor ldpc pilot 3;')\r
+  semilogy(pilot_5.Ebvec, pilot_5.BERvec,'g;QPSK CCIR poor ldpc pilot 5;')\r
+  semilogy(pilot_5.Ebvec, pilot_5.BERldpcvec,'g;QPSK CCIR poor ldpc pilot 5;')\r
+  semilogy(pilot_7.Ebvec, pilot_7.BERvec,'m;QPSK CCIR poor ldpc pilot 7;')\r
+  semilogy(pilot_7.Ebvec, pilot_7.BERldpcvec,'m;QPSK CCIR poor ldpc pilot 7;')\r
+  semilogy(dqpsk.Ebvec, dqpsk.BERvec,'k;DQPSK CCIR poor ldpc;')\r
+  semilogy(dqpsk.Ebvec, dqpsk.BERldpcvec,'k;DQPSK CCIR poor ldpc;')\r
 \r
   hold off;\r
   xlabel('Eb/N0')\r
   ylabel('BER')\r
   grid("minor")\r
   axis([min(Ebvec) max(Ebvec) 1E-2 1])\r
-\r
+end\r
 endfunction\r
 \r
 function phase_est_awgn\r
@@ -788,10 +811,145 @@ function phase_est_awgn
   axis([min(Ebvec) max(Ebvec) 1E-2 1])\r
 endfunction\r
 \r
+function test_dpsk\r
+  sim_in = standard_init();\r
+\r
+  sim_in.Rs               = 100;\r
+  sim_in.Nc               = 8;\r
+\r
+  sim_in.verbose          = 1;\r
+  sim_in.plot_scatter     = 0;\r
+\r
+  sim_in.Esvec            = 5; \r
+  sim_in.Ntrials          = 30;\r
+\r
+  sim_in.newldpc          = 1;\r
+  sim_in.ldpc_code_rate   = 1/2;\r
+  sim_in.ldpc_code        = 1;\r
+\r
+  sim_in.phase_est        = 0;\r
+  sim_in.phase_est_method = 3;\r
+  sim_in.Np               = 41;\r
+  sim_in.phase_offset     = 0;\r
+  sim_in.w_offset         = 0;\r
+  sim_in.sim_coh_dpsk     = 0;\r
+\r
+  sim_in.hf_sim           = 0;\r
+  sim_in.hf_mag_only      = 1;\r
+\r
+  Ebvec = sim_in.Esvec - 10*log10(2);\r
+\r
+  baseline                = ber_test(sim_in, 'qpsk');\r
+  sim_in.phase_est        = 0;\r
+  dqpsk                   = ber_test(sim_in, 'dqpsk');\r
+\r
+  sim_in.phase_est        = 1;\r
+  sim_in.phase_est_method = 3;\r
+  sim_in.sim_coh_dpsk     = 1;\r
+  sim_in.Np               = 41;\r
+  dqpsk_strip_41          = ber_test(sim_in, 'dqpsk');\r
+  \r
+  figure(1); \r
+  clf;\r
+  semilogy(baseline.Ebvec, baseline.BERvec,'r;QPSK CCIR poor;')\r
+  hold on;\r
+  semilogy(baseline.Ebvec, baseline.BERldpcvec,'r;QPSK CCIR poor ldpc;')\r
+  semilogy(dqpsk.Ebvec, dqpsk.BERvec,'c;DQPSK CCIR poor ldpc;')\r
+  semilogy(dqpsk.Ebvec, dqpsk.BERldpcvec,'c;DQPSK CCIR poor ldpc;')\r
+  semilogy(dqpsk_strip_41.Ebvec, dqpsk_strip_41.BERvec,'m;DQPSK CCIR poor ldpc strip 41;')\r
+  semilogy(dqpsk_strip_41.Ebvec, dqpsk_strip_41.BERldpcvec,'m;DQPSK CCIR poor ldpc strip 41;')\r
+\r
+  hold off;\r
+  xlabel('Eb/N0')\r
+  ylabel('BER')\r
+  grid("minor")\r
+  axis([min(Ebvec) max(Ebvec) 1E-2 1])\r
+\r
+endfunction\r
+\r
+function gen_error_pattern_qpsk()\r
+  sim_in = standard_init();\r
+\r
+  % model codec and uncoded streams as 1000 bit/s each\r
+\r
+  sim_in.Rs               = 100;\r
+  sim_in.Nc               = 4;\r
+\r
+  sim_in.verbose          = 1;\r
+  sim_in.plot_scatter     = 0;\r
+\r
+  sim_in.Esvec            = 10; % Eb/No=2dB\r
+  sim_in.Ntrials          = 30;\r
+\r
+  sim_in.newldpc          = 1;\r
+  sim_in.ldpc_code_rate   = 1/2;\r
+  sim_in.ldpc_code        = 1;\r
+\r
+  sim_in.phase_est        = 1;\r
+  sim_in.phase_est_method = 2;\r
+  sim_in.Np               = 5;\r
+  sim_in.phase_offset     = 0;\r
+  sim_in.w_offset         = 0;\r
+  sim_in.sim_coh_dpsk     = 0;\r
+\r
+  sim_in.hf_sim           = 1;\r
+  sim_in.hf_mag_only      = 0;\r
+\r
+  qpsk                    = ber_test(sim_in, 'qpsk');\r
+ \r
+  length(qpsk.errors_log) \r
+  length(qpsk.ldpc_errors_log)\r
+  % multiplex errors into prot and unprot halves of 52 bit codec frames\r
+\r
+  error_pattern = [];\r
+  for i=1:26:length(qpsk.ldpc_errors_log)-52\r
+    error_pattern = [error_pattern qpsk.ldpc_errors_log(i:i+25) qpsk.errors_log(i:i+25)  zeros(1,4)];\r
+    %error_pattern = [error_pattern qpsk.ldpc_errors_log(i:i+25) zeros(1,26)  zeros(1,4)];\r
+    %error_pattern = [error_pattern zeros(1,26) qpsk.errors_log(i:i+25)  zeros(1,4)];\r
+  end\r
+\r
+  fep=fopen("qpsk_errors_2dB.bin","wb"); fwrite(fep, error_pattern, "short"); fclose(fep);\r
+\r
+endfunction\r
+\r
+function gen_error_pattern_dpsk()\r
+  sim_in = standard_init();\r
+\r
+  sim_in.Rs               = 50;\r
+  sim_in.Nc               = 16;\r
+\r
+  sim_in.verbose          = 1;\r
+  sim_in.plot_scatter     = 1;\r
+\r
+  sim_in.Esvec            = 10; % Eb/No=Es/No-3\r
+  sim_in.Ntrials          = 30;\r
+\r
+  sim_in.newldpc          = 1;\r
+  sim_in.ldpc_code_rate   = 1/2;\r
+  sim_in.ldpc_code        = 0;\r
+\r
+  sim_in.phase_est        = 0;\r
+  sim_in.phase_est_method = 3;\r
+  sim_in.Np               = 41;\r
+  sim_in.phase_offset     = 0;\r
+  sim_in.w_offset         = 0;\r
+  sim_in.sim_coh_dpsk     = 0;\r
+\r
+  sim_in.hf_sim           = 1;\r
+  sim_in.hf_mag_only      = 1;\r
+\r
+  dqpsk                    = ber_test(sim_in, 'dqpsk');\r
+ \r
+  fep=fopen("dqpsk_errors_12dB.bin","wb"); fwrite(fep, dqpsk.errors_log, "short"); fclose(fep);\r
+\r
+endfunction\r
+\r
 % Start simulations ---------------------------------------\r
 \r
 more off;\r
 \r
 %ideal();\r
-phase_est_hf();\r
+%phase_est_hf();\r
 %phase_est_awgn();\r
+%test_dpsk();\r
+gen_error_pattern_qpsk\r