1% tfsk.m 2% Brady O'Brien 8 January 2016 3% David Rowe May 2020 4% 5% Automatic testing of C port of FSK modem by comparing to reference 6% Octave version. Currently just a subset of tests enabled in order to 7% run in a reasonable amount of time as ctests, but still trapping any 8% bit-rot. 9 10#{ 11 12 FSK Modem automated test instructions: 13 14 1. Use cmake to build in debug mode to ensure unittest/tfsk is built: 15 16 $ cd ~/codec2 17 $ rm -Rf build_linux && mkdir build_linux 18 $ cd build_linux 19 $ cmake -DCMAKE_BUILD_TYPE=Debug .. 20 $ make 21 22 2 - Change tfsk_location below if required 23 3 - Ensure Octave packages are installed 24 4 - Start Octave and run tfsk.m. It will perform all tests automatically 25 26#} 27 28% tfsk executable path/file 29if getenv("PATH_TO_TFSK") 30 global tfsk_location = getenv("PATH_TO_TFSK") 31 printf("setting tfsk_location from env var: %s\n", getenv("PATH_TO_TFSK")); 32else 33 global tfsk_location = '../build_linux/unittest/tfsk'; 34end 35 36% Set to 1 for verbose printouts 37global print_verbose = 0; 38global mod_pass_fail_maxdiff = 1e-3/5000; 39 40fsk_horus_as_a_lib=1; 41fsk_horus; 42pkg load signal; 43% not needed unless parallel tests running 44%pkg load parallel; 45graphics_toolkit('gnuplot'); 46 47function print_result(test_name, result) 48 printf("%s", test_name); 49 for i=1:(40-length(test_name)) 50 printf("."); 51 end 52 printf(": %s\n", result); 53end 54 55function mod = fsk_mod_c(Fs,Rs,f1,fsp,bits,M) 56 global tfsk_location; 57 %command to be run by system to launch the modulator 58 command = sprintf('%s M %d %d %d %d %d 0 fsk_mod_ut_bitvec fsk_mod_ut_modvec fsk_mod_ut_log.txt',tfsk_location,M,f1,fsp,Fs,Rs); 59 %save input bits into a file 60 bitvecfile = fopen('fsk_mod_ut_bitvec','wb+'); 61 fwrite(bitvecfile,bits,'uint8'); 62 fclose(bitvecfile); 63 64 %run the modulator 65 system(command); 66 67 modvecfile = fopen('fsk_mod_ut_modvec','rb'); 68 mod = fread(modvecfile,'single'); 69 fclose(modvecfile); 70 71endfunction 72 73 74%Compare 2 vectors, fail if they are not close enough 75function pass = vcompare(vc,voct,vname,tname,tol,pnum) 76 global print_verbose; 77 %Get delta of vectors 78 dvec = abs(abs(vc - voct)); 79 80 %Normalize difference 81 dvec = dvec ./ abs(max(abs(voct))+1e-8); 82 83 maxdvec = abs(max(dvec)); 84 pass = maxdvec<tol; 85 if print_verbose == 1 86 printf(' Comparing vectors %s in test %s. Diff is %f\n',vname,tname,maxdvec); 87 end 88 if pass == 0 89 printf('\n*** vcompare failed %s in test %s. Diff: %f Tol: %f\n\n',vname,tname,maxdvec,tol); 90 91 titlestr = sprintf('Diff between C and Octave of %s for %s',vname,tname) 92 figure(10+pnum*2) 93 plot(abs(dvec)) 94 title(titlestr) 95 96 figure(11+pnum*2) 97 plot((1:length(vc)),abs(vc),(1:length(voct)),abs(voct)) 98 99 end 100 101endfunction 102 103% Run C, then Octave verion of demod, and compare results 104function test_stats = fsk_demod_xt(Fs,Rs,f1,fsp,mod,tname,M=2,lock_nin=0) 105 global print_verbose; 106 global tfsk_location; 107 %Name of executable containing the modulator 108 fsk_demod_ex_file = '../build/unittest/tfsk'; 109 modvecfilename = sprintf('fsk_demod_ut_modvec_%d',getpid()); 110 bitvecfilename = sprintf('fsk_demod_ut_bitvec_%d',getpid()); 111 tvecfilename = sprintf('fsk_demod_ut_tracevec_%d.txt',getpid()); 112 113 %command to be run by system to launch the demod 114 command = sprintf('%s D %d %d %d %d %d %d %s %s %s',tfsk_location,M,f1,fsp,Fs,Rs,lock_nin,modvecfilename,bitvecfilename,tvecfilename); 115 116 %save modulated input into a file 117 modvecfile = fopen(modvecfilename,'wb+'); 118 fwrite(modvecfile,mod,'single'); 119 fclose(modvecfile); 120 121 %run the modulator 122 system(command); 123 124 bitvecfile = fopen(bitvecfilename,'rb'); 125 bits = fread(bitvecfile,'uint8'); 126 fclose(bitvecfile); 127 bits = bits!=0; 128 129 %Load test vec dump 130 load(tvecfilename) 131 132 %Clean up files 133 delete(bitvecfilename); 134 delete(modvecfilename); 135 delete(tvecfilename); 136 137 o_f1_dc = []; 138 o_f2_dc = []; 139 o_f3_dc = []; 140 o_f4_dc = []; 141 o_f1_int = []; 142 o_f2_int = []; 143 o_f3_int = []; 144 o_f4_int = []; 145 o_f1 = []; 146 o_f2 = []; 147 o_f3 = []; 148 o_f4 = []; 149 o_EbNodB = []; 150 o_ppm = []; 151 o_Sf = []; 152 o_fest = []; o_mask = []; o_fest2 = []; 153 o_rx_timing = []; 154 o_norm_rx_timing = []; 155 o_nin = []; 156 %Run octave demod, dump some test vectors 157 states = fsk_horus_init(Fs,Rs,M); 158 159 Ts = states.Ts; 160 P = states.P; 161 states.ftx(1) = f1; 162 states.ftx(2) = f1+fsp; 163 states.ftx(3) = f1+fsp*2; 164 states.ftx(4) = f1+fsp*3; 165 states.tx_tone_separation = fsp; 166 states.dF = 0; 167 modin = mod; 168 obits = []; 169 while length(modin)>=states.nin 170 ninold = states.nin; 171 states = est_freq(states, modin(1:states.nin), states.M); 172 [bitbuf,states] = fsk_demod(states, modin(1:states.nin)); 173 if lock_nin states.nin = states.N; end 174 175 modin=modin(ninold+1:length(modin)); 176 obits = [obits bitbuf]; 177 178 %Save other parameters 179 o_f1_dc = [o_f1_dc states.f_dc(1,:)]; 180 o_f2_dc = [o_f2_dc states.f_dc(2,:)]; 181 o_f1_int = [o_f1_int states.f_int(1,:)]; 182 o_f2_int = [o_f2_int states.f_int(2,:)]; 183 o_EbNodB = [o_EbNodB states.EbNodB]; 184 o_ppm = [o_ppm states.ppm]; 185 o_rx_timing = [o_rx_timing states.rx_timing]; 186 o_norm_rx_timing = [o_norm_rx_timing states.norm_rx_timing]; 187 o_Sf = [o_Sf states.Sf']; 188 o_f1 = [o_f1 states.f(1)]; 189 o_f2 = [o_f1 states.f(2)]; 190 o_fest = [o_fest states.f]; 191 o_mask = [o_mask states.mask]; 192 o_fest2 = [o_fest2 states.f2]; 193 o_nin = [o_nin states.nin]; 194 if M==4 195 o_f3_dc = [o_f3_dc states.f_dc(3,:)]; 196 o_f4_dc = [o_f4_dc states.f_dc(4,:)]; 197 o_f3_int = [o_f3_int states.f_int(3,:)]; 198 o_f4_int = [o_f4_int states.f_int(4,:)]; 199 o_f3 = [o_f1 states.f(3)]; 200 o_f4 = [o_f1 states.f(4)]; 201 end 202 end 203 204 assert(vcompare(o_Sf, t_Sf,'fft est',tname,.001,1)); 205 assert(vcompare(o_fest, t_f_est,'f est',tname,.001,2)); 206 assert(vcompare(o_mask, t_mask,'f2 mask',tname,.001,3)); 207 assert(vcompare(o_fest2, t_f2_est,'f2 est',tname,.001,16)); 208 o_fest2(1:12) 209 t_f2_est(1:12) 210 assert(vcompare(o_f1_dc, t_f1_dc, 'f1 dc', tname,.01,8)); 211 assert(vcompare(o_f2_dc, t_f2_dc, 'f2 dc', tname,.01,9)); 212 assert(vcompare(o_f2_int, t_f2_int, 'f2 int', tname,.01,10)); 213 assert(vcompare(o_f1_int, t_f1_int, 'f1 int', tname,.01,11)); 214 if M==4 215 assert(vcompare(o_f3_dc, t_f3_dc, 'f3 dc', tname,.01,4)) 216 assert(vcompare(o_f4_dc, t_f4_dc, 'f4 dc', tname,.01,5)); 217 assert(vcompare(o_f3_int, t_f3_int, 'f3 int', tname,.01,6)); 218 assert(vcompare(o_f4_int, t_f4_int, 'f4 int', tname,.01,7)); 219 end 220 221 assert(vcompare(o_rx_timing, t_rx_timing,'rx timing',tname,.02,3)); 222 223 % Much larger tolerances on unimportant statistics 224 assert(vcompare(o_ppm , t_ppm, 'ppm', tname,.02,12)); 225 assert(vcompare(o_EbNodB, t_EbNodB,'EbNodB', tname,.02,13)); 226 assert(vcompare(o_nin, t_nin, 'nin', tname,.0001,14)); 227 assert(vcompare(o_norm_rx_timing, t_norm_rx_timing,'norm rx timing',tname,.02,15)); 228 diffpass = sum(xor(obits,bits'))<4; 229 diffbits = sum(xor(obits,bits')); 230 231 if diffpass==0 232 printf('\n***bitcompare test failed test %s diff %d\n\n',tname,sum(xor(obits,bits'))) 233 figure(15) 234 plot(xor(obits,bits')) 235 title(sprintf('Bitcompare failure test %s',tname)) 236 end 237 238 assert(diffpass); 239 240 test_stats.pass = 1; 241 test_stats.diff = sum(xor(obits,bits')); 242 test_stats.cbits = bits'; 243 test_stats.obits = obits; 244endfunction 245 246 247function [dmod,cmod,omod,pass] = fsk_mod_test(Fs,Rs,f1,fsp,bits,tname,M=2) 248 global mod_pass_fail_maxdiff; 249 %Run the C modulator 250 cmod = fsk_mod_c(Fs,Rs,f1,fsp,bits,M); 251 %Set up and run the octave modulator 252 states.M = M; 253 states = fsk_horus_init(Fs,Rs,M); 254 255 states.ftx(1) = f1; 256 states.ftx(2) = f1+fsp; 257 258 if states.M == 4 259 states.ftx(3) = f1+fsp*2; 260 states.ftx(4) = f1+fsp*3; 261 end 262 263 states.dF = 0; 264 omod = fsk_mod(states,bits); 265 266 dmod = cmod-omod; 267 pass = max(dmod)<(mod_pass_fail_maxdiff*length(dmod)); 268 if !pass 269 printf('Mod failed test %s!\n',tname); 270 end 271endfunction 272 273% Random bit modulator test 274% Pass random bits through the modulators and compare 275function pass = test_mod_horuscfg_randbits 276 rand('state',1); 277 randn('state',1); 278 bits = rand(1,10000)>.5; 279 [dmod,cmod,omod,pass] = fsk_mod_test(8000,100,1200,1600,bits,"mod horuscfg randbits"); 280 281 if(!pass) 282 figure(1) 283 plot(dmod) 284 title("Difference between octave and C mod impl"); 285 end 286 print_result("test_mod_horuscfg_randbits", "OK"); 287endfunction 288 289% Random bit modulator test 290% Pass random bits through the modulators and compare 291function pass = test_mod_horuscfgm4_randbits 292 rand('state',1); 293 randn('state',1); 294 bits = rand(1,10000)>.5; 295 [dmod,cmod,omod,pass] = fsk_mod_test(8000,100,1200,1600,bits,"mod horuscfg randbits",4); 296 297 if(!pass) 298 figure(1) 299 plot(dmod) 300 title("Difference between octave and C mod impl"); 301 end 302 print_result("test_mod_horuscfgm4_randbits", "OK"); 303 304endfunction 305 306 307% A big ol' channel impairment tester shamelessly taken from fsk_horus 308% This throws some channel imparment or another at the C and octave 309% modem so they may be compared. 310function stats = tfsk_run_sim(test_frame_mode,EbNodB,timing_offset,fading,df,M=2,frames=50,lock_nin=0) 311 #{ 312 timing_offset [0|1] enable a 1000ppm sample clock offset 313 fading [0|1] modulates tx power at 2Hz with 20dB fade depth, 314 e.g. to simulate balloon rotating at end of mission 315 df tx tone freq drift in Hz/s 316 lock_nin [0|1] locks nin to a constant which makes tests much simpler by breaking feeback loop 317 #} 318 global print_verbose; 319 320 more off 321 rand('state',1); 322 randn('state',1); 323 324 % ---------------------------------------------------------------------- 325 326 % sm2000 config ------------------------ 327 %states = fsk_horus_init(96000, 1200); 328 %states.f1_tx = 4000; 329 %states.f2_tx = 5200; 330 331 if test_frame_mode == 2 332 % horus rtty config --------------------- 333 states = fsk_horus_init(8000, 100, M); 334 states.f1_tx = 1200; 335 states.f2_tx = 1600; 336 end 337 338 if test_frame_mode == 4 339 % horus rtty config --------------------- 340 states = fsk_horus_init(8000, 100, M); 341 states.f1_tx = 1200; 342 states.f2_tx = 1600; 343 states.tx_bits_file = "horus_tx_bits_rtty.txt"; % Octave file of bits we FSK modulate 344 345 end 346 347 if test_frame_mode == 5 348 % horus binary config --------------------- 349 states = fsk_horus_init(8000, 100, M); 350 states.f1_tx = 1200; 351 states.f2_tx = 1600; 352 %%%states.tx_bits_file = "horus_tx_bits_binary.txt"; % Octave file of bits we FSK modulate 353 states.tx_bits_file = "horus_payload_rtty.txt"; 354 end 355 356 % ---------------------------------------------------------------------- 357 358 states.verbose = 0; 359 N = states.N; 360 P = states.P; 361 Rs = states.Rs; 362 nsym = states.nsym; 363 Fs = states.Fs; 364 states.df = df; 365 states.M = M; 366 367 EbNo = 10^(EbNodB/10); 368 variance = states.Fs/(states.Rs*EbNo*states.bitspersymbol); 369 370 % set up tx signal with payload bits based on test mode 371 372 if test_frame_mode == 1 373 % test frame of bits, which we repeat for convenience when BER testing 374 test_frame = round(rand(1, states.nsym)); 375 tx_bits = []; 376 for i=1:frames+1 377 tx_bits = [tx_bits test_frame]; 378 end 379 end 380 if test_frame_mode == 2 381 % random bits, just to make sure sync algs work on random data 382 tx_bits = round(rand(1, states.nbit*(frames+1))); 383 end 384 if test_frame_mode == 3 385 % ...10101... sequence 386 tx_bits = zeros(1, states.nsym*(frames+1)); 387 tx_bits(1:2:length(tx_bits)) = 1; 388 end 389 390 if (test_frame_mode == 4) || (test_frame_mode == 5) 391 392 % load up a horus msg from disk and modulate that 393 394 test_frame = load(states.tx_bits_file); 395 ltf = length(test_frame); 396 ntest_frames = ceil((frames+1)*nsym/ltf); 397 tx_bits = []; 398 for i=1:ntest_frames 399 tx_bits = [tx_bits test_frame]; 400 end 401 end 402 403 f1 = states.f1_tx; 404 fsp = states.f2_tx-f1; 405 states.ftx(1) = f1; 406 states.ftx(2) = f1+fsp; 407 408 if states.M == 4 409 states.ftx(3) = f1+fsp*2; 410 states.ftx(4) = f1+fsp*3; 411 end 412 413 tx = fsk_mod(states, tx_bits); 414 415 if timing_offset 416 tx = resample(tx, 1000, 1001); % simulated 1000ppm sample clock offset 417 end 418 419 if fading 420 ltx = length(tx); 421 tx = tx .* (1.1 + cos(2*pi*2*(0:ltx-1)/Fs))'; % min amplitude 0.1, -20dB fade, max 3dB 422 end 423 424 noise = sqrt(variance)*randn(length(tx),1); 425 rx = tx + noise; 426 427 test_name = sprintf("tfsk run sim EbNodB:%d frames:%d timing_offset:%d fading:%d df:%d",EbNodB,frames,timing_offset,fading,df); 428 tstats = fsk_demod_xt(Fs,Rs,states.f1_tx,fsp,rx,test_name,M,lock_nin); 429 430 pass = tstats.pass; 431 obits = tstats.obits; 432 cbits = tstats.cbits; 433 434 % Figure out BER of octave and C modems 435 bitcnt = length(tx_bits); 436 rx_bits = obits; 437 ber = 1; 438 ox = 1; 439 for offset = (1:100) 440 nerr = sum(xor(rx_bits(offset:length(rx_bits)),tx_bits(1:length(rx_bits)+1-offset))); 441 bern = nerr/(bitcnt-offset); 442 if(bern < ber) 443 ox = offset; 444 best_nerr = nerr; 445 end 446 ber = min([ber bern]); 447 end 448 offset = ox; 449 bero = ber; 450 ber = 1; 451 rx_bits = cbits; 452 ox = 1; 453 for offset = (1:100) 454 nerr = sum(xor(rx_bits(offset:length(rx_bits)),tx_bits(1:length(rx_bits)+1-offset))); 455 bern = nerr/(bitcnt-offset); 456 if(bern < ber) 457 ox = offset; 458 best_nerr = nerr; 459 end 460 ber = min([ber bern]); 461 end 462 offset = ox; 463 berc = ber; 464 stats.berc = berc; 465 stats.bero = bero; 466 stats.name = test_name; 467 % coherent BER theory calculation 468 469 stats.thrcoh = .5*(M-1)*erfc(sqrt( (log2(M)/2) * EbNo )); 470 471 % non-coherent BER theory calculation 472 % It was complicated, so I broke it up 473 474 ms = M; 475 ns = (1:ms-1); 476 as = (-1).^(ns+1); 477 bs = (as./(ns+1)); 478 479 cs = ((ms-1)./ns); 480 481 ds = ns.*log2(ms); 482 es = ns+1; 483 fs = exp( -(ds./es)*EbNo ); 484 485 thrncoh = ((ms/2)/(ms-1)) * sum(bs.*((ms-1)./ns).*exp( -(ds./es)*EbNo )); 486 487 stats.thrncoh = thrncoh; 488 stats.pass = pass; 489endfunction 490 491% run a bunch of tests at a range of EbNo's in parallel 492function pass = ebno_battery_test(timing_offset,fading,df,M) 493 %Range of EbNodB over which to test 494 ebnodbrange = (5:2:13); 495 ebnodbs = length(ebnodbrange); 496 497 mode = 2; 498 %Replication of other parameters for parcellfun 499 modev = repmat(mode,1,ebnodbs); 500 timingv = repmat(timing_offset,1,ebnodbs); 501 fadingv = repmat(fading,1,ebnodbs); 502 dfv = repmat(df,1,ebnodbs); 503 mv = repmat(M,1,ebnodbs); 504 505 statv = pararrayfun(floor(1.25*nproc()),@tfsk_run_sim,modev,ebnodbrange,timingv,fadingv,dfv,mv); 506 passv = zeros(1,length(statv)); 507 for ii=(1:length(statv)) 508 passv(ii)=statv(ii).pass; 509 if statv(ii).pass 510 printf("Test %s passed\n",statv(ii).name); 511 else 512 printf("Test %s failed\n",statv(ii).name); 513 end 514 end 515 516 %All pass flags are '1' 517 pass = sum(passv)>=length(passv); 518 %and no tests died 519 pass = pass && length(passv)==ebnodbs; 520 passv; 521 assert(pass) 522endfunction 523 524%Test with and without sample clock offset 525function pass = test_timing_var(df,M) 526 pass = ebno_battery_test(1,0,df,M) 527 assert(pass) 528 pass = pass && ebno_battery_test(0,0,df,M) 529 assert(pass) 530endfunction 531 532%Test with and without 1 Hz/S freq drift 533function pass = test_drift_var(M) 534 pass = test_timing_var(1,1,M) 535 pass = pass && test_timing_var(0,1,M) 536 assert(pass) 537endfunction 538 539function pass = test_fsk_battery() 540 pass = test_mod_horuscfg_randbits; 541 pass = pass && test_mod_horuscfgm4_randbits; 542 pass = pass && test_drift_var(4); 543 assert(pass) 544 pass = pass && test_drift_var(2); 545 assert(pass) 546 if pass 547 printf("***** All tests passed! *****\n"); 548 end 549endfunction 550 551function plot_fsk_bers(M=2) 552 %Range of EbNodB over which to plot 553 ebnodbrange = (4:13); 554 555 berc = ones(1,length(ebnodbrange)); 556 bero = ones(1,length(ebnodbrange)); 557 berinc = ones(1,length(ebnodbrange)); 558 beric = ones(1,length(ebnodbrange)); 559 ebnodbs = length(ebnodbrange) 560 mode = 2; 561 %Replication of other parameters for parcellfun 562 modev = repmat(mode,1,ebnodbs); 563 timingv = repmat(1,1,ebnodbs); 564 fadingv = repmat(0,1,ebnodbs); 565 dfv = repmat(1,1,ebnodbs); 566 Mv = repmat(M,1,ebnodbs); 567 568 statv = pararrayfun(floor(nproc()),@tfsk_run_sim,modev,ebnodbrange,timingv,fadingv,dfv,Mv); 569 %statv = arrayfun(@tfsk_run_sim,modev,ebnodbrange,timingv,fadingv,dfv,Mv); 570 571 for ii = (1:length(statv)) 572 stat = statv(ii); 573 berc(ii)=stat.berc; 574 bero(ii)=stat.bero; 575 berinc(ii)=stat.thrncoh; 576 beric(ii) = stat.thrcoh; 577 end 578 clf; 579 figure(M) 580 581 semilogy(ebnodbrange, berinc,sprintf('r;%dFSK non-coherent theory;',M)) 582 hold on; 583 semilogy(ebnodbrange, beric ,sprintf('g;%dFSK coherent theory;',M)) 584 semilogy(ebnodbrange, bero ,sprintf('b;Octave fsk horus %dFSK Demod;',M)) 585 semilogy(ebnodbrange, berc,sprintf('+;C fsk horus %dFSK Demod;',M)) 586 hold off; 587 grid("minor"); 588 axis([min(ebnodbrange) max(ebnodbrange) 1E-5 1]) 589 legend("boxoff"); 590 xlabel("Eb/No (dB)"); 591 ylabel("Bit Error Rate (BER)") 592 593endfunction 594 595% We kick off tests here ------------------------------------------------------ 596 597pass = 0; ntests = 0; 598pass += test_mod_horuscfg_randbits; ntests++; 599pass += test_mod_horuscfgm4_randbits; ntests++; 600stats = tfsk_run_sim(test_frame_mode=2,EbNodB=5,timing_offset=0,fading=0,df=1,M=4,frames=10,lock_nin=1); ntests++; 601if stats.pass 602 print_result("Demod 10 frames nin locked", "OK"); 603 pass += stats.pass; 604end 605stats = tfsk_run_sim(test_frame_mode=2,EbNodB=5,timing_offset=1,fading=0,df=1,M=4,frames=10,lock_nin=0); ntests++; 606if stats.pass 607 print_result("Demod 10 frames", "OK"); 608 pass += stats.pass; 609end 610printf("tests: %d passed: %d ", ntests, pass); 611if ntests == pass printf("PASS\n"); else printf("FAIL\n"); end 612