From 3d2fba10ae3190797cfed9dfdb4d1fdc2306eaad Mon Sep 17 00:00:00 2001 From: baobrien Date: Thu, 21 Jan 2016 08:29:24 +0000 Subject: [PATCH] fsk.c test framework improvements git-svn-id: https://svn.code.sf.net/p/freetel/code@2637 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/fsk_horus_2fsk.m | 3 +- codec2-dev/octave/tfsk.m | 304 ++++++++++++++++------------- codec2-dev/src/fsk.c | 160 +++++++++++---- codec2-dev/src/fsk_demod.c | 7 +- codec2-dev/src/modem_probe.c | 239 +++++++++++++++++++++++ codec2-dev/src/modem_probe.h | 90 +++++++++ codec2-dev/unittest/CMakeLists.txt | 3 + codec2-dev/unittest/tfsk.c | 90 +++++++++ 8 files changed, 719 insertions(+), 177 deletions(-) create mode 100644 codec2-dev/src/modem_probe.c create mode 100644 codec2-dev/src/modem_probe.h create mode 100644 codec2-dev/unittest/tfsk.c diff --git a/codec2-dev/octave/fsk_horus_2fsk.m b/codec2-dev/octave/fsk_horus_2fsk.m index d4818db6..49637c21 100644 --- a/codec2-dev/octave/fsk_horus_2fsk.m +++ b/codec2-dev/octave/fsk_horus_2fsk.m @@ -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 diff --git a/codec2-dev/octave/tfsk.m b/codec2-dev/octave/tfsk.m index e42586ef..c57c1a31 100644 --- a/codec2-dev/octave/tfsk.m +++ b/codec2-dev/octave/tfsk.m @@ -16,12 +16,13 @@ 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 diff --git a/codec2-dev/src/fsk.c b/codec2-dev/src/fsk.c index 8447edbd..9e04792b 100644 --- a/codec2-dev/src/fsk.c +++ b/codec2-dev/src/fsk.c @@ -34,6 +34,10 @@ /* 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; jphi1_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; iEbNodB = 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);; } diff --git a/codec2-dev/src/fsk_demod.c b/codec2-dev/src/fsk_demod.c index a8eaa310..29e35c3c 100644 --- a/codec2-dev/src/fsk_demod.c +++ b/codec2-dev/src/fsk_demod.c @@ -29,6 +29,7 @@ #include #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 index 00000000..612d3350 --- /dev/null +++ b/codec2-dev/src/modem_probe.c @@ -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 . +*/ + +#include +#include +#include +#include +#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 index 00000000..320acaae --- /dev/null +++ b/codec2-dev/src/modem_probe.h @@ -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 . +*/ + +#ifndef __MODEMPROBE_H +#define __MODEMPROBE_H +#include +#include +#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 diff --git a/codec2-dev/unittest/CMakeLists.txt b/codec2-dev/unittest/CMakeLists.txt index ad4ba7e2..accc7a83 100644 --- a/codec2-dev/unittest/CMakeLists.txt +++ b/codec2-dev/unittest/CMakeLists.txt @@ -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 index 00000000..ae852423 --- /dev/null +++ b/codec2-dev/unittest/tfsk.c @@ -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 . +*/ + + +#define MODEMPROBE_ENABLE + +#include "modem_probe.h" +#include + +/* 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); +} + -- 2.25.1