From: baobrien <baobrien@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Fri, 6 Nov 2015 23:22:47 +0000 (+0000)
Subject: fine timing work done on ya_fsk
X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=05bc86b553814d9e76324bce1692ca1415755189;p=freetel-svn-tracking.git

fine timing work done on ya_fsk

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

diff --git a/codec2-dev/octave/fsk4.m b/codec2-dev/octave/fsk4.m
index 9a4f3491..8891a455 100644
--- a/codec2-dev/octave/fsk4.m
+++ b/codec2-dev/octave/fsk4.m
@@ -174,6 +174,30 @@ function syms = mrd4(bits)
   end
 endfunction
 
+% Minimal Running Disparity, 8 symbol encoder
+% This is a simple 2 bit to 1 symbol encoding for 8fsk modems built
+% on old fashoned FM radios.
+function syms = mrd8(bits)
+  bitlen = length(bits);
+  if mod(bitlen,2) == 1
+    bits = [bits 0]
+  endif
+
+  syms = zeros(1,length(bits)*.5);
+  rd=0;
+  lastsym=0;
+  for n = (1:2:length(bits))
+    bit = (bits(n)*2)+bits(n+1);
+    sp = [1 3 7 5](bit+1); %Map a bit to a +1 or +3
+    [x,v] = min(abs([rd+sp rd-sp])); %Select +n or -n, whichever minimizes disparity
+    ssel = [sp -sp](v);
+    if(ssel == lastsym)ssel = -ssel;endif %never run 2 of the same syms in a row
+    syms((n+1)/2) = ssel; %emit the symbol
+    rd = rd + ssel; %update running disparity
+    lastsym = ssel; %remember this symbol for next time
+  end
+endfunction
+
 % "Manchester 4" encoding
 function syms = mane4(bits)
     syms = zeros(1,floor(bits/2)*2);
