refactored a little and wrote Octave function to generate files of modulated bits

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

diff --git a/codec2-dev/octave/tqpsk.m b/codec2-dev/octave/tqpsk.m
index 4ec8765a0e8ba97b5eb717f3ba34411ab56f786d..f6557c91a871885a639be8f9aa36798ac9f4918d 100644 (file)
--- a/codec2-dev/octave/tqpsk.m
+++ b/codec2-dev/octave/tqpsk.m
@@ -2,7 +2,16 @@
 %
 % David Rowe Sep 2015
 %
-% QPSK modem Octave test sctipt
+% QPSK modem Octave test scripts
+
+% TODO:
+%   correct freq offset
+%   output tx frame to a file of bits for 3rd party modulator
+%   improve phase est impl loss using longer window
+%   separate into mod and demod functions
+%   generate filter coeffs, save to C header file
+%   fine timing estimator and meas IL
+%   coarse timing (frame sync)
 
 % make sure we get same random numbers every time
 
@@ -10,16 +19,8 @@ more off;
 qpsk;
 ldpc;
 
-function sim_out = run_simulation(sim_in)
-  Ntrials      = sim_in.Ntrials;
-  EsNodBvec    = sim_in.EsNodBvec;
-  verbose      = sim_in.verbose;
-  phase_offset = sim_in.phase_offset;
-  phase_est_en = sim_in.phase_est_en;
-  verbose      = sim_in.verbose;
 
-  rand('state',1); 
-  randn('state',1);
+function [aqpsk code_param] = init_modem_ldpc
 
   % Init modem ------------------------------------------------------------------------------
 
@@ -45,9 +46,6 @@ function sim_out = run_simulation(sim_in)
   aqpsk.Np            = 2;                                  % number of pilots used for phase estimation
   aqpsk.rx_pilot_buf  = zeros(1,aqpsk.Np);
 
-  Nsymb = aqpsk.Nsymb;
-  Nsymbpilot = aqpsk.Nsymbpilot;
-
   % Init LDPC --------------------------------------------------------------------------
 
   currentdir = pwd;
@@ -66,6 +64,27 @@ function sim_out = run_simulation(sim_in)
   code_param = ldpc_init(aqpsk.code_rate, aqpsk.framesize, modulation, mod_order, mapping);
   assert(code_param.data_bits_per_frame == aqpsk.Ndatabits);
 
+endfunction
+
+
+function sim_out = run_simulation(sim_in)
+  Ntrials      = sim_in.Ntrials;
+  EsNodBvec    = sim_in.EsNodBvec;
+  verbose      = sim_in.verbose;
+  phase_offset = sim_in.phase_offset;
+  phase_est_en = sim_in.phase_est_en;
+  verbose      = sim_in.verbose;
+
+  rand('state',1); 
+  randn('state',1);
+
+  % Init modem and LPDC ----------------------------------------------------------------
+
+  [aqpsk code_param] = init_modem_ldpc;
+
+  Nsymb = aqpsk.Nsymb;
+  Nsymbpilot = aqpsk.Nsymbpilot;
+
   % Start simulation -------------------------------------------------------------------
 
   for ne = 1:length(EsNodBvec)
@@ -135,32 +154,77 @@ end
 
 % ----------------------------------------------------------------------
 
