fsk.c test framework improvements
authorbaobrien <baobrien@01035d8c-6547-0410-b346-abe4f91aad63>
Thu, 21 Jan 2016 08:29:24 +0000 (08:29 +0000)
committerbaobrien <baobrien@01035d8c-6547-0410-b346-abe4f91aad63>
Thu, 21 Jan 2016 08:29:24 +0000 (08:29 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@2637 01035d8c-6547-0410-b346-abe4f91aad63

codec2-dev/octave/fsk_horus_2fsk.m
codec2-dev/octave/tfsk.m
codec2-dev/src/fsk.c
codec2-dev/src/fsk_demod.c
codec2-dev/src/modem_probe.c [new file with mode: 0644]
codec2-dev/src/modem_probe.h [new file with mode: 0644]
codec2-dev/unittest/CMakeLists.txt
codec2-dev/unittest/tfsk.c [new file with mode: 0644]

index d4818db6b5db1abbe3197afdb34815ba589a06ef..49637c21147dd5d4f17225fd231b51e0031c65bf 100644 (file)
@@ -274,8 +274,6 @@ function [rx_bits states] = fsk_horus_demod(states, sf)
   %%# same as sum of ((abs(f1_int)-abs(f2_int)).^2) .* exp(-j*w*(0:Np-1))
   x = ((abs(f1_int)-abs(f2_int)).^2) * exp(-j*w*(0:Np-1))';
   
-
-  x;
   norm_rx_timing = angle(x)/(2*pi);
   rx_timing = norm_rx_timing*P;
 
@@ -344,6 +342,7 @@ function [rx_bits states] = fsk_horus_demod(states, sf)
 
   x = abs(abs(f1_int_resample) - abs(f2_int_resample));
   states.EbNodB = 20*log10(1E-6+mean(x)/(1E-6+std(x)));
+ % printf("EbNodB %f\n",states.EbNodB);
 endfunction
 
 
index e42586efe137be6d7953e71ce28fc650b3975b3e..c57c1a31582e289aa90db7e73aa16c5e41df6188 100644 (file)
 
 fsk_horus_as_a_lib = 1; % make sure calls to test functions at bottom are disabled
 fsk_horus_2fsk;  
+pkg load signal;
 
 graphics_toolkit('gnuplot');
 
-fsk_demod_ex_file = "../build/src/fsk_demod";
 
 global mod_pass_fail_maxdiff = 1e-3/50000;
+global demod_pass_fail_maxdiff = .1;
 
 function mod = fsk_mod_c(Fs,Rs,f1,f2,bits)
     %Name of executable containing the modulator
@@ -46,12 +47,13 @@ endfunction
 
 function bits = fsk_demod_c(Fs,Rs,f1,f2,mod)
     %Name of executable containing the modulator
-    fsk_demod_ex_file = '../build/src/fsk_demod';
+    fsk_demod_ex_file = '../build/unittest/tfsk';
     modvecfilename = sprintf('fsk_demod_ut_modvec_%d',getpid());
     bitvecfilename = sprintf('fsk_demod_ut_bitvec_%d',getpid());
+    tvecfilename = sprintf('fsk_demod_ut_tracevec_%d.txt',getpid());
     
-    %command to be run by system to launch the modulator
-    command = sprintf('%s %d %d %d %d %s %s',fsk_demod_ex_file,Fs,Rs,f1,f2,modvecfilename,bitvecfilename);
+    %command to be run by system to launch the demod
+    command = sprintf('%s %d %d %d %d %s %s %s',fsk_demod_ex_file,Fs,Rs,f1,f2,modvecfilename,bitvecfilename,tvecfilename);
     
     %save modulated input into a file
     modvecfile = fopen(modvecfilename,'wb+');
@@ -69,6 +71,110 @@ function bits = fsk_demod_c(Fs,Rs,f1,f2,mod)
     %Clean up files
     delete(bitvecfilename);
     delete(modvecfilename);
+    %delete(tvecfilename);
+endfunction
+
+%Compare 2 vectors, fail if they are not close enough
+function pass = vcompare(va,vb,vname,tname)
+    global demod_pass_fail_maxdiff;
+    
+    %Get delta of vectors
+    dvec = abs(va)-abs(vb);     
+    
+    %Normalize difference
+    dvec = dvec ./ max(va);
+    
+    titlestr = sprintf('Diff between C and Octave of %s for %s',vname,tname);
+    pass = max(dvec)<(demod_pass_fail_maxdiff)
+    maxdvec = max(dvec)
+    %figure(12)
+    %title(titlestr)
+    %plot(abs(dvec))
+    
+    %figure(13)
+    %plot(abs(va))
+    
+    %figure(14)
+    %plot(abs(vb))
+    
+    printf('Comparing vectors %s in test %s\n',vname,tname);
+    assert(pass);
+    
+endfunction
+
+function test_stats = fsk_demod_xt(Fs,Rs,f1,f2,mod,tname)
+    global mod_pass_fail_maxdiff;
+    %Name of executable containing the modulator
+    fsk_demod_ex_file = '../build/unittest/tfsk';
+    modvecfilename = sprintf('fsk_demod_ut_modvec_%d',getpid());
+    bitvecfilename = sprintf('fsk_demod_ut_bitvec_%d',getpid());
+    tvecfilename = sprintf('fsk_demod_ut_tracevec_%d.txt',getpid());
+    
+    %command to be run by system to launch the demod
+    command = sprintf('%s %d %d %d %d %s %s %s',fsk_demod_ex_file,Fs,Rs,f1,f2,modvecfilename,bitvecfilename,tvecfilename);
+    
+    %save modulated input into a file
+    modvecfile = fopen(modvecfilename,'wb+');
+    fwrite(modvecfile,mod,'single');
+    fclose(modvecfile);
+    
+    %run the modulator
+    system(command);
+    
+    bitvecfile = fopen(bitvecfilename,'rb');
+    bits = fread(bitvecfile,'uint8');
+    fclose(bitvecfile);
+    bits = bits!=0;
+    
+    %Load test vec dump
+    load(tvecfilename);
+    
+    %Clean up files
+    delete(bitvecfilename);
+    delete(modvecfilename);
+    delete(tvecfilename);
+    
+    
+    o_f1_dc = [];
+    o_f2_dc = [];
+    o_f1_int = [];
+    o_f2_int = [];
+    o_EbNodB = [];
+    o_rx_timing = [];
+    %Run octave demod, dump some test vectors
+    states = fsk_horus_init(Fs,Rs);
+    states.f1_tx = f1;
+    states.f2_tx = f2;
+    states.dA = 1;
+    states.dF = 0;
+    modin = mod;
+    obits = [];
+    while length(modin)>=states.nin
+        ninold = states.nin;
+        [bitbuf,states] = fsk_horus_demod(states, modin(1:states.nin));
+        modin=modin(ninold+1:length(modin));
+        obits = [obits bitbuf];
+        
+        %Save other parameters
+        o_f1_dc = [o_f1_dc states.f1_dc];
+        o_f2_dc = [o_f2_dc states.f2_dc];
+        o_f1_int = [o_f1_int states.f1_int];
+        o_f2_int = [o_f2_int states.f2_int];
+        o_EbNodB = [o_EbNodB states.EbNodB];
+        o_rx_timing = [o_rx_timing states.rx_timing];
+    end
+    
+    pass = 1;
+    pass = and(pass,vcompare(o_rx_timing,t_rx_timing,'rx_timing',tname));
+    pass = and(pass,vcompare(o_f1_int,t_f1_int,'f1_int',tname));
+    pass = and(pass,vcompare(o_f2_int,t_f2_int,'f2_int',tname));
+     
+    diff = sum(xor(obits,bits'))
+    
+    pass = pass && diff<2;
+    test_stats.pass = pass;
+    assert(pass)
+    test_stats.diff = diff;
 endfunction
 
 function [dmod,cmod,omod,pass] = fsk_mod_test(Fs,Rs,f1,f2,bits,tname)
@@ -84,8 +190,6 @@ function [dmod,cmod,omod,pass] = fsk_mod_test(Fs,Rs,f1,f2,bits,tname)
     omod = fsk_horus_mod(states,bits);
     
     dmod = cmod-omod;
-    (mod_pass_fail_maxdiff*length(dmod))
-    max(abs(dmod))
     pass = max(dmod)<(mod_pass_fail_maxdiff*length(dmod))
     if !pass
         printf('Mod failed test %s!\n',tname);
@@ -128,6 +232,8 @@ endfunction
 % Random bit modulator test
 % Pass random bits through the modulators and compare
 function pass = test_mod_horuscfg_randbits
+    rand('state',1); 
+    randn('state',1);
     bits = rand(1,10000)>.5;
     [dmod,cmod,omod,pass] = fsk_mod_test(8000,100,1200,1600,bits,"mod horuscfg randbits");
     
@@ -135,12 +241,12 @@ function pass = test_mod_horuscfg_randbits
     plot(dmod)
     title("Difference between octave and C mod impl");
     
-    figure(2)
-    %plot((1:length(cmod)),cmod,(1:length(omod)),omod,'mod random bits')
 endfunction
 
 % A big ol' channel impairment tester
 % Shamlessly taken from fsk_horus
+% This throws some channel imparment or another at the C and octave modem so they 
+% may be compared.
 function pass = tfsk_run_sim(test_frame_mode,EbNodB,timing_offset,fading,df,dA)
   frames = 60;
   %EbNodB = 10;
@@ -160,6 +266,14 @@ function pass = tfsk_run_sim(test_frame_mode,EbNodB,timing_offset,fading,df,dA)
   %states = fsk_horus_init(96000, 1200);
   %states.f1_tx = 4000;
   %states.f2_tx = 5200;
+  
+  if test_frame_mode == 2
+    % horus rtty config ---------------------
+    states = fsk_horus_init(8000, 100);
+    states.f1_tx = 1200;
+    states.f2_tx = 1600;
+    
+  end
 
   if test_frame_mode == 4
     % horus rtty config ---------------------
@@ -228,8 +342,10 @@ function pass = tfsk_run_sim(test_frame_mode,EbNodB,timing_offset,fading,df,dA)
 
   tx = fsk_horus_mod(states, tx_bits);
 
-  tx = resample(tx, 1000, 1001); % simulated 1000ppm sample clock offset
-
+  if timing_offset
+    tx = resample(tx, 1000, 1001); % simulated 1000ppm sample clock offset
+  end
+  
   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
@@ -237,128 +353,18 @@ function pass = tfsk_run_sim(test_frame_mode,EbNodB,timing_offset,fading,df,dA)
 
   noise = sqrt(variance)*randn(length(tx),1);
   rx    = tx + noise;
-  %rx = real(rx);
-  %b1 = fir2(100, [0 4000 5200 48000]/48000, [1 1 0.5 0.5]);
-  %rx = filter(b1,1,rx);
-  %[b a] = cheby2(6,40,[3000 6000]/(Fs/2));
-  %rx = filter(b,a,rx);
-  %rx = sign(rx);
-  %rx(find (rx > 1)) = 1;
-  %rx(find (rx < -1)) = -1;
-
-  % dump simulated rx file
-  ftx=fopen("fsk_horus_100bd_binary.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);
-  x_log = [];
-  norm_rx_timing_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] = fsk_horus_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];
-
-    if test_frame_mode == 1
-       states = ber_counter(states, test_frame, rx_bits_buf);
-    end
-  end
   
   test_name = sprintf("tfsk_run_sim EbNodB:%d frames:%d timing_offset:%d fading:%d",EbNodB,frames,timing_offset,fading);
-  [obits cbits pass]=fsk_demod_test(Fs,Rs,states.f1_tx,states.f2_tx,rx,test_name); 
+  tstats = fsk_demod_xt(Fs,Rs,states.f1_tx,states.f2_tx,rx,test_name); 
   printf("Test %s done\n",test_name);
   
-  if test_frame_mode == 1
-    printf("frames: %d Tbits: %d Terrs: %d BER %4.3f\n", frames, states.Tbits,states. Terrs, states.Terrs/states.Tbits);
-  end
-
-  if test_frame_mode == 4
-    extract_and_print_rtty_packets(states, rx_bits_log, rx_bits_sd_log)
-  end
-
-  if test_frame_mode == 5
-    extract_and_decode_binary_packets(states, rx_bits_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(states.nerr_log)
-  title('num bit errors each frame')
-
-  figure(4)
-  clf
-  subplot(211)
-  plot(real(rx(1:Fs)))
-  title('rx signal at demod input')
-  subplot(212)
-  plot(abs(fft(rx(1:Fs))))
-
-  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')
-  if(pass)
-       close all;
-  end
+  pass = tstats.pass
+  
 endfunction
 
 function pass = ebno_battery_test(fading,df,dA)
     pass = 1
-    ebnodbrange = fliplr(3:13)
+    ebnodbrange = (2:13)
     timing_offset=0;
   
     npass=zeros(1,length(ebnodbrange));
@@ -369,27 +375,51 @@ function pass = ebno_battery_test(fading,df,dA)
     end
 endfunction
 
+function pass = ebno_battery_test(timing_offset,fading,df,dA)
+    %Range of EbNodB over which to test
+    ebnodbrange = fliplr(3:13);
+    ebnodbs = length(ebnodbrange);
+    
+    mode = 2;
+    %Replication of other parameters for parcellfun
+    modev   = repmat(mode,1,ebnodbs);
+    timingv = repmat(timing_offset,1,ebnodbs);
+    fadingv = repmat(fading,1,ebnodbs);
+    dfv     = repmat(df,1,ebnodbs);
+    dav     = repmat(dA,1,ebnodbs);
+    
+    passv = pararrayfun(floor(.7*nproc()),@tfsk_run_sim,modev,ebnodbrange,timingv,fadingv,dfv,dav);
+    %passv = arrayfun(@tfsk_run_sim,modev,ebnodbrange,timingv,fadingv,dfv,dav);
+    
+    pass = sum(passv)>=length(passv)
+    passv
+    assert(pass)
+endfunction
+
+%Test with and without channel fading
 function pass = test_fading_var(timing_offset,df,dA)
-    pass = 1
-    npass = ebno_battery_test(0,df,dA)
-    assert(npass)
-    pass = npass && pass
-    npass = ebno_battery_test(1,df,dA)
-    assert(npass)
-    pass = npass && pass
+    pass = ebno_battery_test(timing_offset,0,df,dA)
+    assert(pass)
+    pass = pass && ebno_battery_test(timing_offset,1,df,dA)
+    assert(pass)
 endfunction
 
+%Test with and without sample clock offset
+function pass = test_timing_var(df,dA)
+    pass = test_fading_var(0,df,dA)
+    assert(pass)
+    pass = pass && test_fading_var(1,df,dA)
+    assert(pass)
+endfunction
 
 function pass = test_fsk_battery()
-    pass = 1
-    npass = test_mod_horuscfg_randbits
-    pass = npass && pass
+    pass = test_mod_horuscfg_randbits
     assert(pass)
-    npass = test_fading_var(0,0,1);
-    pass = npass && pass
+    pass = pass && test_timing_var(0,1);
     assert(pass)
+    
     if pass
+        close all;
         printf("***** All tests passed! *****\n");
     end
 endfunction
index 8447edbdcc7eb587ea6959944bbe7a664ab6a376..9e04792bdc3a6333004eaea27e1fbe939b43b12e 100644 (file)
 /* P oversampling rate constant -- should probably be init-time configurable */
 #define ct_P 8
 
+/* Define this to enable EbNodB estimate */
+/* This needs square roots, may take more cpu time than it's worth */
+#define EST_EBNO
+
 /*---------------------------------------------------------------------------*\
 
                                INCLUDES
@@ -48,6 +52,7 @@
 #include "fsk.h"
 #include "comp_prim.h"
 #include "kiss_fftr.h"
+#include "modem_probe.h"
 
 /*---------------------------------------------------------------------------*\
 
@@ -323,17 +328,20 @@ void fsk_demod(struct FSK *fsk, uint8_t rx_bits[], float fsk_in[]){
     int nstash = fsk->nstash;
     COMP *f1_int, *f2_int;
     COMP t1,t2;
-    float phi1_d = fsk->phi1_d;
-    float phi2_d = fsk->phi2_d;
-    float phi_ft = 0; /* Phase of fine timing estimator */
+    COMP phi1_c = fsk->phi1_c;
+    COMP phi2_c = fsk->phi2_c;
+    COMP phi_ft;                 /* Phase of fine timing estimator */
     int nold = Nmem-nin;
-    float dphi1,dphi2,dphift;
+    COMP dphi1,dphi2;
+    COMP dphift;
     float f1,f2;
     float rx_timing,norm_rx_timing;//,old_norm_rx_timing;//,d_norm_rx_timing;
     int using_old_samps;
     float *sample_src;
     COMP *f1_intbuf,*f2_intbuf;
     
+    float meanebno,stdebno;
+    
     /* Estimate tone frequencies */
     fsk_demod_freq_est(fsk,fsk_in,&f1,&f2,&twist);
     
@@ -350,14 +358,9 @@ void fsk_demod(struct FSK *fsk, uint8_t rx_bits[], float fsk_in[]){
     f1_intbuf = (COMP*) alloca(sizeof(COMP)*Ts);
     f2_intbuf = (COMP*) alloca(sizeof(COMP)*Ts);
 
-    
-    /* Note: This would all be quite a bit faster with complex oscillators, like
-     * TX. */
-    /* TODO: change these to complex oscillators */
     /* Figure out how much to nudge each sample downmixer for every sample */
-    dphi1 = 2*M_PI*((float)(f1)/(float)(Fs));
-    dphi2 = 2*M_PI*((float)(f2)/(float)(Fs));
-
+    dphi1 = comp_exp_j(-2*M_PI*((float)(f1)/(float)(Fs)));
+    dphi2 = comp_exp_j(-2*M_PI*((float)(f2)/(float)(Fs)));
 
     dc_i = 0;
     cbuf_i = 0;
@@ -373,17 +376,21 @@ void fsk_demod(struct FSK *fsk, uint8_t rx_bits[], float fsk_in[]){
             using_old_samps = 0;
         }
         /* Downconvert and place into integration buffer */
-        f1_intbuf[dc_i]=fcmult(sample_src[dc_i],comp_exp_j(phi1_d));;
-        f2_intbuf[dc_i]=fcmult(sample_src[dc_i],comp_exp_j(phi2_d));;
+        f1_intbuf[dc_i]=fcmult(sample_src[dc_i],phi1_c);
+        f2_intbuf[dc_i]=fcmult(sample_src[dc_i],phi2_c);
 
         /* Spin downconversion phases */
-        phi1_d -= dphi1;
-        phi2_d -= dphi2;
-        if(phi1_d<0) phi1_d+=2*M_PI;
-        if(phi2_d<0) phi2_d+=2*M_PI;
+        phi1_c = cmult(phi1_c,dphi1);
+        phi2_c = cmult(phi2_c,dphi2);
     }
+    modem_probe_samp_c("t_f1_dc",&f1_intbuf[0],Ts-(Ts/P));
+    modem_probe_samp_c("t_f2_dc",&f2_intbuf[0],Ts-(Ts/P));
     cbuf_i = dc_i;
     
+    float f1_strs,f1_stis,f2_strs,f2_stis;
+    float f1_strc,f1_stic,f2_strc,f2_stic;
+    float y,t;
+    
     /* Integrate over Ts at offsets of Ts/P */
     for(i=0; i<(nsym+1)*P; i++){
         /* Downconvert and Place Ts/P samples in the integration buffers */
@@ -395,31 +402,61 @@ void fsk_demod(struct FSK *fsk, uint8_t rx_bits[], float fsk_in[]){
                 using_old_samps = 0;
             }
             /* Downconvert and place into integration buffer */
-            f1_intbuf[cbuf_i+j]=fcmult(sample_src[dc_i],comp_exp_j(phi1_d));;
-            f2_intbuf[cbuf_i+j]=fcmult(sample_src[dc_i],comp_exp_j(phi2_d));;
-
+            f1_intbuf[cbuf_i+j]=fcmult(sample_src[dc_i],phi1_c);
+            f2_intbuf[cbuf_i+j]=fcmult(sample_src[dc_i],phi2_c);
+            
             /* Spin downconversion phases */
-            phi1_d -= dphi1;
-            phi2_d -= dphi2;
-            if(phi1_d<0) phi1_d+=2*M_PI;
-            if(phi2_d<0) phi2_d+=2*M_PI;
+            phi1_c = cmult(phi1_c,dphi1);
+            phi2_c = cmult(phi2_c,dphi2);
         }
+        
+        /* Dump internal samples */
+        modem_probe_samp_c("t_f1_dc",&f1_intbuf[cbuf_i],Ts/P);
+        modem_probe_samp_c("t_f2_dc",&f2_intbuf[cbuf_i],Ts/P);
+        
         cbuf_i += Ts/P;
         if(cbuf_i>=Ts) cbuf_i = 0;
         
         /* Integrate over the integration buffers, save samples */
+        /* This uses Kahan summation to reduce floating point funnyness */
         t1 = t2 = comp0();
+        f1_strc = f1_stic = 0;
+        f2_strc = f2_stic = 0;
+        f1_strs = f1_stis = 0;
+        f2_strs = f2_stis = 0;
         for(j=0; j<Ts; j++){
-            t1 = cadd(t1,f1_intbuf[j]);
-            t2 = cadd(t2,f2_intbuf[j]);
+            y = f1_intbuf[j].real - f1_strc;
+            t = f1_strs + y;
+            f1_strc = (t - f1_strs) - y;
+            f1_strs = t;
+            
+            y = f1_intbuf[j].imag - f1_stic;
+            t = f1_stis + y;
+            f1_stic = (t - f1_stis) - y;
+            f1_stis = t;
+            
+            y = f2_intbuf[j].real - f2_strc;
+            t = f2_strs + y;
+            f2_strc = (t - f2_strs) - y;
+            f2_strs = t;
+            
+            y = f2_intbuf[j].imag - f2_stic;
+            t = f2_stis + y;
+            f2_stic = (t - f2_stis) - y;
+            f2_stis = t;
+            
+            //t1 = cadd(t1,f1_intbuf[j]);
+            //t2 = cadd(t2,f2_intbuf[j]);
         }
-        f1_int[i] = t1;
-        f2_int[i] = t2;
+        f1_int[i].real = f1_strs;
+        f1_int[i].imag = f1_stis;
+        f2_int[i].real = f2_strs;
+        f2_int[i].imag = f2_stis;
         
     }
 
-    fsk->phi1_d = phi1_d;
-    fsk->phi2_d = phi2_d;
+    fsk->phi1_c = comp_normalize(phi1_c);
+    fsk->phi2_c = comp_normalize(phi2_c);
 
     /* Stash samples away in the old sample buffer for the next round of bit getting */
     memcpy((void*)&(fsk->samp_old[0]),(void*)&(fsk_in[nin-nstash]),sizeof(float)*nstash);
@@ -429,20 +466,23 @@ void fsk_demod(struct FSK *fsk, uint8_t rx_bits[], float fsk_in[]){
      * exract angle */
      
     /* Figure out how much to spin the oscillator to extract magic spectral line */
-    dphift = 2*M_PI*((float)(Rs)/(float)(P*Rs));
+    dphift = comp_exp_j(-2*M_PI*((float)(Rs)/(float)(P*Rs)));
+    phi_ft.real = 1;
+    phi_ft.imag = 0;
     t1=comp0();
     for(i=0; i<(nsym+1)*P; i++){
         /* Get abs of f1_int[i] and f2_int[i] */
-        ft1 = sqrtf( (f1_int[i].real*f1_int[i].real) + (f1_int[i].imag*f1_int[i].imag) );
-        ft2 = sqrtf( (f2_int[i].real*f2_int[i].real) + (f2_int[i].imag*f2_int[i].imag) );
+        //ft1 = sqrtf( (f1_int[i].real*f1_int[i].real) + (f1_int[i].imag*f1_int[i].imag) );
+        //ft2 = sqrtf( (f2_int[i].real*f2_int[i].real) + (f2_int[i].imag*f2_int[i].imag) );
+        ft1 = cabsolute(f1_int[i]);
+        ft2 = cabsolute(f2_int[i]);
         /* Add and square 'em */
         ft1 = ft1-ft2;
         ft1 = ft1*ft1;
         /* Spin the oscillator for the magic line shift */
         /* Down shift and accumulate magic line */
-        t1 = cadd(t1,fcmult(ft1,comp_exp_j(phi_ft)));
-        phi_ft -= dphift;
-        if(phi_ft<0) phi_ft+=2*M_PI;
+        t1 = cadd(t1,fcmult(ft1,phi_ft));
+        phi_ft = cmult(phi_ft,dphift);
     }
     /* Get the magic angle */
     norm_rx_timing =  -atan2f(t1.imag,t1.real)/(2*M_PI);
@@ -465,10 +505,16 @@ void fsk_demod(struct FSK *fsk, uint8_t rx_bits[], float fsk_in[]){
     else
         fsk->nin = N;
     
+    modem_probe_samp_f("t_norm_rx_timing",&(norm_rx_timing),1);;
+    
     /* Re-sample the integrators with linear interpolation magic */
     int low_sample = (int)floorf(rx_timing);
     float fract = rx_timing - (float)low_sample;
     int high_sample = (int)ceilf(rx_timing);
+    
+    #ifdef EST_EBNO
+    meanebno = 0;
+    #endif
   
     /* FINALLY, THE BITS */
     /* also, resample f1_int,f2_int */
@@ -479,12 +525,52 @@ void fsk_demod(struct FSK *fsk, uint8_t rx_bits[], float fsk_in[]){
         t2 =         fcmult(1-fract,f2_int[st+ low_sample]);
         t2 = cadd(t2,fcmult(  fract,f2_int[st+high_sample]));
         
+        /* Accumulate resampled int magnitude for EbNodB estimation */
+        #ifdef EST_EBNO
+        ft1 = sqrtf(t1.real*t1.real + t1.imag*t1.imag);
+        ft2 = sqrtf(t2.real*t2.real + t2.imag*t2.imag);
+        ft1 = fabsf(ft1-ft2);
+        meanebno += ft1;
+        #endif
+        
         /* THE BIT! */
         rx_bits[i] = comp_mag_gt(t2,t1)?1:0;
         /* Soft output goes here */
     }
     
-    printf("rx_timing: %3.2f low_sample: %d high_sample: %d fract: %3.3f nin_next: %d\n", rx_timing, low_sample, high_sample, fract, fsk->nin);
+    #ifdef EST_EBNO
+    /* Calculate mean for EbNodB estimation */
+    meanebno = meanebno/(float)nsym;
+    stdebno = 0;
+    /* Go back through the data and figure the std dev */
+    for(i=0; i<nsym; i++){
+        int st = (i+1)*P;
+        t1 =         fcmult(1-fract,f1_int[st+ low_sample]);
+        t1 = cadd(t1,fcmult(  fract,f1_int[st+high_sample]));
+        t2 =         fcmult(1-fract,f2_int[st+ low_sample]);
+        t2 = cadd(t2,fcmult(  fract,f2_int[st+high_sample]));
+        
+        /* Accumulate resampled int magnitude for EbNodB estimation */
+        ft1 = sqrtf(t1.real*t1.real + t1.imag*t1.imag);
+        ft2 = sqrtf(t2.real*t2.real + t2.imag*t2.imag);
+        ft1 = fabsf(ft1-ft2);
+        ft2 = abs(meanebno-ft1);
+        stdebno += ft2*ft2;
+    }
+    /* Finish figuring std. dev. */
+    stdebno = sqrtf(stdebno/(float)(nsym-1));
+    fsk->EbNodB = 20*log10f((1e-6+meanebno)/(1e-6+stdebno));
+    #else
+    fsk->EbNodB = 0;
+    #endif
+    
+    /* Dump some internal samples */
+    modem_probe_samp_f("t_EbNodB",&(fsk->EbNodB),1);
+    modem_probe_samp_f("t_f1",&f1,1);
+    modem_probe_samp_f("t_f2",&f2,1);
+    modem_probe_samp_c("t_f1_int",f1_int,(nsym+1)*P);
+    modem_probe_samp_c("t_f2_int",f2_int,(nsym+1)*P);
+    modem_probe_samp_f("t_rx_timing",&(rx_timing),1);;
 }
 
 
index a8eaa3102660915c89ce4a410e3f0430db953e61..29e35c3c8e5bc4504c6560a35da7cd8c4f9f58f1 100644 (file)
@@ -29,6 +29,7 @@
 
 #include <stdio.h>
 #include "fsk.h"
+#include "modem_probe.h"
 
 int main(int argc,char *argv[]){
     struct FSK *fsk;
@@ -38,7 +39,7 @@ int main(int argc,char *argv[]){
     float *modbuf;
     
     if(argc<7){
-        printf("Usage: %s samplerate bitrate f1 f2 rffile bitfile\n",argv[0]);
+        printf("Usage: %s samplerate bitrate f1 f2 rffile bitfile [probefile]\n",argv[0]);
         exit(1);
     }
     
@@ -52,6 +53,9 @@ int main(int argc,char *argv[]){
     fin = fopen(argv[5],"r");
     fout = fopen(argv[6],"w");
     
+    if(argc>7)
+               modem_probe_init("fsk2",argv[7]);
+       
     /* set up FSK */
     fsk = fsk_create(Fs,Rs,f1,f2);
     
@@ -70,6 +74,7 @@ int main(int argc,char *argv[]){
         fwrite(bitbuf,sizeof(uint8_t),fsk->Nsym,fout);
     }
     
+    modem_probe_close();
     cleanup:
     fclose(fin);
     fclose(fout);
diff --git a/codec2-dev/src/modem_probe.c b/codec2-dev/src/modem_probe.c
new file mode 100644 (file)
index 0000000..612d335
--- /dev/null
@@ -0,0 +1,239 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: modem_probe.c
+  AUTHOR......: Brady O'Brien
+  DATE CREATED: 9 January 2016
+
+  Library to easily extract debug traces from modems during development
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2016 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program is
+  distributed in the hope that it will be useful, but WITHOUT ANY
+  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+  License for more details.
+
+  You should have received a copy of the GNU Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include "comp.h"
+#include "octave.h"
+
+#define TRACE_I 1
+#define TRACE_F 2
+#define TRACE_C 3
+
+
+typedef struct probe_trace_info_s probe_trace_info;
+typedef struct datlink_s datlink;
+
+struct datlink_s{
+       void * data;
+       size_t len;
+       datlink * next;
+};
+
+struct probe_trace_info_s{
+       int type;
+       char name[255];
+       datlink * data;
+       datlink * last;
+       probe_trace_info *next;
+};
+
+static char *run = NULL;
+static char *mod = NULL;
+static probe_trace_info *first_trace = NULL;
+
+/* Init the probing library */
+void modem_probe_init_int(char *modname, char *runname){
+       mod = malloc((strlen(modname)+1)*sizeof(char));
+       run = malloc((strlen(runname)+1)*sizeof(char));
+       strcpy(run,runname);
+       strcpy(mod,modname);
+}
+
+/* 
+ * Gather the data stored in the linked list into a single blob,
+ * freeing links and buffers as it goes
+ */
+void * gather_data(datlink * d,size_t * len){
+       size_t size = 0;
+       datlink * cur = d;
+       datlink * next;
+       while(cur!=NULL){
+               size += d->len;
+               cur = cur->next;
+       }
+       cur = d;
+       size_t i = 0;
+       void * newbuf = malloc(size);
+       
+       while(cur!=NULL){
+               memcpy(newbuf+i,cur->data,cur->len);
+               i += cur->len;
+               free(cur->data);
+               next = cur->next;
+               free(cur);
+               cur = next;
+       }
+       *len = size;
+       return newbuf;
+}
+
+/* Dump all of the traces into a nice octave-able dump file */
+void modem_probe_close_int(){
+       if(run==NULL)
+               return;
+       
+       probe_trace_info *cur,*next;
+       cur = first_trace;
+       FILE * dumpfile = fopen(run,"w");
+       void * dbuf;
+       size_t len;
+       
+       while(cur != NULL){
+               dbuf = gather_data(cur->data,&len);
+               switch(cur->type){
+                       case TRACE_I:
+                               octave_save_int(dumpfile,cur->name,(int32_t*)dbuf,1,len/sizeof(int32_t));
+                               break;
+                       case TRACE_F:
+                               octave_save_float(dumpfile,cur->name,(float*)dbuf,1,len/sizeof(float),10);
+                               break;
+                       case TRACE_C:
+                               octave_save_complex(dumpfile,cur->name,(COMP*)dbuf,1,len/sizeof(COMP),10);
+                               break;
+               }
+               next = cur->next;
+               free(cur);
+               free(dbuf);
+               cur = next;
+       }
+       
+       fclose(dumpfile);
+       free(run);
+       free(mod);
+}
+
+/* Look up or create a trace by name */
+probe_trace_info * modem_probe_get_trace(char * tracename){
+       probe_trace_info *cur,*npti;
+       
+       /* Make sure probe session is open */
+       if(run==NULL)
+               return NULL;
+       
+       cur = first_trace;
+       /* Walk through list, find trace with matching name */
+       while(cur != NULL){
+               /* We got one! */
+               if(strcmp( cur->name, tracename) == 0){
+                       return cur;
+               }
+               cur = cur->next;
+       }
+       /* None found, open a new trace */
+       
+       npti = (probe_trace_info *) malloc(sizeof(probe_trace_info));
+       npti->next = first_trace;
+       npti->data = NULL;
+       npti->last = NULL;
+       strcpy(npti->name,tracename);
+       first_trace = npti;
+       
+       return npti;
+       
+}
+
+void modem_probe_samp_i_int(char * tracename,int32_t samp[],size_t cnt){
+       probe_trace_info *pti;
+       datlink *ndat;
+       
+       pti = modem_probe_get_trace(tracename);
+       if(pti == NULL)
+               return;
+       
+       pti->type = TRACE_I;
+       
+       ndat = (datlink*) malloc(sizeof(datlink));
+       ndat->data = malloc(sizeof(int32_t)*cnt);
+       
+       ndat->len = cnt*sizeof(int32_t);
+       ndat->next = NULL;
+       memcpy(ndat->data,(void*)&(samp[0]),sizeof(int32_t)*cnt);
+       
+       if(pti->last!=NULL){
+               pti->last->next = ndat;
+               pti->last = ndat;
+       } else {
+               pti->data = ndat;
+               pti->last = ndat;
+       }
+       
+}
+
+void modem_probe_samp_f_int(char * tracename,float samp[],size_t cnt){
+       probe_trace_info *pti;
+       datlink *ndat;
+       
+       pti = modem_probe_get_trace(tracename);
+       if(pti == NULL)
+               return;
+       
+       pti->type = TRACE_F;
+       
+       ndat = (datlink*) malloc(sizeof(datlink));
+       ndat->data = malloc(sizeof(float)*cnt);
+       
+       ndat->len = cnt*sizeof(float);
+       ndat->next = NULL;
+       memcpy(ndat->data,(void*)&(samp[0]),sizeof(float)*cnt);
+       
+       if(pti->last!=NULL){
+               pti->last->next = ndat;
+               pti->last = ndat;
+       } else {
+               pti->data = ndat;
+               pti->last = ndat;
+       }
+}
+       
+void modem_probe_samp_c_int(char * tracename,COMP samp[],size_t cnt){
+       probe_trace_info *pti;
+       datlink *ndat;
+       
+       pti = modem_probe_get_trace(tracename);
+       if(pti == NULL)
+               return;
+       
+       pti->type = TRACE_C;
+       
+       ndat = (datlink*) malloc(sizeof(datlink));
+       ndat->data = malloc(sizeof(COMP)*cnt);
+       
+       ndat->len = cnt*sizeof(COMP);
+       ndat->next = NULL;
+       memcpy(ndat->data,(void*)&(samp[0]),sizeof(COMP)*cnt);
+       
+       if(pti->last!=NULL){
+               pti->last->next = ndat;
+               pti->last = ndat;
+       } else {
+               pti->data = ndat;
+               pti->last = ndat;
+       }
+}
diff --git a/codec2-dev/src/modem_probe.h b/codec2-dev/src/modem_probe.h
new file mode 100644 (file)
index 0000000..320acaa
--- /dev/null
@@ -0,0 +1,90 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: modem_probe.h
+  AUTHOR......: Brady O'Brien
+  DATE CREATED: 9 January 2016
+
+  Library to easily extract debug traces from modems during development
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2016 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program is
+  distributed in the hope that it will be useful, but WITHOUT ANY
+  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+  License for more details.
+
+  You should have received a copy of the GNU Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef __MODEMPROBE_H
+#define __MODEMPROBE_H
+#include <stdint.h>
+#include <stdlib.h>
+#include "comp.h"
+
+
+#ifdef MODEMPROBE_ENABLE
+
+void modem_probe_init_int(char *modname, char *runname);
+void modem_probe_close_int();
+
+void modem_probe_samp_i_int(char * tracename,int samp[],size_t cnt);
+void modem_probe_samp_f_int(char * tracename,float samp[],size_t cnt);
+void modem_probe_samp_c_int(char * tracename,COMP samp[],size_t cnt);
+
+
+static inline void modem_probe_init(char *modname,char *runname){
+       modem_probe_init_int(modname,runname);
+}
+
+static inline void modem_probe_close(){
+       modem_probe_close_int();
+}
+
+static inline void modem_probe_samp_i(char *tracename,int samp[],size_t cnt){
+       modem_probe_samp_i_int(tracename,samp,cnt);
+}
+
+
+static inline void modem_probe_samp_f(char *tracename,float samp[],size_t cnt){
+       modem_probe_samp_f_int(tracename,samp,cnt);
+}      
+
+static inline void modem_probe_samp_c(char *tracename,COMP samp[],size_t cnt){
+       modem_probe_samp_c_int(tracename,samp,cnt);
+}
+
+#else
+
+static inline void modem_probe_init(char *modname,char *runname){
+       return;
+}
+
+static inline void modem_probe_close(){
+       return;
+}
+
+static inline void modem_probe_samp_i(char *name,int samp[],size_t sampcnt){
+       return;
+}
+
+static inline void modem_probe_samp_f(char *name,float samp[],size_t cnt){
+       return;
+}
+
+static inline void modem_probe_samp_c(char *name,COMP samp[],size_t cnt){
+       return;
+}
+
+#endif
+
+#endif
index ad4ba7e283648d6f695dfcb888fcc9e99fe2c74b..accc7a83f99bb2e0dc3d3752dcb4f6db4b9b57f4 100644 (file)
@@ -63,6 +63,9 @@ target_link_libraries(t16_8 codec2)
 add_executable(t16_8_short t16_8_short.c ../src/fdmdv.c ../src/kiss_fft.c)
 target_link_libraries(t16_8_short codec2)
 
+add_executable(tfsk tfsk.c ../src/kiss_fft.c ../src/kiss_fftr.c ../src/octave.c ../src/modem_probe.c)
+target_link_libraries(tfsk m)
+
 #add_executable(t48_8 t48_8.c ../src/fdmdv.c ../src/kiss_fft.c)
 #target_link_libraries(t48_8 codec2)
 
diff --git a/codec2-dev/unittest/tfsk.c b/codec2-dev/unittest/tfsk.c
new file mode 100644 (file)
index 0000000..ae85242
--- /dev/null
@@ -0,0 +1,90 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: tfsk.c
+  AUTHOR......: Brady O'Brien
+  DATE CREATED: 20 January 2016
+
+  C test driver for fsk_mod and fsk_demod in fsk.c. Reads a file with input
+  bits/rf and spits out modulated/demoduladed samples and a dump of internal
+  state 
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2016 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program is
+  distributed in the hope that it will be useful, but WITHOUT ANY
+  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+  License for more details.
+
+  You should have received a copy of the GNU Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#define MODEMPROBE_ENABLE
+
+#include "modem_probe.h"
+#include <stdio.h>
+
+/* Note: This is a dirty hack to force fsk.c to compile with modem probing enabled */
+#include "fsk.c"
+
+
+int main(int argc,char *argv[]){
+    struct FSK *fsk;
+    int Fs,Rs,f1,f2;
+    FILE *fin,*fout;
+    uint8_t *bitbuf;
+    float *modbuf;
+    
+    if(argc<7){
+        printf("Usage: %s samplerate bitrate f1 f2 infile outfile [probefile]\n",argv[0]);
+        exit(1);
+    }
+    
+    /* Extract parameters */
+    Fs = atoi(argv[1]);
+    Rs = atoi(argv[2]);
+    f1 = atoi(argv[3]);
+    f2 = atoi(argv[4]);
+    
+    /* Open files */
+    fin = fopen(argv[5],"r");
+    fout = fopen(argv[6],"w");
+    
+    if(argc>7)
+               modem_probe_init("fsk2",argv[7]);
+       
+    /* set up FSK */
+    fsk = fsk_create(Fs,Rs,f1,f2);
+    
+    if(fin==NULL || fout==NULL || fsk==NULL){
+        printf("Couldn't open test vector files\n");
+        goto cleanup;
+    }
+    
+    /* allocate buffers for processing */
+    bitbuf = (uint8_t*)alloca(sizeof(uint8_t)*fsk->Nsym);
+    modbuf = (float*)alloca(sizeof(float)*(fsk->N+fsk->Ts*2));
+    
+    /* Demodulate! */
+    while( fread(modbuf,sizeof(float),fsk_nin(fsk),fin) == fsk_nin(fsk) ){
+        fsk_demod(fsk,bitbuf,modbuf);
+        fwrite(bitbuf,sizeof(uint8_t),fsk->Nsym,fout);
+    }
+    
+    modem_probe_close();
+    cleanup:
+    fclose(fin);
+    fclose(fout);
+    fsk_destroy(fsk);
+    exit(0);
+}
+