plotting SNR contributions from various sources

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

diff --git a/codec2-dev/octave/test_dqpsk.m b/codec2-dev/octave/test_dqpsk.m
index 713ae3964ea9f6f08a15d28287fa458dd7068337..3d33bc376aab1e474cb580fddf3452e48c1dfca5 100644 (file)
--- a/codec2-dev/octave/test_dqpsk.m
+++ b/codec2-dev/octave/test_dqpsk.m
@@ -24,6 +24,7 @@ function sim_out = ber_test(sim_in)
     hf_phase_only    = sim_in.hf_phase_only;\r
     hf_mag_only      = sim_in.hf_mag_only;\r
     Nc               = sim_in.Nc;\r
+    symbol_amp       = sim_in.symbol_amp;\r
 \r
     bps              = 2;\r
     Nsymb            = framesize/bps;\r
@@ -88,10 +89,11 @@ function sim_out = ber_test(sim_in)
         \r
         Terrs = 0;  Tbits = 0;\r
 \r
-        tx_symb_log      = [];\r
-        rx_symb_log      = [];\r
-        noise_log        = [];\r
-        \r
+        tx_symb_log        = [];\r
+        rx_symb_log        = [];\r
+        noise_log          = [];\r
+        sim_out.errors_log = [];\r
+\r
         % init HF channel\r
 \r
         hf_n = 1;\r
@@ -100,7 +102,11 @@ function sim_out = ber_test(sim_in)
         hf_model = ones(Ntrials*Nsymb/Nc, Nc); % defaults for plotting surface\r
 \r
         sim_out.errors_log = [];\r
-        sim_out.ldpc_errors_log = [];\r
+        sim_out.Nerrs = [];\r
+        sim_out.snr_log = [];\r
+        sim_out.hf_model_pwr = [];\r
+\r
+        symbol_amp_index = 1;\r
 \r
         for nn = 1: Ntrials\r
                   \r
@@ -114,10 +120,11 @@ function sim_out = ber_test(sim_in)
                 tx_symb = qpsk_mod(tx_bits(2*(i-1+k-1)+1:2*(i+k-1)));\r
                 tx_symb *= prev_sym_tx(k);\r
                 prev_sym_tx(k) = tx_symb;\r
-                s(i+k-1) = tx_symb;\r
+                s(i+k-1) = symbol_amp(symbol_amp_index)*tx_symb;\r
               end\r
             end\r
             s_ch = s;\r
+            symbol_amp_index++;\r
 \r
             % HF channel simulation  ------------------------------------\r
             \r
@@ -156,6 +163,12 @@ function sim_out = ber_test(sim_in)
             \r
             tx_symb_log = [tx_symb_log s_ch];\r
 \r
