file based GMSK routines
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Sun, 25 Jan 2015 08:14:48 +0000 (08:14 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Sun, 25 Jan 2015 08:14:48 +0000 (08:14 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@2027 01035d8c-6547-0410-b346-abe4f91aad63

codec2-dev/octave/gmsk.m

index cb040fa1ef10d1c312212d6dc0e306a3134057ff..5a6be8d4087d8fcbcfb50e72105f4afe5101d5bf 100644 (file)
@@ -22,7 +22,7 @@
 %     + Maybe based on test tone/carrier from the other side?
 %     + think about process ... total signal plus noise power?  Increase power until S+N doubles?
 % [ ] generate curves for baseline modem and with sync algorithms
-%     + used coarse sync code to remove need for knowing delays
+%     [ ] used coarse sync code to remove need for knowing delays
 
 % Filter coeffs From:
 % https://github.com/on1arf/gmsk/blob/master/gmskmodem_codec2/API/a_dspstuff.h,
@@ -85,6 +85,7 @@ function gmsk_states = gmsk_init(gmsk_states, Rs)
 
   % general 
 
+  verbose = gmsk_states.verbose;
   gmsk_states.Fs = 48000;
   gmsk_states.Rs = Rs;
   M = gmsk_states.M = gmsk_states.Fs/gmsk_states.Rs;
@@ -92,12 +93,14 @@ function gmsk_states = gmsk_init(gmsk_states, Rs)
   global gmsk_demod_coeff;
   gmsk_states.mod_coeff = (Rs/4800)*resample(gmsk_mod_coeff, 4800, Rs);
 
-  figure(12)
-  plot(gmsk_mod_coeff,'r;original 4800;')
-  hold on;
-  plot(gmsk_states.mod_coeff,'g;interpolated;')
-  hold off;
-  title('GMSK pulse shaping filter')
+  if verbose
+    figure(12)
+    plot(gmsk_mod_coeff,'r;original 4800;')
+    hold on;
+    plot(gmsk_states.mod_coeff,'g;interpolated;')
+    hold off;
+    title('GMSK pulse shaping filter')
+  end
 
   % set up FM modulator
 
@@ -117,6 +120,7 @@ function [tx tx_filt tx_symbols] = gmsk_mod(gmsk_states, tx_bits)
   M = gmsk_states.M;
   nsym = length(tx_bits);
   nsam = nsym*M;
+  verbose = gmsk_states.verbose;
 
   % NRZ sequence of symbols
 
@@ -127,10 +131,12 @@ function [tx tx_filt tx_symbols] = gmsk_mod(gmsk_states, tx_bits)
 
   tx_filt = filter(gmsk_states.mod_coeff, 1, tx_symbols);
   
-  figure(13)
-  clf
-  plot(tx_filt(1:M*10))
-  title('tx signal after filtering, before FM mod')
+  if verbose
+    figure(13)
+    clf
+    plot(tx_filt(1:M*10))
+    title('tx signal after filtering, before FM mod')
+  end
 
   tx = analog_fm_mod(gmsk_states.fm_states, tx_filt);
 endfunction
@@ -506,6 +512,8 @@ function run_gmsk_single
 endfunction
 
 
+% Generate a bunch of BER versus Eb/No curves for various demods
+
 function run_gmsk_curves
   sim_in.coherent_demod = 1;
   sim_in.nsym = 2400;
@@ -557,13 +565,73 @@ function run_gmsk_curves
 
 endfunction
 
-% TODO: [ ] test over range of freq, phase, coarse timing and Eb/No values
-%       [ ] Modify for Rs=1200, e.g. mod filter above
-%       [ ] extra stuff at begining/end for filter delays
+
+function [preamble_location freq_offset_est] = find_preamble(gmsk_states, M, npreamble, rx)
+  verbose = gmsk_states.verbose;
+
+  % look through rx buffer and determine if there is a valid preamble.  Use steps of half the
+  % preamble size in samples to try to bracket the pre-amble.
+
+  preamble_step = npreamble*M/2;
+  ratio = 0; freq_offset_est = 0; preamble_location = 0;
+  ratio_log = [];
+  for i=1:preamble_step:length(rx)-preamble_step
+    [afreq_offset_est aratio] = gmsk_est_freq_offset(gmsk_states, rx(i:i+preamble_step-1), verbose);
+    ratio_log = [ratio_log aratio];
+    if aratio > ratio
+      preamble_location = i;
+      ratio = aratio;
+      freq_offset_est = afreq_offset_est;
+    end
+  end
+  if verbose
+    printf("preamble location: %d frequency offset: %f ratio: %f\n", 
+    preamble_location, freq_offset_est, ratio);   
+    figure(9)
+    plot(ratio_log);
+    title('Preamble ratio');
+  end
+endfunction
+
+
+% attempt to perform "coarse sync" sync with the received frames, we
+% check each frame for the best coarse sync position.  Brute force
+% approach, that would be changed for a real demod which has some
+% sort of unique word.  Start looking for valid frames 1 frame
+% after start of pre-amble to give PLL time to lock
+
+function [total_errors total_bits Nerrs_log Nerrs_all_log] = coarse_sync_ber(nframes_rx, tx_frame, rx_bits)
+
+  Nerrs_log = zeros(1, nframes_rx);
+  Nerrs_all_log = zeros(1, nframes_rx);
+  total_errors = 0;
+  total_bits   = 0;
+  framesize = length(tx_frame);
+
+  for f=2:nframes_rx-1
+    Nerrs_min = framesize;
+    for i=1:framesize;
+      st = (f-1)*framesize+i; en = st+framesize-1;
+      Nerrs = sum(xor(rx_bits(st:en), tx_frame));
+      if Nerrs < Nerrs_min
+        Nerrs_min = Nerrs;
+      end
+    end
+    Nerrs_all_log(f) = Nerrs_min;
+    if Nerrs_min/framesize < 0.1
+      Nerrs_log(f) = Nerrs_min;
+      total_errors += Nerrs_min;
+      total_bits   += framesize;
+    end
+  end
+endfunction
+
+
+% Give the demod a hard time: frequency, phase, time offsets, sample clock difference
    
-function run_test_freq_offset
-  Rs = 4800;
-  verbose = 1;
+function run_test_channel_impairments
+  Rs = 1200;
+  verbose = 2;
   aEbNodB = 6;
   phase_offset = pi/2;
   freq_offset  = -104;
@@ -572,6 +640,7 @@ function run_test_freq_offset
   nsym = 4800*2;
   npreamble = 480;
 
+  gmsk_states.verbose = verbose;
   gmsk_states.coherent_demod = 1;
   gmsk_states.phase_track    = 1;
   gmsk_states = gmsk_init(gmsk_states, Rs);
@@ -580,10 +649,10 @@ function run_test_freq_offset
   M  = gmsk_states.M;
 
   % A frame consists of nsym random data bits.  Some experimentation
-  % has shown they need to be random for the sync algorithms to work
-
-  % note must be random-ish data (not say 11001100...) for timing estimator to work.
-  % However initial freq offset estimation is a lot easier with a 01010 type sequence
+  % has shown they must be random-ish data (not say 11001100...) for
+  % timing estimator to work.  However initial freq offset estimation
+  % is a lot easier with a 01010 type sequence, so we construct a 
+  % frame with a pre-amble followed by frames of random data.
 
   framesize = 480;
   nframes = floor(nsym/framesize);
@@ -600,6 +669,11 @@ function run_test_freq_offset
   tx = [zeros(1,timing_offset) tx];
   nsam = length(tx);
 
+  fc = 1500; gain = 2;
+  wc = 2*pi*fc/Fs;
+  w1 = exp(j*wc*(1:nsam));
+  tx = gain*real(tx .* w1);
+
   if verbose
     figure(11);
     subplot(211)
@@ -618,29 +692,14 @@ function run_test_freq_offset
 
   rx = tx.*exp(j*w) + noise;
 
-  % look through rx buffer and determine if there is a valid preamble.  Use steps of half the
-  % preamble size in samples to try to bracket the pre-amble.
+  rx = rx .* conj(w1);
 
-  preamble_step = npreamble*M/2;
-  ratio = 0; freq_offset_est = 0; preamble_location = 0;
-  ratio_log = [];
-  for i=1:preamble_step:length(rx)-preamble_step
-    [afreq_offset_est aratio] = gmsk_est_freq_offset(gmsk_states, rx(i:i+preamble_step-1), verbose);
-    ratio_log = [ratio_log aratio];
-    if aratio > ratio
-      preamble_location = i;
-      ratio = aratio;
-      freq_offset_est = afreq_offset_est;
-    end
-  end
-  if verbose
-    printf("preamble location: %d frequency offset: %f ratio: %f\n", 
-    preamble_location, freq_offset_est, ratio);   
-    figure(9)
-    plot(ratio_log);
-    title('Preamble ratio');
-  end
+  figure(3)
+  clf
+  Rx=20*log10(abs(fft(rx)));
+  plot(Rx);
 
+  [preamble_location freq_offset_est] = find_preamble(gmsk_states, M, npreamble, rx);
   w_est  = (0:nsam-1)*2*pi*freq_offset_est/Fs;
   rx = rx.*exp(-j*w_est);
 
@@ -651,34 +710,7 @@ function run_test_freq_offset
 
   % printf("ntx: %d nrx: %d ntx_bits: %d nrx_bits: %d\n", length(tx), length(rx), length(tx_bits), length(rx_bits));
 
-  % attempt to perform "coarse sync" sync with the received frames, we
-  % check each frame for the best coarse sync position.  Brute force
-  % approach, that would be changed for a real demod which has some
-  % sort of unique word.  Start looking for valid frames 1 frame
-  % after start of pre-amble to give PLL time to lock
-
-  Nerrs_log = zeros(1,nframes_rx);
-  Nerrs_all_log = zeros(1,nframes_rx);
-  total_errors = 0;
-  total_bits   = 0;
-  
-  for f=2:nframes_rx-1
-    Nerrs_min = framesize;
-    for i=1:framesize;
-      st = (f-1)*framesize+i; en = st+framesize-1;
-      Nerrs = sum(xor(rx_bits(st:en), tx_frame));
-      %printf("nframes: %d f: %f st: %d en: %d Nerrs: %d\n", nframes, f, st, en, Nerrs);
-      if Nerrs < Nerrs_min
-        Nerrs_min = Nerrs;
-      end
-    end
-    Nerrs_all_log(f) = Nerrs_min;
-    if Nerrs_min/framesize < 0.1
-      Nerrs_log(f) = Nerrs_min;
-      total_errors += Nerrs_min;
-      total_bits   += framesize;
-    end
-  end
+  [total_errors total_bits Nerrs_log Nerrs_all_log] = coarse_sync_ber(nframes_rx, tx_frame, rx_bits);
 
   ber = total_errors/total_bits;
 
@@ -699,8 +731,122 @@ function run_test_freq_offset
 
 endfunction
     
+
+% Generates a Fs=48kHz raw file of 16 bit samples centred on 1500Hz,
+% Suitable for transmitting with a SSB tx
+
+function gmsk_tx(tx_file_name)
+  Rs = 1200;
+  nsym =  Rs*10;
+  framesize = 480;
+  npreamble = 480;
+  gain      = 10000;
+  fc        = 1500;
+
+  gmsk_states.verbose = 0;
+  gmsk_states.coherent_demod = 1;
+  gmsk_states.phase_track    = 1;
+  gmsk_states = gmsk_init(gmsk_states, Rs);
+  Fs = gmsk_states.Fs;
+  Rs = gmsk_states.Rs;
+  M  = gmsk_states.M;
+
+  % generate frame with preamble
+
+  nframes = floor(nsym/framesize);
+  tx_frame = round(rand(1, framesize));
+  tx_bits = zeros(1,npreamble);
+  tx_bits(1:2:npreamble) = 1;
+  for i=1:nframes
+    tx_bits = [tx_bits tx_frame];
+  end
+  
+  [tx tx_filt tx_symbols] = gmsk_mod(gmsk_states, tx_bits);
+  nsam = length(tx);
+
+  wc = 2*pi*fc/Fs;
+  w = exp(j*wc*(1:nsam));
+  tx = gain*real(tx .* w);
+  figure(1)
+  plot(tx(1:4000))
+  fout = fopen(tx_file_name,"wb");
+  fwrite(fout, tx, "short");
+  fclose(fout);
+
+endfunction
+
+
+% Reads a file of Fs=48kHz 16 bit samples centred on 1500Hz, and
+% measures the BER.
+
+function gmsk_rx(rx_file_name)
+  Rs = 1200;
+  framesize = 480;
+  npreamble = 480;
+  fc        = 1500;
+
+  gmsk_states.verbose = 2;
+  gmsk_states.coherent_demod = 1;
+  gmsk_states.phase_track    = 1;
+  gmsk_states = gmsk_init(gmsk_states, Rs);
+  Fs = gmsk_states.Fs;
+  Rs = gmsk_states.Rs;
+  M  = gmsk_states.M;
+
+  tx_frame = round(rand(1, framesize));
+
+  % get real signal at fc offset and convert to baseband complex
+  % signal
+  
+  fin = fopen(rx_file_name,"rb");
+  rx = fread(fin,"short")';
+  fclose(fin);
+  nsam = length(rx);
+  wc = 2*pi*fc/Fs;
+  w = exp(-j*wc*(1:nsam));
+  rx = rx .* w;
+
+  figure(3)
+  clf
+  Rx=20*log10(abs(fft(rx)));
+  plot(Rx);
+
+  % find preamble and demodulate
+
+  [preamble_location freq_offset_est] = find_preamble(gmsk_states, M, npreamble, rx);
+  w_est  = (0:nsam-1)*2*pi*freq_offset_est/Fs;
+  rx = rx.*exp(-j*w_est);
+  [rx_bits rx_out rx_filt] = gmsk_demod(gmsk_states, rx(preamble_location+framesize:nsam));
+  nframes_rx = length(rx_bits)/framesize;
+
+  % count errors
+
+  [total_errors total_bits Nerrs_log Nerrs_all_log] = coarse_sync_ber(nframes_rx, tx_frame, rx_bits);
+
+  ber = total_errors/total_bits;
+
+  printf("f_off: %4.1f Nframes: %d Nbits: %d Nerrs: %d BER: %f\n", 
+         freq_offset_est, nframes_rx, total_bits, total_errors, ber);
+
+  figure(2)
+  clf
+  subplot(211)
+  plot(Nerrs_log);
+  hold on;
+  plot(Nerrs_all_log,'g');
+  hold off;
+  title('Bit Errors')
+  subplot(212)
+  stem(real(cumsum(Nerrs_log)))
+  title('Cumulative Bit Errors')
+
+endfunction
+
+
 %run_gmsk_single
 %run_gmsk_curves
 %run_gmsk_init
-run_test_freq_offset
+%run_test_channel_impairments
+%gmsk_tx("test_gmsk.raw")
+gmsk_rx("test_gmsk.raw")