HF channel model and carrier power normalisation, but dqpsk broken

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

diff --git a/codec2-dev/octave/test_qpsk.m b/codec2-dev/octave/test_qpsk.m
index c84c10153c47c8702b0cfddb071f59b752224ced..fdb22d34324268dd7986cb88b232d5066f261e09 100644 (file)
--- a/codec2-dev/octave/test_qpsk.m
+++ b/codec2-dev/octave/test_qpsk.m
@@ -1,17 +1,20 @@
 % test_qpsk.m\r
 % David Rowe Feb 2014\r
 %\r
-% Single sample per symbol QPSK modem simulation, based on code by Bill Cowley\r
-% Generates curves BER versus E/No curves for different modems.  Design to\r
-% test perform initial tests on coherent demodulation for HF channels without\r
-% building a full blown modem.  Lacks filtering, timing estimation, frame sync.\r
-% \r
+% QPSK modem simulation, initially based on code by Bill Cowley\r
+% Generates curves BER versus E/No curves for different modems.\r
+% Design to test coherent demodulation ideas on HF channels without\r
+% building a full blown modem.  Lacks timing estimation, frame sync.\r
+\r
+% TODO\r
+%   [ ] put spreading sample files sonewhere useful\r
 \r
 1;\r
 \r
 % main test function \r
 \r
 function sim_out = ber_test(sim_in, modulation)\r
+    Fs = 8000;\r
 \r
     framesize        = sim_in.framesize;\r
     Ntrials          = sim_in.Ntrials;\r
@@ -22,8 +25,7 @@ function sim_out = ber_test(sim_in, modulation)
     plot_scatter     = sim_in.plot_scatter;\r
     Rs               = sim_in.Rs;\r
     hf_sim           = sim_in.hf_sim;\r
-    hf_spread_hz     = sim_in.hf_spread_hz;\r
-    hf_delay_samples = sim_in.hf_delay_samples;\r
+    Nhfdelay         = floor(sim_in.hf_delay_ms*1000/Fs);\r
 \r
     bps              = 2;\r
     Nsymb            = framesize/bps;\r
@@ -34,26 +36,32 @@ function sim_out = ber_test(sim_in, modulation)
     Np               = 5;\r
     r_delay_line     = zeros(1,Np);\r
     s_delay_line     = zeros(1,Np);\r
-    hf_delay_line    = zeros(1,Nsymb+hf_delay_samples);\r
     spread_main_phi  = 0;\r
     spread_delay_phi = 0;\r
     spread_main_phi_log = [];\r
 \r
-    % convert "spreading" samples from 1kHz carrier at Fs to complex baseband at Rs\r
+    % convert "spreading" samples from 1kHz carrier at Fs to complex baseband\r
 \r
-    Fs = 8000; Fc = 1000;\r
+    Fc = 1000;\r
     fspread = fopen("../unittest/sine1k_2Hz_spread.raw","rb");\r
     spread1k = fread(fspread, "int16")/10000;\r