diff --git a/codec2-dev/octave/ya_fsk.m b/codec2-dev/octave/ya_fsk.m
deleted file mode 100644
index f22923b2..00000000
--- a/codec2-dev/octave/ya_fsk.m
+++ /dev/null
@@ -1,161 +0,0 @@
-% yafsk.m
-% Yet-another-FSK
-% Brady O'Brien 20 October 2015
-%
-% This is a model for the first attempt at a C FSK modem. Based on fsk_horus and maybe a little fsk4.
-% First revision will just be 2fsk
-% [x] - Modem init and struct def
-% [x] - Direct SDR modulator, probably not FM based
-% [.] - Direct SDR non-coherent demodulator
-%    [ ] - Core demodulation routine
-%    [ ] - Timing offset estimation
-%    [ ] - Freq. offset estimation
-%    [ ] - Bit slip, maybe
-% [ ] - The C port
-% [ ] - Some stuff to verify the C port
-
-%clear all;
-graphics_toolkit('gnuplot');
-fm
-
-%Basic parameters for a simple FSK modem
-fsk_setup_info.Rs = 2400;  % Symbol rate
-fsk_setup_info.nfsk = 2;      % Number of unique symbols. Must be 2.
-fsk_setup_info.Fs = 48000; % Sample frequency
-fsk_setup_info.Fsym = fsk_setup_info.Rs; %Symbol spacing
-fsk_setup_info.txmap = @(bits) bits+1; %Map TX bits to 2fsk symbols
-fsk_setup_info.rxmap = @(syms) syms==2; %Map 2fsk RX symbols to bits
-
-
-function states = yafsk_init(fsk_config)
-  Fs = states.Fs = fsk_config.Fs;
-  Rs = states.Rs = fsk_config.Rs;
-  nfsk = states.nfsk = fsk_config.nfsk;
-  Ts = states.Ts = Fs/Rs;
-  Fsym = states.Fsym = fsk_config.Fsym;
-  states.config = fsk_config;
-
-  if nfsk != 2
-    error("Gotta be 2fsk")
-  endif
- 
-  %Symbol frequencies. Fixed to intervals of Fsym.
-  states.fsyms = [-(Fsym/2) (Fsym/2)];
-  states.tx_phase = 0;
-
-  states.dc = zeros(1,nfsk);
-  states.rx_phi = ones(1,nfsk);
-  states.isamp = 0;
-  states.sums = zeros(1,nfsk);
-endfunction
-
-function [tx states] = yafsk_mod(states,bits)
-  Ts = states.Ts;
-  Fs = states.Fs;
-  fsyms = states.fsyms;
-  tx_phase = states.tx_phase;
-  %Map bits into symbols
-  syms = states.config.txmap(bits);
-  tx = zeros(1,Ts*length(syms));
-  
-  for ii = (1:length(syms))
-    cur_sym_f = fsyms(syms(ii));
-    tx_phase_i = tx_phase;
-    for jj = (1:Ts)
-        tx_phase_i = tx_phase + jj*2*pi*cur_sym_f/Fs;
-        tx((ii-1)*Ts+jj) = exp(j*tx_phase_i);   
-    end
-    tx_phase = tx_phase + Ts*2*pi*cur_sym_f/Fs;
-    if tx_phase>2*pi
-        tx_phase = tx_phase-2*pi;
-    elseif tx_phase<-2*pi
-        tx_phase = tx_phase+2*pi;
-    endif
-    %tx_phase_vec = tx_phase + (1:Ts)*2*pi*cur_sym_f/Fs;
-    %tx((ii-1)*Ts+1:ii*Ts) = exp(j*tx_phase_vec);
-  end
-  states.tx_phase = tx_phase;
-endfunction
-
-function d = idmp(data, M)
-    d = zeros(1,length(data)/M);
-    for i = 1:length(d)
-      d(i) = sum(data(1+(i-1)*M:i*M));
-    end
-endfunction
-
-function [bits states] = yafsk_demod_2a(states,rx)
-  fine_timing = 1;
-  Fs = states.Fs;
-  Rs = states.Rs;
-  Ts = states.Ts;
-  nfsk = states.nfsk;
-
-  phy_f1 = states.rx_phi(1);
-  phy_f2 = states.rx_phi(2);
-
-  dphase_f1 = exp(states.fsyms(1)*-j*2*pi/Fs);
-  dphase_f2 = exp(states.fsyms(2)*-j*2*pi/Fs);
-
-  sum_f1 = states.sums(1);
-  sum_f2 = states.sums(2);
-
-  isamp = states.isamp;
-  symcnt = 1;
-  syms = [0];
-  sums1 = [0];
-  sums2 = [0];
-  isamp = 1;
-  for ii = (1:length(rx))
-    phy_f1 *= dphase_f1;   %Spin the oscillators
-    phy_f2 *= dphase_f2;
-
-    dcs_f1 = rx(ii)*phy_f1; %Figure out the DC
-    dcs_f2 = rx(ii)*phy_f2;
-
-    sum_f1 += dcs_f1; %Integrate
-    sum_f2 += dcs_f2;
-
-    isamp += 1;
-    if isamp==Ts %If it's time to take a sample and spit out a symbol..
-        syms(symcnt) = (abs(sum_f1)>abs(sum_f2))+1; %Spit out a symbol
-        symcnt += 1;
-
-        %sums1(symcnt) = abs(sum_f1);
-	%sums2(symcnt) = abs(sum_f2);
-
-        sum_f1 = 0;   %Reset integrators
-        sum_f2 = 0;
-        isamp = 0;    %Reset integrator count
-    endif
-  end
-  plot((1:length(sums1)),sums1,(1:length(sums2)),sums2);
-
-  bits = states.config.rxmap(syms);
-endfunction
-function [bits states] = yafsk_demod_2(states,rx)
-  fine_timing = 1;
-  Fs = states.Fs;
-  Rs = states.Rs;
-  Ts = states.Ts;
-  nfsk = states.nfsk;
-
-  rx = rx(fine_timing:length(rx));
-  sym_phases = (1:length(rx)).*rot90(states.fsyms)*2*pi/Fs;
-
-  sym_mixers = exp(-j*sym_phases);
-  rx_mixd = repmat(rx,nfsk,1).*sym_mixers;
-  
-  dc1 = abs(idmp(rx_mixd(1,1:length(rx_mixd)),Ts));
-  dc2 = abs(idmp(rx_mixd(2,1:length(rx_mixd)),Ts));
-
-  t=(1:length(dc1));
-  
-  plot(t,dc1,t,dc2)
-
-endfunction
-
-
-
-
-
diff --git a/codec2-dev/octave/yafsk.m b/codec2-dev/octave/yafsk.m
new file mode 100644
index 00000000..fe7745c3
--- /dev/null
+++ b/codec2-dev/octave/yafsk.m
@@ -0,0 +1,554 @@
+% yafsk.m
+% Yet-another-FSK
+% Brady O'Brien 20 October 2015
+%
+% This is a model for the first attempt at a C FSK modem. Based on fsk_horus and maybe a little fsk4.
+% First revision will just be 2fsk
+% [x] - Modem init and struct def
+% [x] - Direct SDR modulator, probably not FM based
+% [x] - Direct SDR non-coherent demodulator
+%    [x] - Core demodulation routine
+%    [ ] - Timing offset estimation
+%    [ ] - Freq. offset estimation
+%    [ ] - Bit slip, maybe
+% [ ] - The C port
+% [ ] - Some stuff to verify the C port
+
+%clear all;
+graphics_toolkit('gnuplot');
+%fm
+
+%Basic parameters for a simple FSK modem
+fsk_setup_info.Rs = 2400;  % Symbol rate
+fsk_setup_info.nfsk = 2;      % Number of unique symbols. Must be 2.
+fsk_setup_info.Fs = 48000; % Sample frequency
+fsk_setup_info.Fsym = fsk_setup_info.Rs; %Symbol spacing
+fsk_setup_info.txmap = @(bits) bits+1; %Map TX bits to 2fsk symbols
+fsk_setup_info.rxmap = @(syms) syms==1; %Map 2fsk RX symbols to bits
+global fsk_setup = fsk_setup_info
+
+function states = yafsk_init(fsk_config)
+  Fs = states.Fs = fsk_config.Fs;
+  Rs = states.Rs = fsk_config.Rs;
+  nfsk = states.nfsk = fsk_config.nfsk;
+  Ts = states.Ts = Fs/Rs;
+  Fsym = states.Fsym = fsk_config.Fsym;
+  states.config = fsk_config;
+
+  if nfsk != 2
+    error("Gotta be 2fsk")
+  endif
+ 
+  %Symbol frequencies. Fixed to intervals of Fsym.
+  states.fsyms = [-(Fsym/2) (Fsym/2)];
+  states.tx_phase = 0;
+
+  states.dc = zeros(1,nfsk);
+  states.rx_phi = ones(1,nfsk);
+  states.isamp = 0;
+  states.sums = zeros(1,nfsk);
+endfunction
+
+function [tx states] = yafsk_mod(states,bits)
+  Ts = states.Ts;
+  Fs = states.Fs;
+  fsyms = states.fsyms;
+  tx_phase = states.tx_phase;
+  %Map bits into symbols
+  syms = states.config.txmap(bits);
+  tx = zeros(1,Ts*length(syms));
+  
+  for ii = (1:length(syms))
+    cur_sym_f = fsyms(syms(ii));
+    tx_phase_i = tx_phase;
+    for jj = (1:Ts)
+        tx_phase_i = tx_phase + jj*2*pi*cur_sym_f/Fs;
+        tx((ii-1)*Ts+jj) = exp(j*tx_phase_i);   
+    end
+    tx_phase = tx_phase + Ts*2*pi*cur_sym_f/Fs;
+    if tx_phase>2*pi
+        tx_phase = tx_phase-2*pi;
+    elseif tx_phase<-2*pi
+        tx_phase = tx_phase+2*pi;
+    endif
+    %tx_phase_vec = tx_phase + (1:Ts)*2*pi*cur_sym_f/Fs;
+    %tx((ii-1)*Ts+1:ii*Ts) = exp(j*tx_phase_vec);
+  end
+  states.tx_phase = tx_phase;
+endfunction
+
+function d = idmp(data, M)
+    d = zeros(1,length(data)/M);
+    for i = 1:length(d)
+      d(i) = sum(data(1+(i-1)*M:i*M));
+    end
+endfunction
+
+function [bits states phi phis] = yafsk_demod_2a(states,rx)
+  fine_timing = 1;
+  Fs = states.Fs;
+  Rs = states.Rs;
+  Ts = states.Ts;
+  nfsk = states.nfsk;
+
+
+
+  phy_f1 = states.rx_phi(1);
+  phy_f2 = states.rx_phi(2);
+
+  dphase_f1 = exp(states.fsyms(1)*-j*2*pi/Fs);
+  dphase_f2 = exp(states.fsyms(2)*-j*2*pi/Fs);
+
+  sum_f1 = states.sums(1);
+  sum_f2 = states.sums(2);
+
+  isamp = states.isamp;
+  symcnt = 1;
+  subcnt = 1;
+  syms = [0];
+  sums1 = [0];
+  sums2 = [0];
+  phis = [0];
+  isamp = 1;
+  isub = 1;
+
+  ssum_f1 = 0;
+  ssum_f2 = 0;
+  ssamp=0;
+ 
+  timing_syms = 10;
+  timing_nudge = .085; %How far to 'nudge' the sampling point
+  timing_samps = timing_syms*(Ts/4);
+  timing_db = zeros(1,timing_samps);  
+
+  for ii = (1:length(rx))
+    phy_f1 *= dphase_f1;   %Spin the oscillators
+    phy_f2 *= dphase_f2;
+
+    dcs_f1 = rx(ii)*phy_f1; %Figure out the DC
+    dcs_f2 = rx(ii)*phy_f2;
+
+    sum_f1 += dcs_f1; %Integrate
+    sum_f2 += dcs_f2;
+    ssum_f1 += dcs_f1;
+    ssum_f2 += dcs_f2;
+
+    ssamp += 1;
+    isamp += 1;
+    if isamp>=Ts %If it's time to take a sample and spit out a symbol..
+      syms(symcnt) = (abs(sum_f1)>abs(sum_f2))+1; %Spit out a symbol
+      symcnt += 1;
+      tdmd = timing_db.*exp(-j*2*pi*(1:length(timing_db))*24000/Fs);
+
+      timing_phi = angle(sum(tdmd));
+      phis(symcnt) = timing_phi;
+      sums1(symcnt) = sum_f1;
+      sums2(symcnt) = sum_f2;
+      sum_f1 = 0;   %Reset integrators
+      sum_f2 = 0;
+      isamp -= Ts;    %Reset integrator count
+      phi_tgt = 1.252;
+      if(timing_phi>-1.88 && timing_phi<1.252)
+          isamp += timing_nudge;
+          ssamp += timing_nudge;
+      else
+          isamp -= timing_nudge;
+          ssamp -= timing_nudge;
+      endif
+    endif
+
+    if ssamp>= Ts/4
+      %sums1(subcnt) = ssum_f1;
+      %sums2(subcnt) = ssum_f2;
+      timing_db(1:timing_samps-1)=timing_db(2:timing_samps);
+      timing_db(timing_samps) = (abs(ssum_f1)-ssum_f2).^2;
+
+      subcnt += 1;
+      ssum_f1 = 0;
+      ssum_f2 = 0;
+      ssamp -= Ts/4;
+    endif
+
+  end
+
+  nlin = (abs(sums1)-sums2).^2;
+  %nlin = nlin.*exp(-j*2*pi*(1:length(nlin))*000/Fs);
+  nlen = length(nlin)
+  fto = sum(nlin)
+  fta = angle(fto)
+  phi = fto;
+
+  figure(1);
+  nlf = fft(nlin);
+  length(nlf)
+  plot(20*log10(abs(nlf)));
+
+  states.rx_phy(1) = phy_f1;
+  states.rx_phy(2) = phy_f2;
+ 
+  states.sums(1) = sum_f1;
+  states.sums(2) = sum_f2;
+
+  states.isamp = isamp;
+
+  figure(2)
+  plot((1:length(sums1)),sums1,(1:length(sums2)),sums2);
+
+  bits = states.config.rxmap(syms);
+endfunction
+function [bits states] = yafsk_demod_2(states,rx)
+  fine_timing = 1;
+  Fs = states.Fs;
+  Rs = states.Rs;
+  Ts = states.Ts;
+  nfsk = states.nfsk;
+
+  rx = rx(fine_timing:length(rx));
+  sym_phases = (1:length(rx)).*rot90(states.fsyms)*2*pi/Fs;
+
+  sym_mixers = exp(-j*sym_phases);
+  rx_mixd = repmat(rx,nfsk,1).*sym_mixers;
+  
+  dc1 = abs(idmp(rx_mixd(1,1:length(rx_mixd)),Ts));
+  dc2 = abs(idmp(rx_mixd(2,1:length(rx_mixd)),Ts));
+
+  t=(1:length(dc1));
+  
+  plot(t,dc1,t,dc2)
+
+endfunction
+
+% Simulation thing, shamelessly taken from fsk_horus.m
+% simulation of tx and rx side, add noise, channel impairments ----------------------
+
+function run_sim
+  global fsk_setup;
+  frames = 10;
+  EbNodB = 26;
+  timing_offset = 0.0; % see resample() for clock offset below
+  test_frame_mode = 2;
+  fading = 0;          % modulates tx power at 2Hz with 20dB fade depth, 
+                       % to simulate balloon rotating at end of mission
+  df     = 0;          % tx tone freq drift in Hz/s
+
+  more off
+  rand('state',1); 
+  randn('state',1);
+  states = yafsk_init(fsk_setup);
+  N = states.N;
+  %P = states.P;
+  Rs = states.Rs;
+  nsym = states.nsym = 4096;
+  Fs = states.Fs;
+  states.df = df;
+
+  EbNo = 10^(EbNodB/10);
+  variance = states.Fs/(states.Rs*EbNo);
+
+  % set up tx signal with payload bits based on test mode
+
+  if test_frame_mode == 1
+     % test frame of bits, which we repeat for convenience when BER testing
+    test_frame = round(rand(1, states.nsym));
+    tx_bits = [];
+    for i=1:frames+1
+      tx_bits = [tx_bits test_frame];
+    end
+  end
+  if test_frame_mode == 2
+    % random bits, just to make sure sync algs work on random data
+    tx_bits = round(rand(1, states.nsym*(frames+1)));
+  end
+  if test_frame_mode == 3
+    % ...10101... sequence
+    tx_bits = zeros(1, states.nsym*(frames+1));
+    tx_bits(1:2:length(tx_bits)) = 1;
+  end
+
+  if test_frame_mode == 4
+
+    % load up a horus msg from disk and modulate that
+
+    test_frame = load("horus_msg.txt");
+    ltf = length(test_frame);
+    ntest_frames = ceil((frames+1)*nsym/ltf);
+    tx_bits = [];
+    for i=1:ntest_frames
+      tx_bits = [tx_bits test_frame];
+    end
+  end
+
+  tx = yafsk_mod(states, tx_bits);
+
+  %tx = resample(tx, 1000, 1001); % simulated 1000ppm sample clock offset
+
+  if fading
+     ltx = length(tx);
+     tx = tx .* (1.1 + cos(2*pi*2*(0:ltx-1)/Fs))'; % min amplitude 0.1, -20dB fade, max 3dB
+  end
+
+  noise = sqrt(variance)*randn(length(tx),1);
+  rx    = tx + noise;
+  
+  % dump simulated rx file
+  ftx=fopen("ya_fsk_rx.raw","wb"); rxg = rx*1000; fwrite(ftx, rxg, "short"); fclose(ftx);
+
+  timing_offset_samples = round(timing_offset*states.Ts);
+  st = 1 + timing_offset_samples;
+  rx_bits_buf = zeros(1,2*nsym);
+  Terrs = Tbits = 0;
+  state = 0;
+  x_log = [];
+  norm_rx_timing_log = [];
+  nerr_log = [];
+  f1_int_resample_log = [];
+  f2_int_resample_log = [];
+  f1_log = f2_log = [];
+  EbNodB_log = [];
+  rx_bits_log = [];
+  rx_bits_sd_log = [];
+
+  for f=1:frames
+
+    % extract nin samples from input stream
+
+    nin = states.nin;
+    en = st + states.nin - 1;
+    sf = rx(st:en);
+    st += nin;
+
+    % demodulate to stream of bits
+
+    [rx_bits states] = yafsk_demod(states, sf);
+    rx_bits_buf(1:nsym) = rx_bits_buf(nsym+1:2*nsym);
+    rx_bits_buf(nsym+1:2*nsym) = rx_bits;
+    rx_bits_log = [rx_bits_log rx_bits];
+    rx_bits_sd_log = [rx_bits_sd_log states.rx_bits_sd];
+
+    norm_rx_timing_log = [norm_rx_timing_log states.norm_rx_timing];
+    x_log = [x_log states.x];
+    f1_int_resample_log = [f1_int_resample_log abs(states.f1_int_resample)];
+    f2_int_resample_log = [f2_int_resample_log abs(states.f2_int_resample)];
+    f1_log = [f1_log states.f1];
+    f2_log = [f2_log states.f2];
+    EbNodB_log = [EbNodB_log states.EbNodB];
+
+    % frame sync based on min BER
+
+    if test_frame_mode == 1
+      nerrs_min = nsym;
+      next_state = state;
+      if state == 0
+        for i=1:nsym
+          error_positions = xor(rx_bits_buf(i:nsym+i-1), test_frame);
+          nerrs = sum(error_positions);
+          if nerrs < nerrs_min
+            nerrs_min = nerrs;
+            coarse_offset = i;
+          end
+        end
+        if nerrs_min < 3
+          next_state = 1;
+          %printf("%d %d\n", coarse_offset, nerrs_min);
+        end
+      end
+
+      if state == 1  
+        error_positions = xor(rx_bits_buf(coarse_offset:coarse_offset+nsym-1), test_frame);
+        nerrs = sum(error_positions);
+        Terrs += nerrs;
+        Tbits += nsym;
+        nerr_log = [nerr_log nerrs];
+      end
+
+      state = next_state;
+    end
+  end
+
+  if test_frame_mode == 1
+    printf("frames: %d Tbits: %d Terrs: %d BER %4.3f\n", frames, Tbits, Terrs, Terrs/Tbits);
+  end
+
+  if test_frame_mode == 4
+    extract_and_print_packets(states, rx_bits_log, rx_bits_sd_log)
+  end
+
+  figure(1);
+  plot(f1_int_resample_log,'+')
+  hold on;
+  plot(f2_int_resample_log,'g+')
+  hold off;
+
+  figure(2)
+  clf
+  m = max(abs(x_log));
+  plot(x_log,'+')
+  axis([-m m -m m])
+  title('fine timing metric')
+
+  figure(3)
+  clf
+  subplot(211)
+  plot(norm_rx_timing_log);
+  axis([1 frames -1 1])
+  title('norm fine timing')
+  subplot(212)
+  plot(nerr_log)
+  title('num bit errors each frame')
+
+  figure(4)
+  clf
+  plot(real(rx(1:Fs)))
+  title('rx signal at demod input')
+
+  figure(5)
+  clf
+  plot(f1_log)
+  hold on;
+  plot(f2_log,'g');
+  hold off;
+  title('tone frequencies')
+  axis([1 frames 0 Fs/2])
+
+  figure(6)
+  clf
+  plot(EbNodB_log);
+  title('Eb/No estimate')
+
+endfunction
+
+
+% Bit error rate test ----------------------------------------------------------
+% Params - aEsNodB - EbNo in decibels
+%        - timing_offset - how far the fine timing is offset
+%        - bitcnt - how many bits to check
+%        - demod_fx - demodulator function
+% Returns - ber - teh measured BER
+%         - thrcoh - theory BER of a coherent demod
+%         - thrncoh - theory BER of non-coherent demod
+function [ber thrcoh thrncoh rxphis] = nfbert_2(aEsNodB,modem_config, bitcnt=100000, timing_offset = 1)
+
+  rand('state',1); 
+  randn('state',1);
+  
+  %How many bits should this test run?
+  bitcnt = 12000;
+  
+  test_bits = [repmat([1 0],1,10) rand(1,bitcnt)>.5]; %Random bits. Pad with zeros to prime the filters
+  states.M = 1;
+  states = yafsk_init(modem_config);
+  
+  %Set this to 0 to cut down on the plotting
+  states.verbose = 1; 
+  Fs = states.Fs;
+  Rb = states.Rs;  % Multiply symbol rate by 2, since we have 2 bits per symbol
+  
+  tx = yafsk_mod(states,test_bits);
+
+  %add noise here
+  %shamelessly copied from gmsk.m
+  EsNo = 10^(aEsNodB/10);
+  EbNo = EsNo
+  variance = Fs/(Rb*EbNo);
+  nsam = length(tx);
+  noise = sqrt(variance/2)*(randn(1,nsam) + j*randn(1,nsam));
+  rx    = tx*exp(j*pi/2) + noise;
+
+  rx    = rx(timing_offset:length(rx));
+
+  [rx_bits states a rxphis] = yafsk_demod_2a(states,rx);
+  ber = 1;
+  
+  %thing to account for offset from input data to output data
+  %No preamble detection yet
+  plot(rxphis);
+  ox = 1;
+  for offset = (1:100)
+    nerr = sum(xor(test_bits(offset:length(rx_bits)),rx_bits(1:length(rx_bits)+1-offset)));
+    bern = nerr/(bitcnt-offset);
+    if(bern < ber)
+      ox = offset;
+      best_nerr = nerr;
+    end
+    ber = min([ber bern]);
+  end
+  offset = ox;
+  printf("\ncoarse timing: %d nerr: %d\n", offset, best_nerr);
+
+  % Coherent BER theory
+  thrcoh = .5*erfc(sqrt(.5*EbNo));
+
+  % non-coherent BER theory calculation
+  % It was complicated, so I broke it up
+
+  ms = 2;
+  ns = (1:ms-1);
+  as = (-1).^(ns+1);
+  bs = (as./(ns+1));
+  
+  cs = ((ms-1)./ns);
+
+  ds = ns.*log2(ms);
+  es = ns+1;
+  fs = exp( -(ds./es)*EbNo );
+  
+  thrncoh = ((ms/2)/(ms-1)) * sum(bs.*((ms-1)./ns).*exp( -(ds./es)*EbNo ));
+
+endfunction
+
+
+function rxphi = fine_ex(timing_offset = 1,fsk_info)
+
+  rand('state',1); 
+  randn('state',1);
+
+  bitcnt = 12051;
+  test_bits = [zeros(1,100) rand(1,bitcnt)>.5]; %Random bits. Pad with zeros to prime the filters
+  t_vec = [0 0 1 1];
+  %test_bits = repmat(t_vec,1,ceil(24000/length(t_vec)));
+
+  states = yafsk_init(fsk_info);
+  Fs = states.Fs;
+  Rb = states.Rs; 
+  
+  tx = yafsk_mod(states,test_bits);
+
+  rx    = tx;
+  rx    = rx(timing_offset:length(rx));
+
+  [rx_bits states rxphi] = yafsk_demod_2a(states,rx);
+  ber = 1;
+  
+  %thing to account for offset from input data to output data
+  %No preamble detection yet
+  ox = 1;
+  for offset = (1:100)
+    nerr = sum(xor(rx_bits(offset:length(rx_bits)),test_bits(1:length(rx_bits)+1-offset)));
+    bern = nerr/(bitcnt-offset);
+    if(bern < ber)
+      ox = offset;
+      best_nerr = nerr;
+    end
+    ber = min([ber bern]);
+  end
+  offset = ox;
+  printf("\ncoarse timing: %d nerr: %d\n", offset, best_nerr);
+
+endfunction
+
+function yafsk_rx_phi
+  global fsk_setup
+  pkg load parallel
+  offrange = [1:2:101];
+
+  setups = repmat(fsk_setup,1,length(offrange));
+  phi = pararrayfun(nproc(),@nfbert_2,100,setups,1,offrange);
+  
+  close all;
+  figure(1);
+  clf;
+  plot(offrange,real(phi),offrange,imag(phi));
+  figure(2);
+  plotyy(offrange,angle(phi),offrange,abs(phi))
+endfunction
+