aboutsummaryrefslogtreecommitdiff
path: root/octave/tfsk_2400a.m
diff options
context:
space:
mode:
Diffstat (limited to 'octave/tfsk_2400a.m')
-rw-r--r--octave/tfsk_2400a.m648
1 files changed, 0 insertions, 648 deletions
diff --git a/octave/tfsk_2400a.m b/octave/tfsk_2400a.m
deleted file mode 100644
index 87a1e25..0000000
--- a/octave/tfsk_2400a.m
+++ /dev/null
@@ -1,648 +0,0 @@
-% 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.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/>.
-
-
-% 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_2400a.m. It will perform all tests automatically
-
-#}
-
-
-
-%tfsk executable path/file
-global tfsk_location = '../build_linux/unittest/tfsk';
-
-%Set to 1 for verbose printouts
-global print_verbose = 0;
-
-
-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)
- global print_verbose;
- %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<tol;
-
- if print_verbose == 1
- printf(' Comparing vectors %s in test %s. Diff is %f\n',vname,tname,maxdvec);
- end
-
- if pass == 0
- printf('\n*** vcompare failed %s in test %s. Diff: %f Tol: %f\n\n',vname,tname,maxdvec,tol);
-
- titlestr = sprintf('Diff between C and Octave of %s for %s',vname,tname)
- figure(10+pnum*2)
- plot(abs(dvec))
- title(titlestr)
-
- figure(11+pnum*2)
- plot((1:length(vc)),abs(vc),(1:length(voct)),abs(voct))
-
- end
-
-endfunction
-
-function test_stats = fsk_demod_xt(Fs,Rs,f1,fsp,mod,tname,M=2)
- global tfsk_location;
- global print_verbose;
- %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 DX %d %d %d %d %d %s %s %s',tfsk_location,M,f1,fsp,Fs,Rs,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_f3_dc = [];
- o_f4_dc = [];
- o_f1_int = [];
- o_f2_int = [];
- o_f3_int = [];
- o_f4_int = [];
- o_f1 = [];
- o_f2 = [];
- o_f3 = [];
- o_f4 = [];
- o_EbNodB = [];
- o_ppm = [];
- o_Sf = [];
- o_fest = [];
- o_rx_timing = [];
- o_norm_rx_timing = [];
- o_nin = [];
- %Run octave demod, dump some test vectors
- states = fsk_horus_init_hbr(Fs,10,Rs,M);
- Ts = states.Ts;
- P = states.P;
- states.ftx(1) = f1;
- states.ftx(2) = f1+fsp;
- states.ftx(3) = f1+fsp*2;
- states.ftx(4) = f1+fsp*3;
- states.dA = 1;
- states.dF = 0;
- modin = mod;
- obits = [];
- while length(modin)>=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'))<5;
- diffbits = sum(xor(obits,bits'));
-
- if print_verbose == 1
- printf('%d bit diff in test %s\n',diffbits,tname);
- end
- 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)
- global print_verbose;
- frames = 190;
- %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',10);
- randn('state',10);
-
- % ----------------------------------------------------------------------
-
- % 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;
- states.ftx(1) = 1200;
- states.ftx(2) = 2400;
- states.ftx(3) = 3600;
- states.ftx(4) = 4800;
-
- end
-
- if test_frame_mode == 4
- % horus rtty 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;
- states.ftx(1) = 1200;
- states.ftx(2) = 2400;
- states.ftx(3) = 3600;
- states.ftx(4) = 4800;
-
- 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_hbr(48000,10, 1200, M);
- states.f1_tx = 1200;
- states.f2_tx = 2400;
- states.f3_tx = 3600;
- states.f4_tx = 4800;
- states.ftx(1) = 1200;
- states.ftx(2) = 2400;
- states.ftx(3) = 3600;
- states.ftx(4) = 4800;
- %%%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 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);
-
- pass = tstats.pass;
- obits = tstats.obits;
- cbits = tstats.cbits;
- stats.name = test_name;
-
- if tstats.pass
- printf("Test %s passed\n",test_name);
- else
- printf("Test %s failed\n",test_name);
- end
-
- % 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
-
- if print_verbose == 1
- printf("C BER: %f Oct BER: %f Test %s\n",berc,bero,test_name);
- end
-
- 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;
- if statv(ii).pass
- printf("Test %s passed\n",statv(ii).name);
- else
- printf("Test %s failed\n",statv(ii).name);
- end
-
- 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
-
-
-xpass = test_fsk_battery
-%plot_fsk_bers(2)
-plot_fsk_bers(4)
-
-if xpass
- printf("***** All tests passed! *****\n");
-else
- printf("***** Some test failed! Look back through output to find failed test *****\n");
-end
generated by cgit v1.2.3 (git 2.46.0) at 2026-01-26 22:07:40 +0000