From: baobrien Date: Fri, 8 Apr 2016 01:17:54 +0000 (+0000) Subject: Octave unit tests now cover 2400a fsk configuration X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=49a77a9b742883a1b59767f8d29d455479ab00a9;p=freetel-svn-tracking.git Octave unit tests now cover 2400a fsk configuration git-svn-id: https://svn.code.sf.net/p/freetel/code@2775 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/octave/fsk_horus.m b/codec2-dev/octave/fsk_horus.m index 5b864543..323d4653 100644 --- a/codec2-dev/octave/fsk_horus.m +++ b/codec2-dev/octave/fsk_horus.m @@ -93,7 +93,7 @@ endfunction % Small, processing window nsym rather than nsym=Fs (1 second window) % Wider freq est limits -function states = fsk_horus_init_hbr(Fs,P,Rs,M=2,nsym=96) +function states = fsk_horus_init_hbr(Fs,P,Rs,M=2,nsym=48) assert((M==2) || (M==4), "Only M=2 and M=4 FSK supported"); states.M = M; @@ -106,7 +106,7 @@ function states = fsk_horus_init_hbr(Fs,P,Rs,M=2,nsym=96) states.nsym = N/Ts; % number of symbols in one processing frame states.nbit = states.nsym*states.bitspersymbol; % number of bits per processing frame - states.Ndft = 2.^ceil(log2(N))/2; % find nearest power of 2 for efficient FFT + states.Ndft = (2.^ceil(log2(N)))/2; % find nearest power of 2 for efficient FFT Nmem = states.Nmem = N+2*Ts; % two symbol memory in down converted signals to allow for timing adj @@ -125,7 +125,7 @@ function states = fsk_horus_init_hbr(Fs,P,Rs,M=2,nsym=96) states.fest_fmax = (Fs/2)-Rs; states.fest_fmin = Rs/2; - states.fest_min_spacing = Rs-(Rs/5); + states.fest_min_spacing = 2*(Rs-(Rs/5)); % BER stats diff --git a/codec2-dev/octave/tfsk.m b/codec2-dev/octave/tfsk.m index 841f076e..2013e9ce 100644 --- a/codec2-dev/octave/tfsk.m +++ b/codec2-dev/octave/tfsk.m @@ -550,9 +550,9 @@ function pass = test_drift_var(M) endfunction function pass = test_fsk_battery() - pass = test_mod_horuscfg_randbits; + %pass = test_mod_horuscfg_randbits; assert(pass) - pass = pass && test_mod_horuscfgm4_randbits; + %pass = pass && test_mod_horuscfgm4_randbits; assert(pass) pass = pass && test_drift_var(4); assert(pass) diff --git a/codec2-dev/octave/tfsk_2400a.m b/codec2-dev/octave/tfsk_2400a.m new file mode 100644 index 00000000..dd1a9ab8 --- /dev/null +++ b/codec2-dev/octave/tfsk_2400a.m @@ -0,0 +1,600 @@ +% tfsk.m +% Author: Brady O'Brien 8 January 2016 + + + +% Copyright 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, 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 . + + +% Octave script to check c port of fsk_horus against the fsk_horus.m +% +% [X] - Functions to wrap around fsk_mod and fsk_demod executables +% [X] - fsk_mod +% [X] - fsk_demod +% [X] - Functions to wrap around octave and c implementations, pass +% same dataset, compare outputs, and give clear go/no-go +% [X] - fsk_mod_test +% [X] - fsk_demod_test +% [X] - Port of run_sim and EbNodB curve test battery +% [X] - Extract and compare more parameters from demod +% [X] - Run some tests in parallel + +#{ + + FSK Modem automated test instructions: + + 1. Use cmake to build in debug mode to ensure unittest/tfsk is built: + + $ cd ~/codec2 + $ rm -Rf build_linux && mkdir build_linux + $ cd build_linux + $ cmake -DCMAKE_BUILD_TYPE=Debug .. + $ make + + 2 - Change tfsk_location below if required + 3 - Ensure Octave packages signal and parallel are installed + 4 - Start Octave and run tfsk.m. It will perform all tests automatically + +#} + +%tfsk executable path/file +global tfsk_location = '../build_linux/unittest/tfsk'; + + + +fsk_horus_as_a_lib = 1; % make sure calls to test functions at bottom are disabled +%fsk_horus_2fsk; +fsk_horus +pkg load signal; +pkg load parallel; +graphics_toolkit('gnuplot'); + + +global mod_pass_fail_maxdiff = 1e-3/5000; + +function mod = fsk_mod_c(Fs,Rs,f1,fsp,bits,M) + global tfsk_location; + %command to be run by system to launch the modulator + command = sprintf('%s MX %d %d %d %d %d fsk_mod_ut_bitvec fsk_mod_ut_modvec fsk_mod_ut_log.txt',tfsk_location,M,f1,fsp,Fs,Rs); + %save input bits into a file + bitvecfile = fopen('fsk_mod_ut_bitvec','wb+'); + fwrite(bitvecfile,bits,'uint8'); + fclose(bitvecfile); + + %run the modulator + system(command); + + modvecfile = fopen('fsk_mod_ut_modvec','rb'); + mod = fread(modvecfile,'single'); + fclose(modvecfile); + +endfunction + + +%Compare 2 vectors, fail if they are not close enough +function pass = vcompare(vc,voct,vname,tname,tol,pnum) + + %Get delta of vectors + dvec = abs(abs(vc)-abs(voct)); + + %Normalize difference + dvec = dvec ./ abs(max(abs(voct))+1e-8); + + maxdvec = abs(max(dvec)); + pass = maxdvec=states.nin + ninold = states.nin; + states = est_freq(states, modin(1:states.nin), states.M); + [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.f_dc(1,1:states.Nmem-Ts/P)]; + o_f2_dc = [o_f2_dc states.f_dc(2,1:states.Nmem-Ts/P)]; + o_f1_int = [o_f1_int states.f_int(1,:)]; + o_f2_int = [o_f2_int states.f_int(2,:)]; + o_EbNodB = [o_EbNodB states.EbNodB]; + o_ppm = [o_ppm states.ppm]; + o_rx_timing = [o_rx_timing states.rx_timing]; + o_norm_rx_timing = [o_norm_rx_timing states.norm_rx_timing]; + o_Sf = [o_Sf states.Sf']; + o_f1 = [o_f1 states.f(1)]; + o_f2 = [o_f1 states.f(2)]; + o_fest = [o_fest states.f]; + o_nin = [o_nin states.nin]; + if M==4 + o_f3_dc = [o_f3_dc states.f_dc(3,1:states.Nmem-Ts/P)]; + o_f4_dc = [o_f4_dc states.f_dc(4,1:states.Nmem-Ts/P)]; + o_f3_int = [o_f3_int states.f_int(3,:)]; + o_f4_int = [o_f4_int states.f_int(4,:)]; + o_f3 = [o_f1 states.f(3)]; + o_f4 = [o_f1 states.f(4)]; + end + end + + %close all + + pass = 1; + + pass = vcompare(o_Sf, t_fft_est(1:length(o_Sf)),'fft est',tname,1,1) && pass; + pass = vcompare(o_fest, t_f_est,'f est',tname,1,2) && pass; + pass = vcompare(o_rx_timing, t_rx_timing,'rx timing',tname,1,3) && pass; + + if M==4 + pass = vcompare(o_f3_dc, t_f3_dc, 'f3 dc', tname,1,4) && pass; + pass = vcompare(o_f4_dc, t_f4_dc, 'f4 dc', tname,1,5) && pass; + pass = vcompare(o_f3_int, t_f3_int, 'f3 int', tname,1,6) && pass; + pass = vcompare(o_f4_int, t_f4_int, 'f4 int', tname,1,7) && pass; + end + + pass = vcompare(o_f1_dc, t_f1_dc, 'f1 dc', tname,1,8) && pass; + pass = vcompare(o_f2_dc, t_f2_dc, 'f2 dc', tname,1,9) && pass; + pass = vcompare(o_f2_int, t_f2_int, 'f2 int', tname,1,10) && pass; + pass = vcompare(o_f1_int, t_f1_int, 'f1 int', tname,1,11) && pass; + + pass = vcompare(o_ppm , t_ppm, 'ppm', tname,1,12) && pass; + pass = vcompare(o_EbNodB, t_EbNodB,'EbNodB', tname,1,13) && pass; + pass = vcompare(o_nin, t_nin, 'nin', tname,1,14) && pass; + pass = vcompare(o_norm_rx_timing, t_norm_rx_timing,'norm rx timing',tname,1,15) && pass; + + + diffpass = sum(xor(obits,bits'))<4; + diffbits = sum(xor(obits,bits')); + + + printf('%d bit diff in test %s\n',diffbits,tname); + if diffpass==0 + printf('\n***bitcompare test failed test %s diff %d\n\n',tname,sum(xor(obits,bits'))) + figure(15) + plot(xor(obits,bits')) + title(sprintf('Bitcompare failure test %s',tname)) + end + + pass = pass && diffpass; + + assert(pass); + + test_stats.pass = pass; + test_stats.diff = sum(xor(obits,bits')); + test_stats.cbits = bits'; + test_stats.obits = obits; + +endfunction + +function [dmod,cmod,omod,pass] = fsk_mod_test(Fs,Rs,f1,fsp,bits,tname,M=2) + global mod_pass_fail_maxdiff; + %Run the C modulator + cmod = fsk_mod_c(Fs,Rs,f1,fsp,bits,M); + %Set up and run the octave modulator + states.M = M; + states = fsk_horus_init_hbr(Fs,10,Rs,M); + + states.ftx(1) = f1; + states.ftx(2) = f1+fsp; + + if states.M == 4 + states.ftx(3) = f1+fsp*2; + states.ftx(4) = f1+fsp*3; + end + + states.dA = [1 1 1 1]; + states.dF = 0; + omod = fsk_horus_mod(states,bits); + + dmod = cmod-omod; + pass = max(dmod)<(mod_pass_fail_maxdiff*length(dmod)); + if !pass + printf('Mod failed test %s!\n',tname); + end +endfunction + +% Random bit modulator test +% Pass random bits through the modulators and compare +function pass = test_mod_2400a_randbits + rand('state',1); + randn('state',1); + bits = rand(1,96000)>.5; + [dmod,cmod,omod,pass] = fsk_mod_test(48000,1200,1200,1200,bits,"mod 2400a randbits",4); + + if(!pass) + figure(1) + plot(dmod) + title("Difference between octave and C mod impl"); + end + +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 stats = tfsk_run_sim(test_frame_mode,EbNodB,timing_offset,fading,df,dA,M=2) + frames = 90; + %EbNodB = 10; + %timing_offset = 2.0; % see resample() for clock offset below + %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 + %dA = 1; % amplitude imbalance of tones (note this affects Eb so not a gd idea) + + more off + rand('state',1); + randn('state',1); + + % ---------------------------------------------------------------------- + + % sm2000 config ------------------------ + %states = fsk_horus_init(96000, 1200); + %states.f1_tx = 4000; + %states.f2_tx = 5200; + + if test_frame_mode == 2 + % 2400A config + states = fsk_horus_init_hbr(48000,10, 1200, M); + states.f1_tx = 1200; + states.f2_tx = 2400; + states.f3_tx = 3600; + states.f4_tx = 4800; + + end + + if test_frame_mode == 4 + % horus rtty config --------------------- + states = fsk_horus_init(8000, 100, M); + states.f1_tx = 1200; + states.f2_tx = 1600; + states.tx_bits_file = "horus_tx_bits_rtty.txt"; % Octave file of bits we FSK modulate + + end + + if test_frame_mode == 5 + % horus binary config --------------------- + states = fsk_horus_init(8000, 100, M); + states.f1_tx = 1200; + states.f2_tx = 1600; + %%%states.tx_bits_file = "horus_tx_bits_binary.txt"; % Octave file of bits we FSK modulate + states.tx_bits_file = "horus_payload_rtty.txt"; + end + + % ---------------------------------------------------------------------- + + states.verbose = 0;%x1; + N = states.N; + P = states.P; + Rs = states.Rs; + nsym = states.nsym; + Fs = states.Fs; + states.df = df; + states.dA = [dA dA dA dA]; + states.M = M; + + EbNo = 10^(EbNodB/10); + variance = states.Fs/(states.Rs*EbNo*states.bitspersymbol); + + % 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.nbit*(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) || (test_frame_mode == 5) + + % load up a horus msg from disk and modulate that + + test_frame = load(states.tx_bits_file); + 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 + + + + f1 = states.f1_tx; + fsp = states.f2_tx-f1 + states.dA = [dA dA dA dA]; + states.ftx(1) = f1; + states.ftx(2) = f1+fsp; + + if states.M == 4 + states.ftx(3) = f1+fsp*2; + states.ftx(4) = f1+fsp*3; + end + + tx = fsk_horus_mod(states, tx_bits); + + 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 + end + + noise = sqrt(variance)*randn(length(tx),1); + rx = tx + noise; + + test_name = sprintf("tfsk run sim EbNodB:%d frames:%d timing_offset:%d fading:%d df:%d",EbNodB,frames,timing_offset,fading,df); + tstats = fsk_demod_xt(Fs,Rs,states.f1_tx,fsp,rx,test_name,M); + printf("Test %s done\n",test_name); + + pass = tstats.pass; + obits = tstats.obits; + cbits = tstats.cbits; + + % Figure out BER of octave and C modems + bitcnt = length(tx_bits); + rx_bits = obits; + ber = 1; + ox = 1; + for offset = (1:100) + nerr = sum(xor(rx_bits(offset:length(rx_bits)),tx_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; + bero = ber; + ber = 1; + rx_bits = cbits; + ox = 1; + for offset = (1:100) + nerr = sum(xor(rx_bits(offset:length(rx_bits)),tx_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; + berc = ber; + stats.berc = berc; + stats.bero = bero; + % coherent BER theory calculation + + stats.thrcoh = .5*(M-1)*erfc(sqrt( (log2(M)/2) * EbNo )); + + % non-coherent BER theory calculation + % It was complicated, so I broke it up + + ms = M; + 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 )); + + stats.thrncoh = thrncoh; + stats.pass = pass; +endfunction + + +function pass = ebno_battery_test(timing_offset,fading,df,dA,M) + %Range of EbNodB over which to test + ebnodbrange = fliplr(5:2: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); + mv = repmat(M,1,ebnodbs); + %statv = pararrayfun(floor(1.25*nproc()),@tfsk_run_sim,modev,ebnodbrange,timingv,fadingv,dfv,dav,mv); + statv = arrayfun(@tfsk_run_sim,modev,ebnodbrange,timingv,fadingv,dfv,dav,mv); + + passv = zeros(1,length(statv)); + for ii=(1:length(statv)) + passv(ii)=statv(ii).pass; + end + + %All pass flags are '1' + pass = sum(passv)>=length(passv); + %and no tests died + pass = pass && length(passv)==ebnodbs; + passv; + assert(pass) +endfunction + +%Test with and without sample clock offset +function pass = test_timing_var(df,dA,M) + pass = ebno_battery_test(1,0,df,dA,M) + assert(pass) + pass = pass && ebno_battery_test(0,0,df,dA,M) + assert(pass) +endfunction + +%Test with and without 1 Hz/S freq drift +function pass = test_drift_var(M) + pass = test_timing_var(1,1,M) + assert(pass) + pass = pass && test_timing_var(0,1,M) + assert(pass) +endfunction + +function pass = test_fsk_battery() + pass = 1; + pass = pass && test_mod_2400a_randbits; + assert(pass) + pass = pass && test_drift_var(4); + assert(pass) + if pass + printf("***** All tests passed! *****\n"); + end +endfunction + +function plot_fsk_bers(M=2) + %Range of EbNodB over which to plot + ebnodbrange = (4:13); + + berc = ones(1,length(ebnodbrange)); + bero = ones(1,length(ebnodbrange)); + berinc = ones(1,length(ebnodbrange)); + beric = ones(1,length(ebnodbrange)); + ebnodbs = length(ebnodbrange) + mode = 2; + %Replication of other parameters for parcellfun + modev = repmat(mode,1,ebnodbs); + timingv = repmat(1,1,ebnodbs); + fadingv = repmat(0,1,ebnodbs); + dfv = repmat(1,1,ebnodbs); + dav = repmat(1,1,ebnodbs); + Mv = repmat(M,1,ebnodbs); + + + statv = pararrayfun(floor(nproc()),@tfsk_run_sim,modev,ebnodbrange,timingv,fadingv,dfv,dav,Mv); + %statv = arrayfun(@tfsk_run_sim,modev,ebnodbrange,timingv,fadingv,dfv,dav,Mv); + + for ii = (1:length(statv)) + stat = statv(ii); + berc(ii)=stat.berc; + bero(ii)=stat.bero; + berinc(ii)=stat.thrncoh; + beric(ii) = stat.thrcoh; + end + clf; + figure(M) + + semilogy(ebnodbrange, berinc,sprintf('r;%dFSK non-coherent theory;',M)) + hold on; + semilogy(ebnodbrange, beric ,sprintf('g;%dFSK coherent theory;',M)) + semilogy(ebnodbrange, bero ,sprintf('b;Octave fsk horus %dFSK Demod;',M)) + semilogy(ebnodbrange, berc,sprintf('+;C fsk horus %dFSK Demod;',M)) + hold off; + grid("minor"); + axis([min(ebnodbrange) max(ebnodbrange) 1E-5 1]) + legend("boxoff"); + xlabel("Eb/No (dB)"); + ylabel("Bit Error Rate (BER)") + +endfunction + + +test_fsk_battery +%plot_fsk_bers(2) +plot_fsk_bers(4) diff --git a/codec2-dev/src/freedv_api.c b/codec2-dev/src/freedv_api.c index d8706ee7..7eff7c1f 100644 --- a/codec2-dev/src/freedv_api.c +++ b/codec2-dev/src/freedv_api.c @@ -154,7 +154,7 @@ struct freedv *freedv_open(int mode) { if(f->deframer == NULL) return NULL; - f->fsk = fsk_create_hbr(48000,1200,40,4,1200,1200); + f->fsk = fsk_create_hbr(48000,1200,10,4,1200,1200); /* Note: fsk expects tx/rx bits as an array of uint8_ts, not ints */ f->tx_bits = (int*)malloc(f->fsk->Nbits*sizeof(uint8_t)); diff --git a/codec2-dev/src/fsk.c b/codec2-dev/src/fsk.c index 4b411194..be1eaeba 100644 --- a/codec2-dev/src/fsk.c +++ b/codec2-dev/src/fsk.c @@ -408,7 +408,7 @@ void fsk_demod_freq_est(struct FSK *fsk, float fsk_in[],float *freqs,int M){ for(j=0; j= max){ + if(fftout[j].i > max){ max = fftout[j].i; imax = j; } @@ -461,7 +461,7 @@ void fsk_demod_freq_est(struct FSK *fsk, float fsk_in[],float *freqs,int M){ f_min = f_min < 0 ? 0 : f_min; f_max = imax + f_zero; f_max = f_max > Ndft ? Ndft : f_max; - for(j=f_min+1; j