-sim_in.phase_est_en = 1;
-sim_in.phase_offset = 0;
-sim_in.verbose      = 1;
-sim_in.EsNodBvec    = -2:10;
-sim_in.Ntrials      = 10;
-
-sim_out = run_simulation(sim_in);
-
-EbNodB = sim_in.EsNodBvec - 10*log10(2);  % QPSK has two bits/symbol
-uncoded_BER_theory = 0.5*erfc(sqrt(10.^(EbNodB/10)));
-
-figure(1)
-clf
-semilogy(EbNodB, uncoded_BER_theory,'r;uncoded QPSK theory;')
-hold on;
-semilogy(EbNodB, sim_out.BER_uc,'r+;uncoded QPSK simulated;')
-semilogy(EbNodB-10*log10(sim_out.code_rate), sim_out.BER+1E-10,'g;LDPC coded QPSK simulated;');
-hold off;
-grid('minor')
-xlabel('Eb/No (dB)')
-ylabel('BER')
-axis([min(EbNodB) max(EbNodB) min(uncoded_BER_theory) 1])
-
-% correct freq offset
-% output "demodulated" tx frame to a file of bits for 3rd party modulator
-% actually estimate Es/No for LDPC
-% improve phase est impl loss
-% separate into enc and demod functions
-% generate filter coeffs, save to C header file
+
+function symbol_rate_curves
+  sim_in.phase_est_en = 1;
+  sim_in.phase_offset = 0;
+  sim_in.verbose      = 1;
+  sim_in.EsNodBvec    = -2:10;
+  sim_in.Ntrials      = 10;
+
+  sim_out = run_simulation(sim_in);
+
+  EbNodB = sim_in.EsNodBvec - 10*log10(2);  % QPSK has two bits/symbol
+  uncoded_BER_theory = 0.5*erfc(sqrt(10.^(EbNodB/10)));
+
+  figure(1)
+  clf
+  semilogy(EbNodB, uncoded_BER_theory,'r;uncoded QPSK theory;')
+  hold on;
+  semilogy(EbNodB, sim_out.BER_uc,'r+;uncoded QPSK simulated;')
+  semilogy(EbNodB-10*log10(sim_out.code_rate), sim_out.BER+1E-10,'g;LDPC coded QPSK simulated;');
+  hold off;
+  grid('minor')
+  xlabel('Eb/No (dB)')
+  ylabel('BER')
+  axis([min(EbNodB) max(EbNodB) min(uncoded_BER_theory) 1])
+
+endfunction
+
+
+% We are using a 3rd party modulator, this function generates bits
+% that can be passed to it.  It's organised as a file of 8 bit bytes
+% MSB is sent first.  Each pair of bits is mapped to phases as:
+%
+% MSB LSB   Phase
+%   0   0       0
+%   0   1      90
+%   1   1     180
+%   1   0     270 
+
+function generate_modulator_input_file(filename, seconds)
+         
+  [aqpsk code_param] = init_modem_ldpc;
+
+  Nframes = floor(seconds/aqpsk.Rf);
+  Nsecs = Nframes*aqpsk.Rf;
+  Nbitsframe = aqpsk.Nsymbpilot*aqpsk.bits_per_symb;
+  Nbytesframe = Nbitsframe/8;
+
+  printf("Frames: %d Seconds: %3.1f QPSK Symb Rate: %5.2f\n", Nframes, Nsecs, aqpsk.Rs);
+  printf("Bit/frame: %d  Packed Bytes/frame: %4.1f\n", Nbitsframe, Nbytesframe);
+
+  fbit = fopen(filename,"wb"); 
+ 
+  for nn = 1:Nframes
+    tx_bits = round(rand(1,  aqpsk.Ndatabits));
+    [tx_codeword, tx_coded_symb] = ldpc_enc(tx_bits, code_param);
+    symbpilot_tx = insert_pilots(tx_coded_symb, aqpsk.pilots, aqpsk.Npilotstep);
+
+    % demodulate QPSK symbols back to bits
+   
+    tx_frame = zeros(1, aqpsk.Nsymbpilot*aqpsk.bits_per_symb);
+    for s=1:aqpsk.Nsymbpilot
+      tx_frame(2*(s-1)+1:2*s) = qpsk_demod(symbpilot_tx(s));
+    end
+    packed_tx_frame = packmsb(tx_frame);
+    fwrite(fbit, packed_tx_frame, "uchar");
+  end
+  
+  fclose(fbit);
+
+endfunction
+
+% ---------------------------------------------------------------------------------
+
+generate_modulator_input_file("modinputbits_v0.1.bin", 300)