+    fclose(fspread);\r
+    fspread = fopen("../unittest/sine1k_2ms_delay_2Hz_spread.raw","rb");\r
+    spread1k_2ms = fread(fspread, "int16")/10000;\r
+    fclose(fspread);\r
 \r
     % down convert to complex baseband\r
     spreadbb = spread1k.*exp(-j*(2*pi*Fc/Fs)*(1:length(spread1k))');\r
+    spreadbb_2ms = spread1k_2ms.*exp(-j*(2*pi*Fc/Fs)*(1:length(spread1k_2ms))');\r
 \r
     % remove -2000 Hz image\r
     b = fir1(50, 5/Fs);\r
-    spreadlpf = filter(b,1,spreadbb);\r
+    spread = filter(b,1,spreadbb);\r
+    spread_2ms = filter(b,1,spreadbb_2ms);\r
+\r
+    % Determine "gain" of HF channel model, so we can normalise later\r
 \r
-    % decimate to symbol rate\r
-    spread = spreadlpf(1:floor(Fs/Rs):length(spreadlpf));\r
+    hf_gain = 1.0/sqrt(var(spread)+var(spread_2ms));\r
 \r
     sc = 1;\r
 \r
@@ -69,8 +77,7 @@ function sim_out = ber_test(sim_in, modulation)
     tx_filter_memory = zeros(1, Nfilter);\r
     rx_filter_memory = zeros(1, Nfilter);\r
     s_delay_line_filt = zeros(1,Nfiltsym);\r
- \r
-    wc = 2*pi*1500/Fs;\r
+    hf_sim_delay_line = zeros(1,M+Nhfdelay);\r
 \r
     for ne = 1:length(Esvec)\r
         Es = Esvec(ne);\r
@@ -97,14 +104,16 @@ function sim_out = ber_test(sim_in, modulation)
         rx_baseband_log = [];\r
         tx_baseband_log = [];\r
         noise_log = [];\r
+        c_est_log = [];\r
+        c_est = 0;\r
 \r
         tx_phase = rx_phase = 0;\r
 \r
         for nn = 1: Ntrials\r
                   \r
-           tx_bits = round( rand( 1, framesize ) );\r
+            tx_bits = round( rand( 1, framesize ) );\r
 \r
-           % modulate\r
+            % modulate\r
 \r
             s = zeros(1, Nsymb);\r
             for i=1:Nsymb\r
@@ -138,10 +147,27 @@ function sim_out = ber_test(sim_in, modulation)
                end\r
                tx_filter_memory(1:Nfilter-M) = tx_filter_memory(M+1:Nfilter);\r
                tx_filter_memory(Nfilter-M+1:Nfilter) = zeros(1,M);\r
-               tx_filt_log = [tx_filt_log tx_filt];\r
                \r
+               % HF channel simulation\r
+\r
+               if hf_sim\r
+                   hf_sim_delay_line(1:Nhfdelay) = hf_sim_delay_line(M+1:M+Nhfdelay);\r
+                   hf_sim_delay_line(Nhfdelay+1:M+Nhfdelay) = tx_filt;\r
+\r
+                   if (sc > (length(spread)-M))\r
+                       sc =1 ;\r
+                   end\r
+                   tx_filt = tx_filt.*spread(sc:sc+M-1)' + hf_sim_delay_line(1:M).*spread_2ms(sc+sc+M-1);\r
+                   sc += M;\r
+\r
+                   % normalise so average HF power C=1\r
+\r
+                   tx_filt *= hf_gain;\r
+               end\r
+               tx_filt_log = [tx_filt_log tx_filt];\r
+\r
                % AWGN noise and phase/freq offset channel simulation\r
-               % 0.5 factor ensures var(noise) == variance , i.e. splits power betwen Re & Im\r
+               % 0.5 factor ensures var(noise) == variance , i.e. splits power between Re & Im\r
 \r
                noise = sqrt(variance*0.5)*( randn(1,M) + j*randn(1,M) );\r
                noise_log = [noise_log noise];\r
@@ -165,13 +191,6 @@ function sim_out = ber_test(sim_in, modulation)
                s_ch(k) = rx_filt;               \r
             end\r
 \r
-            % Channel simulation\r
-\r
-            if hf_sim\r
-                s_ch = s_ch.*spread(sc:sc+Nsymb-1)';\r
-                sc += Nsymb;\r
-            end\r
-\r
             % coherent demod phase estimation and correction\r
             % need sliding window\r
             % use known (pilot) symbols\r
@@ -249,6 +268,8 @@ function sim_out = ber_test(sim_in, modulation)
         clf;\r
         scat = rx_symb_log(2*Nfiltsym:length(rx_symb_log)) .* exp(j*pi/4);\r
         plot(real(scat), imag(scat),'+');\r
+        figure(3);\r
+        plot(c_est_log);\r
     end\r
 endfunction\r
 \r
@@ -278,28 +299,22 @@ sim_in.Esvec            = 1:10;
 sim_in.Ntrials          = 100;\r
 sim_in.framesize        = 30;\r
 sim_in.phase_offset     = 0;\r
-sim_in.phase_est        = 0;\r
+sim_in.phase_est        = 1;\r
 sim_in.w_offset         = 0;\r
 sim_in.plot_scatter     = 0;\r
 sim_in.Rs               = 100;\r
-sim_in.hf_sim           = 0;\r
-sim_in.hf_spread_hz     = 2;\r
-sim_in.hf_delay_samples = 5;\r
+sim_in.hf_sim           = 1;\r
+sim_in.hf_delay_ms      = 2;\r
 \r
 sim_qpsk                = ber_test(sim_in, 'qpsk');\r
 \r
 sim_in.phase_offset     = 0;\r
 sim_in.phase_est        = 0;\r
 sim_in.w_offset         = 0;  \r
-%sim_qpsk_coh            = ber_test(sim_in, 'qpsk');\r
-\r
-sim_in.phase_offset     = 0;\r
-sim_in.phase_est        = 1;\r
-sim_in.w_offset         = 0;  \r
 sim_in.plot_scatter     = 1;\r
-sim_in.Esvec            = 7;\r
+sim_in.Esvec            = 70;\r
 sim_in.hf_sim           = 0;\r
-sim_qpsk_scatter        = ber_test(sim_in, 'qpsk');\r
+sim_qpsk_scatter        = ber_test(sim_in, 'dqpsk');\r
 \r
 figure(1); \r
 clf;\r