+            % "genie" SNR estimate \r
+            \r
+            snr = (s_ch*s_ch')/(Nsymb*variance);\r
+            sim_out.snr_log = [sim_out.snr_log snr];\r
+            sim_out.hf_model_pwr = [sim_out.hf_model_pwr mean(hf_fading.^2)];\r
+\r
             % AWGN noise and phase/freq offset channel simulation\r
             % 0.5 factor ensures var(noise) == variance , i.e. splits power between Re & Im\r
 \r
@@ -186,9 +199,12 @@ function sim_out = ber_test(sim_in)
             end\r
 \r
             error_positions = xor(rx_bits, tx_bits);\r
+            sim_out.errors_log = [sim_out.errors_log error_positions];\r
             Nerrs = sum(error_positions);\r
+            sim_out.Nerrs = [sim_out.Nerrs Nerrs];\r
             Terrs += Nerrs;\r
             Tbits += length(tx_bits);\r
+\r
         end\r
 \r
         TERvec(ne) = Terrs;\r
@@ -228,6 +244,7 @@ function sim_out = ber_test(sim_in)
         grid\r
         title('HF Channel Es/No');\r
 \r
+        if 0 \r
         figure(4);\r
         clf;\r
         subplot(211)\r
@@ -236,6 +253,7 @@ function sim_out = ber_test(sim_in)
         subplot(212)\r
         plot(y,angle(hf_model(y,1)))\r
         title('HF Channel Carrier 1 Phase');\r
+        end\r
     end\r
 \r
 endfunction\r
@@ -284,23 +302,87 @@ function sim_in = standard_init
   sim_in.hf_mag_only      = 0;\r
 endfunction\r
 \r
+function awgn_hf_ber_curves()\r
+  sim_in = standard_init();\r
+\r
+  Ebvec = sim_in.Esvec - 10*log10(2);\r
+  BER_theory = 0.5*erfc(sqrt(10.^(Ebvec/10)));\r
+\r
+  dpsk_awgn = ber_test(sim_in);\r
+  sim_in.hf_sim           = 1;\r
+  dpsk_hf   = ber_test(sim_in);\r
+\r
+  figure(1); \r
+  clf;\r
+  semilogy(Ebvec, BER_theory,'r;QPSK theory;')\r
+  hold on;\r
+  semilogy(dpsk_awgn.Ebvec, dpsk_awgn.BERvec,'g;DQPSK;')\r
+  semilogy(dpsk_hf.Ebvec, dpsk_hf.BERvec,'g;DQPSK HF;')\r
+  hold off;\r
+  xlabel('Eb/N0')\r
+  ylabel('BER')\r
+  grid("minor")\r
+  axis([min(Ebvec) max(Ebvec) 1E-3 1])\r
+end\r
+\r
 sim_in = standard_init();\r
 \r
-Ebvec = sim_in.Esvec - 10*log10(2);\r
-BER_theory = 0.5*erfc(sqrt(10.^(Ebvec/10)));\r
+% energy file sampled every 10ms\r
+\r
+load ../src/ve9qrp.txt\r
+pdB=10*log10(ve9qrp);\r
+for i=1:length(pdB)\r
+  if pdB(i) < 0\r
+    pdB(i) = 0;\r
+  end\r
+end\r
+\r
+% Down sample to 40ms rate used for 1300 bit/s codec, every 4th sample is transmitted\r
+\r
+pdB = pdB(4:4:length(pdB));\r
+\r
+% Use linear mapping function in dB domain to map to symbol power\r
 \r
-dpsk_awgn = ber_test(sim_in);\r
-sim_in.hf_sim           = 1;\r
-dpsk_hf   = ber_test(sim_in);\r
+power_map_x  = [ 0 20 24 40 50 ];\r
+power_map_y  = [-3 -3  0 6  6];\r
+mapped_pdB = interp1(power_map_x, power_map_y, pdB);\r
 \r
-figure(1); \r
+sim_in.symbol_amp = 10 .^ (mapped_pdB/20);\r
+%sim_in.symbol_amp = ones(1,length(pdB));\r
+sim_in.plot_scatter = 1;\r
+sim_in.verbose      = 2;\r
+sim_in.hf_sim       = 1;\r
+sim_in.Esvec        = 10;\r
+\r
+dqpsk_pwr_hf = ber_test(sim_in);\r
+\r
+% note: need way to test that power is aligned with speech\r
+\r
+figure(4)\r
 clf;\r
-semilogy(Ebvec, BER_theory,'r;QPSK theory;')\r
+plot((1:sim_in.Ntrials)*80*4, pdB(1:sim_in.Ntrials));\r
 hold on;\r
-semilogy(dpsk_awgn.Ebvec, dpsk_awgn.BERvec,'g;DQPSK;')\r
-semilogy(dpsk_hf.Ebvec, dpsk_hf.BERvec,'g;DQPSK HF;')\r
+plot((1:sim_in.Ntrials)*80*4, mapped_pdB(1:sim_in.Ntrials),'r');\r
 hold off;\r
-xlabel('Eb/N0')\r
-ylabel('BER')\r
-grid("minor")\r
-axis([min(Ebvec) max(Ebvec) 1E-3 1])\r
+\r
+figure(5)\r
+clf;\r
+\r
+s = load_raw("../raw/ve9qrp.raw");\r
+M=320; M_on_2 = M/2; % processing delay between input speech and centre of analysis window\r
+plot(M_on_2:(M_on_2-1+sim_in.Ntrials*M),s(1:sim_in.Ntrials*M))\r
+hold on;\r
+plot((1:sim_in.Ntrials)*M, 5000*sim_in.symbol_amp(1:sim_in.Ntrials),'r');\r
+hold off;\r
+axis([1 sim_in.Ntrials*M -3E4 3E4]);\r
+\r
+figure(6)\r
+clf;\r
+plot((1:sim_in.Ntrials)*M, 20*log10(sim_in.symbol_amp(1:sim_in.Ntrials)),'b;Es (dB);');\r
+hold on;\r
+plot((1:sim_in.Ntrials)*M, 10*log10(dqpsk_pwr_hf.hf_model_pwr),'g;Fading (dB);');\r
+plot((1:sim_in.Ntrials)*M, 10*log10(dqpsk_pwr_hf.snr_log),'r;Es/No (dB);');\r
+hold off;\r
+\r
+fep=fopen("dqpsk_errors_pwr.bin","wb"); fwrite(fep, dqpsk_pwr_hf.errors_log, "short"); fclose(fep);\r
+fmute=fopen("dqpsk_mute.bin","wb"); fwrite(fmute, mute, "short"); fclose(fmute);\r