1% test_ldpc_fsk_lib 2% David Rowe 16 April 2016 3% 4% A series of tests for ldpc_fsk_lib, and C versions ldpc_enc and ldpc_dec. 5% Gradually builds up complete C command line for SSTV balloon system, 6% using Octave versions of LDPC and FSK modem as reference points. 7 81; 9 10% encodes and decodes one frame, also writes codeword.bin for testing 11% decode_from_file() below, and can optionally generate include file for 12% C version of decoder. 13 14function [data code_param] = simple_ut(c_include_file) 15 load('H2064_516_sparse.mat'); 16 HRA = full(HRA); 17 max_iterations = 100; 18 decoder_type = 0; 19 EsNodB = 3; 20 mod_order = 2; 21 22 code_param = ldpc_init(HRA, mod_order); 23 data = round( rand( 1, code_param.data_bits_per_frame ) ); 24 codeword = ldpc_encode(code_param, data); 25 f = fopen("codeword.bin","wt"); fwrite(f, codeword, "uint8"); fclose(f); 26 s = 1 - 2 * codeword; 27 code_param.symbols_per_frame = length( s ); 28 29 EsNo = 10^(EsNodB/10); 30 variance = 1/(2*EsNo); 31 noise = sqrt(variance)* randn(1,code_param.symbols_per_frame); 32 rx = s + noise; 33 34 if nargin == 1 35 code_param.c_include_file = c_include_file; 36 end 37 [detected_data Niters] = ldpc_decode(rx, code_param, max_iterations, decoder_type); 38 39 error_positions = xor(detected_data(1:code_param.data_bits_per_frame), data); 40 Nerrs = sum(error_positions); 41 42 printf("Nerrs = %d\n", Nerrs); 43end 44 45 46% This version decodes from a file of bits 47 48function detected_data = decode_from_file(filename) 49 max_iterations = 100; 50 decoder_type = 0; 51 load('H2064_516_sparse.mat'); 52 HRA = full(HRA); 53 mod_order = 2; 54 55 f = fopen(filename,"rb"); codeword = fread(f, "uint8")'; fclose(f); 56 r = 1 - 2 * codeword; 57 code_param = ldpc_init(HRA, mod_order); 58 [detected_data Niters] = ldpc_decode(r, code_param, max_iterations, decoder_type); 59end 60 61 62% plots a BER curve for the LDPC decoder. Takes a while to run, uses parallel cores 63 64function plot_curve 65 num_cores = 4; % set this to the number of cores you have 66 67 load('H2064_516_sparse.mat'); 68 HRA = full(HRA); 69 [Nr Nc] = size(HRA); 70 sim_in.rate = (Nc-Nr)/Nc; 71 72 sim_in.HRA = HRA; 73 sim_in.mod_order = 2; 74 sim_in.framesize = Nc; 75 sim_in.mod_order = 2; 76 sim_in.Lim_Ferrs = 100; 77 78 % note we increase number of trials as BER goes down 79 80 Esvec = [ 0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 ]; 81 Ntrials = [ 1E4 1E4 1E4 1E4 1E5 1E5 1E5 1E5 1E5 ]; 82 num_runs = length(Esvec) 83 84 sim_in_vec(1:num_runs) = sim_in; 85 for i = 1:num_runs 86 sim_in_vec(i).Esvec = Esvec(i); 87 sim_in_vec(i).Ntrials = Ntrials(i); 88 end 89 90 %sim_out = ldpc5(sim_in_vec(1)); 91 tstart = time(); 92 sim_out = pararrayfun(num_cores, @ldpc5, sim_in_vec); 93 tend = time(); 94 95 total_bits = sum(Ntrials)*sim_in.framesize; 96 total_secs = tend - tstart; 97 printf("%d bits in %4.1f secs, or %5f bits/s\n", total_bits, total_secs, total_bits/total_secs); 98 99 for i=1:num_runs 100 Ebvec(i) = sim_out(i).Ebvec; 101 BERvec(i) = sim_out(i).BERvec; 102 end 103 semilogy(Ebvec, BERvec, '+-') 104 xlabel('Eb/N0') 105 ylabel('BER') 106 title(['H2064 516 sparse.mat' ' ' date]) 107 108end 109 110 111% Test C encoder 112 113function test_c_encoder 114 load('H2064_516_sparse.mat'); 115 HRA = full(HRA); 116 max_iterations = 100; 117 decoder_type = 0; 118 EsNodB = 3; 119 mod_order = 2; 120 frames = 100; 121 122 EsNo = 10^(EsNodB/10); 123 variance = 1/(2*EsNo); 124 125 code_param = ldpc_init(HRA, mod_order); 126 127 data = round(rand(1,frames*code_param.data_bits_per_frame)); 128 f = fopen("data.bin","wt"); fwrite(f, data, "uint8"); fclose(f); 129 130 % Outboard C encoder 131 132 system("../src/ldpc_enc data.bin codewords.bin"); 133 134 % Test with Octave decoder 135 136 f = fopen("codewords.bin","rb"); codewords = fread(f, "uint8")'; fclose(f); 137 138 Nerrs = 0; 139 for i=1:frames 140 st = (i-1)*code_param.symbols_per_frame+1; en = st+code_param.symbols_per_frame-1; 141 tx = 1 - 2 * codewords(st:en); 142 143 noise = sqrt(variance)*randn(1,code_param.symbols_per_frame); 144 rx = tx + noise; 145 146 [detected_data Niters] = ldpc_decode(rx, code_param, max_iterations, decoder_type); 147 148 st = (i-1)*code_param.data_bits_per_frame+1; en = st+code_param.data_bits_per_frame-1; 149 error_positions = xor(detected_data(1:code_param.data_bits_per_frame), data(st:en)); 150 Nerrs += sum(error_positions); 151 end 152 153 printf("Nerrs = %d\n", Nerrs); 154end 155 156 157function test_c_decoder 158 load('H2064_516_sparse.mat'); 159 HRA = full(HRA); 160 max_iterations = 100; 161 decoder_type = 0; 162 mod_order = 2; 163 frames = 10; 164 EsNodB = 2; 165 sdinput = 1; 166 167 EsNo = 10^(EsNodB/10); 168 variance = 1/(2*EsNo); 169 170 code_param = ldpc_init(HRA, mod_order); 171 data = round(rand(1,code_param.data_bits_per_frame*frames)); 172 173 f = fopen("data.bin","wt"); fwrite(f, data, "uint8"); fclose(f); 174 system("../src/ldpc_enc data.bin codewords.bin"); 175 f = fopen("codewords.bin","rb"); codewords = fread(f, "uint8")'; fclose(f); 176 177 s = 1 - 2 * codewords; 178 noise = sqrt(variance)*randn(1,code_param.symbols_per_frame*frames); 179 r = s + noise; 180 181 % calc LLRs frame by frame 182 183 for i=1:frames 184 st = (i-1)*code_param.symbols_per_frame+1; 185 en = st + code_param.symbols_per_frame-1; 186 llr(st:en) = sd_to_llr(r(st:en)); 187 end 188 189 % Outboard C decoder 190 191 if sdinput 192 f = fopen("sd.bin","wb"); fwrite(f, r, "double"); fclose(f); 193 system("../src/ldpc_dec sd.bin data_out.bin --sd"); 194 else 195 f = fopen("llr.bin","wb"); fwrite(f, llr, "double"); fclose(f); 196 system("../src/ldpc_dec llr.bin data_out.bin"); 197 end 198 199 f = fopen("data_out.bin","rb"); data_out = fread(f, "uint8")'; fclose(f); 200 201 Nerrs = Nerrs2 = zeros(1,frames); 202 for i=1:frames 203 204 % Check C decoder 205 206 data_st = (i-1)*code_param.data_bits_per_frame+1; 207 data_en = data_st+code_param.data_bits_per_frame-1; 208 st = (i-1)*code_param.symbols_per_frame+1; 209 en = st+code_param.data_bits_per_frame-1; 210 data_out_c = data_out(st:en); 211 error_positions = xor(data_out_c, data(data_st:data_en)); 212 Nerrs(i) = sum(error_positions); 213 214 % Octave decoder 215 216 st = (i-1)*code_param.symbols_per_frame+1; en = st+code_param.symbols_per_frame-1; 217 [detected_data Niters] = ldpc_decode(r(st:en), code_param, max_iterations, decoder_type); 218 st = (i-1)*code_param.data_bits_per_frame+1; en = st+code_param.data_bits_per_frame-1; 219 data_out_octave = detected_data(1:code_param.data_bits_per_frame); 220 error_positions = xor(data_out_octave, data(st:en)); 221 Nerrs2(i) = sum(error_positions); 222 %printf("%4d ", Niters); 223 end 224 printf("Errors per frame:\nC.....:"); 225 for i=1:frames 226 printf("%4d ", Nerrs(i)); 227 end 228 printf("\nOctave:"); 229 for i=1:frames 230 printf("%4d ", Nerrs2(i)); 231 end 232 printf("\n"); 233 234end 235 236% Saves a complex vector s to a file "filename" of IQ unsigned 8 bit 237% chars, same as RTLSDR format. 238 239function save_rtlsdr(filename, s) 240 mx = max(abs(s)); 241 re = real(s); im = imag(s); 242 l = length(s); 243 iq = zeros(1,2*l); 244 %iq(1:2:2*l) = 127 + re*(127/mx); 245 %iq(2:2:2*l) = 127 + im*(127/mx); 246 iq(1:2:2*l) = 127 + 32*re; 247 iq(2:2:2*l) = 127 + 32*im; 248 figure(3); clf; plot(iq); title('simulated IQ signal from RTL SDR'); 249 fs = fopen(filename,"wb"); 250 fwrite(fs,iq,"uchar"); 251 fclose(fs); 252endfunction 253 254 255% Oversamples by a factor of 2 using Octaves resample() function then 256% uses linear interpolation to achive fractional sample rate 257 258function rx_resample_fract = fractional_resample(rx, resample_rate); 259 assert(resample_rate < 2, "keep resample_rate between 0 and 2"); 260 rx_resample2 = resample(rx, 2, 1); 261 l = length(rx_resample2); 262 rx_resample_fract = zeros(1,l); 263 k = 1; 264 step = 2/resample_rate; 265 for i=1:step:l-1 266 i_low = floor(i); 267 i_high = ceil(i); 268 f = i - i_low; 269 rx_resample_fract(k) = (1-f)*rx_resample2(i_low) + f*rx_resample2(i_high); 270 %printf("i: %f i_low: %d i_high: %d f: %f\n", i, i_low, i_high, f); 271 k++; 272 end 273 rx_resample_fract = rx_resample_fract(1:k-1); 274endfunction 275 276 277% Using simulated SSTV packet, generate complex fsk mod signals, 8-bit 278% unsigned IQ for feeding into C demod chain. Can also be used to 279% generate BER curves. Found bugs in UW size and our use of csdr 280% re-sampler using this function, and by gradually and carefully 281% building up the C command line. 282 283#{ 284todo: [X] uncoded BER 285 [X] octave fsk demod 286 [X] use C demod 287 [X] compare uncoded BER to unsigned 8 bit IQ to regular 16-bit 288 [X] generate complex rx signal with noise 289 [X] used cmd line utils to drive demod 290 [X] test with resampler 291 [X] measure effect on PER with coding 292#} 293 294function [n_uncoded_errs n_uncoded_bits] = run_sstv_sim(sim_in, EbNodB) 295 296 frames = sim_in.frames; 297 demod_type = sim_in.demod_type; 298 299 % init LDPC code 300 301 load('H2064_516_sparse.mat'); 302 HRA = full(HRA); 303 max_iterations = 100; 304 decoder_type = 0; 305 mod_order = 2; 306 307 code_param = ldpc_init(HRA, mod_order); 308 309 % note fixed frame of bits used for BER testing 310 311 tx_codeword = gen_sstv_frame; 312 313 % init FSK modem 314 315 fsk_horus_as_a_lib = 1; 316 fsk_horus; 317 states = fsk_horus_init_hbr(9600, 8, 1200, 2, length(tx_codeword)); 318 states.df(1:states.M) = 0; 319 states.dA(1:states.M) = 1; 320 states.tx_real = 0; % Octave fsk_mod generates complex valued output 321 % so we can simulate rtl_sdr complex ouput 322 323 % Set up simulated tx tones to sit in the middle of cdsr passband 324 325 filt_low_norm = 0.1; filt_high_norm = 0.4; 326 fc = states.Fs*(filt_low_norm + filt_high_norm)/2; 327 %fc = 1800; 328 f1 = fc - states.Rs/2; 329 f2 = fc + states.Rs/2; 330 states.ftx = [f1 f2]; 331 332 % set up AWGN channel 333 334 EbNo = 10^(EbNodB/10); 335 variance = states.Fs/(states.Rs*EbNo*states.bitspersymbol); 336 337 % start simulation ---------------------------------------- 338 339 tx_bit_stream = []; 340 for i=1:frames 341 % uncomment for different data on each frame 342 %tx_codeword = gen_sstv_frame; 343 tx_bit_stream = [tx_bit_stream tx_codeword]; 344 end 345 346 printf("%d bits at %d bit/s is a %3.1f second run\n", length(tx_bit_stream), 115200,length(tx_bit_stream)/115200); 347 348 % modulate and channel model 349 350 tx = fsk_horus_mod(states, tx_bit_stream); 351 noise_real = sqrt(variance)*randn(length(tx),1); 352 noise_complex = sqrt(variance/2)*(randn(length(tx),1) + j*randn(length(tx),1)); 353 354 % demodulate ----------------------------------------------------- 355 356 if demod_type == 1 357 358 % Octave demod 359 360 if states.tx_real 361 rx = tx + noise_real; 362 else 363 rx = tx + noise_complex; 364 end 365 SNRdB = 10*log10(var(tx)/var(noise_complex)); 366 367 % demodulate frame by frame using Octave demod 368 369 st = 1; 370 run_frames = floor(length(rx)/states.N); 371 rx_bit_stream = []; 372 rx_sd_stream = []; 373 for f=1:run_frames 374 375 % extract nin samples from rx sample stream 376 377 nin = states.nin; 378 en = st + states.nin - 1; 379 380 if en <= length(rx) % due to nin variations its possible to overrun buffer 381 sf = rx(st:en); 382 st += nin; 383 384 % demodulate to stream of bits 385 386 states.f = [f1 f2]; % note that for Octave demod we cheat and use known tone frequencies 387 % allows us to determine if freq offset estimation in C demod is a problem 388 389 [rx_bits states] = fsk_horus_demod(states, sf); 390 rx_bit_stream = [rx_bit_stream rx_bits]; 391 rx_sd_stream = [rx_sd_stream states.rx_bits_sd]; 392 end 393 end 394 end 395 396 if demod_type == 2 397 % baseline C demod 398 399 if states.tx_real 400 rx = tx + noise_real; 401 else 402 rx = 2*real(tx) + noise_real; 403 end 404 SNRdB = 10*log10(var(tx)/var(noise_real)); 405 rx_scaled = 1000*real(rx); 406 f = fopen("fsk_demod.raw","wb"); fwrite(f, rx_scaled, "short"); fclose(f); 407 system("../build_linux/src/fsk_demod 2X 8 9600 1200 fsk_demod.raw fsk_demod.bin"); 408 f = fopen("fsk_demod.bin","rb"); rx_bit_stream = fread(f, "uint8")'; fclose(f); 409 end 410 411 if demod_type == 3 412 % C demod driven by csdr command line kung fu 413 414 assert(states.tx_real == 0, "need complex signal for this test"); 415 rx = tx + noise_complex; 416 SNRdB = 10*log10(var(tx)/var(noise_real)); 417 save_rtlsdr("fsk_demod.iq", rx); 418 system("cat fsk_demod.iq | csdr convert_u8_f | csdr bandpass_fir_fft_cc 0.1 0.4 0.05 | csdr realpart_cf | csdr convert_f_s16 | ../build_linux/src/fsk_demod 2X 8 9600 1200 - fsk_demod.bin"); 419 f = fopen("fsk_demod.bin","rb"); rx_bit_stream = fread(f, "uint8")'; fclose(f); 420 end 421 422 if demod_type == 4 423 % C demod with resampler ....... getting closer to Mark's real time cmd line 424 425 assert(states.tx_real == 0, "need complex signal for this test"); 426 rx = tx + noise_complex; 427 SNRdB = 10*log10(var(tx)/var(noise_real)); 428 429 printf("resampling ...\n"); 430 rx_resample_fract = fractional_resample(rx, 1.08331); 431 %rx_resample_fract = fractional_resample(rx_resample_fract, 1/1.08331); 432 save_rtlsdr("fsk_demod_resample.iq", rx_resample_fract); 433 434 printf("run C cmd line chain ...\n"); 435% system("cat fsk_demod_resample.iq | csdr convert_u8_f | csdr bandpass_fir_fft_cc 0.1 0.4 0.05 | csdr realpart_cf | csdr convert_f_s16 | ../build_linux/src/fsk_demod 2X 8 9600 1200 - fsk_demod.bin"); 436 system("cat fsk_demod_resample.iq | csdr convert_u8_f | csdr bandpass_fir_fft_cc 0.1 0.4 0.05 | csdr realpart_cf | csdr convert_f_s16 | ../unittest/tsrc - - 0.9230968 | ../build_linux/src/fsk_demod 2X 8 9600 1200 - fsk_demod.bin"); 437% system("cat fsk_demod_resample.iq | csdr convert_u8_f | csdr bandpass_fir_fft_cc 0.1 0.4 0.05 | csdr realpart_cf | csdr fractional_decimator_ff 1.08331 | csdr convert_f_s16 | ../build_linux/src/fsk_demod 2X 8 9600 1200 - fsk_demod.bin"); 438 f = fopen("fsk_demod.bin","rb"); rx_bit_stream = fread(f, "uint8")'; fclose(f); 439 end 440 441 442 if demod_type == 5 443 444 % C demod with resampler and use C code to measure PER, in this 445 % test we don't need to run state machine below as C code gives us 446 % the ouputs we need 447 448 assert(states.tx_real == 0, "need complex signal for this test"); 449 rx = tx + noise_complex; 450 SNRdB = 10*log10(var(tx)/var(noise_real)); 451 452 printf("fract resampling ...\n"); 453 rx_resample_fract = fractional_resample(rx, 1.08331); 454 save_rtlsdr("fsk_demod_resample.iq", rx_resample_fract); 455 456 % useful for HackRF 457 %printf("10X resampling ...\n"); 458 %rx_resample_10M = resample(rx_resample_fract, 10, 1); 459 %save_rtlsdr("fsk_demod_10M.iq", rx_resample_10M); 460 461 printf("run C cmd line chain - uncoded PER\n"); 462 system("cat fsk_demod_resample.iq | csdr convert_u8_f | csdr bandpass_fir_fft_cc 0.1 0.4 0.05 | csdr realpart_cf | csdr convert_f_s16 | ../unittest/tsrc - - 0.9230968 | ../build_linux/src/fsk_demod 2X 8 9600 1200 - - | ../src/drs232 - /dev/null -v"); 463 464 printf("run C cmd line chain - LDPC coded PER\n"); 465 system("cat fsk_demod_resample.iq | csdr convert_u8_f | csdr bandpass_fir_fft_cc 0.1 0.4 0.05 | csdr realpart_cf | csdr convert_f_s16 | ../unittest/tsrc - - 0.9230968 | ../build_linux/src/fsk_demod 2XS 8 9600 1200 - - | ../src/drs232_ldpc - /dev/null -v"); 466 end 467 468 if demod_type == 6 469 % C demod with complex input driven simplfied csdr command line, just measure BER of demod 470 471 assert(states.tx_real == 0, "need complex signal for this test"); 472 rx = tx + noise_complex; 473 SNRdB = 10*log10(var(tx)/var(noise_real)); 474 save_rtlsdr("fsk_demod.iq", rx); 475 system("cat fsk_demod.iq | csdr convert_u8_f | csdr convert_f_s16 | ../build_linux/src/fsk_demod 2X 8 9600 1200 - fsk_demod.bin C"); 476 477 f = fopen("fsk_demod.bin","rb"); rx_bit_stream = fread(f, "uint8")'; fclose(f); 478 end 479 480 if demod_type == 7 481 % C demod with complex input, measure uncoded and uncoded PER 482 483 assert(states.tx_real == 0, "need complex signal for this test"); 484 rx = tx + noise_complex; 485 SNRdB = 10*log10(var(tx)/var(noise_real)); 486 save_rtlsdr("fsk_demod.iq", rx); 487 488 printf("run C cmd line chain - uncoded PER\n"); 489 system("cat fsk_demod.iq | csdr convert_u8_f | csdr convert_f_s16 | ../build_linux/src/fsk_demod 2X 8 9600 1200 - - C | ../src/drs232 - /dev/null -v"); 490 491 printf("run C cmd line chain - LDPC coded PER\n"); 492 %system("cat fsk_demod.iq | csdr convert_u8_f | csdr convert_f_s16 | ../build_linux/src/fsk_demod 2XS 8 9600 1200 - - C | ../src/drs232_ldpc - /dev/null -v"); 493 system("cat fsk_demod.iq | ../build_linux/src/fsk_demod 2XS 8 9600 1200 - - CU8 | ../src/drs232_ldpc - /dev/null -v"); 494 end 495 496 if (demod_type != 5) && (demod_type != 7) 497 % state machine. Look for SSTV UW. When found count bit errors over one frame of bits 498 499 state = "wait for uw"; 500 start_uw_ind = 16*10+1; end_uw_ind = start_uw_ind + 5*10 - 1; 501 uw_rs232 = tx_codeword(start_uw_ind:end_uw_ind); luw = length(uw_rs232); 502 start_frame_ind = end_uw_ind + 1; 503 nbits = length(rx_bit_stream); 504 uw_thresh = 5; 505 n_uncoded_errs = 0; 506 n_uncoded_bits = 0; 507 n_packets_rx = 0; 508 last_i = 0; 509 510 % might as well include RS232 framing bits in uncoded error count 511 512 nbits_frame = code_param.data_bits_per_frame*10/8; 513 514 uw_errs = zeros(1, nbits); 515 for i=luw:nbits 516 uw_errs(i) = sum(xor(rx_bit_stream(i-luw+1:i), uw_rs232)); 517 end 518 519 for i=luw:nbits 520 next_state = state; 521 if strcmp(state, 'wait for uw') 522 if uw_errs(i) <= uw_thresh 523 next_state = 'count errors'; 524 tx_frame_ind = start_frame_ind; 525 rx_frame_ind = i + 1; 526 n_uncoded_errs_this_frame = 0; 527 %printf("%d %s %s\n", i, state, next_state); 528 if last_i 529 printf("i: %d i-last_i: %d ", i, i-last_i); 530 end 531 end 532 end 533 if strcmp(state, 'count errors') 534 n_uncoded_errs_this_frame += xor(rx_bit_stream(i), tx_codeword(tx_frame_ind)); 535 n_uncoded_bits++; 536 tx_frame_ind++; 537 if tx_frame_ind == (start_frame_ind+nbits_frame) 538 n_uncoded_errs += n_uncoded_errs_this_frame; 539 printf("n_uncoded_errs_this_frame: %d\n", n_uncoded_errs_this_frame); 540 frame_rx232_rx = rx_bit_stream(rx_frame_ind:rx_frame_ind+nbits_frame-1); 541 %tx_codeword(start_frame_ind+1:start_frame_ind+10) 542 %frame_rx232_rx(1:10) 543 sstv_checksum(frame_rx232_rx); 544 last_i = i; 545 n_packets_rx++; 546 next_state = 'wait for uw'; 547 end 548 end 549 state = next_state; 550 end 551 552 uncoded_ber = n_uncoded_errs/n_uncoded_bits; 553 printf("EbNodB: %4.1f SNRdB: %4.1f pkts: %d bits: %d errs: %d BER: %4.3f\n", 554 EbNodB, SNRdB, n_packets_rx, n_uncoded_bits, n_uncoded_errs, uncoded_ber); 555 556 figure(2); 557 plot(uw_errs); 558 title('Unique Word Hamming Distance') 559 end 560 561endfunction 562 563 564% Function to test flight mode software. Takes a rx stream of 565% demodulated bits, and locates frames using UW detection. Extracts 566% data and parity bits. Uses data bits to generate parity bits here 567% and compare. 568 569function compare_parity_bits(rx_bit_stream) 570 nframes = 500; 571 572 % init LDPC code 573 574 load('H2064_516_sparse.mat'); 575 HRA = full(HRA); 576 max_iterations = 100; 577 decoder_type = 0; 578 mod_order = 2; 579 580 code_param = ldpc_init(HRA, mod_order); 581 582 % generate frame, this will have random bits not related to 583 % rx_stream, however we just use it for the UW 584 585 tx_codeword = gen_sstv_frame; 586 l = length(tx_codeword); 587 printf("expected rs232 frames codeword length: %d\n", l); 588 589 % state machine. Look for SSTV UW. When found count bit errors over one frame of bits 590 591 state = "wait for uw"; 592 start_uw_ind = 16*10+1; end_uw_ind = start_uw_ind + 4*10 - 1; 593 uw_rs232 = tx_codeword(start_uw_ind:end_uw_ind); luw = length(uw_rs232); 594 start_frame_ind = end_uw_ind + 1; 595 nbits = nframes*l; 596 uw_thresh = 5; 597 n_uncoded_errs = 0; 598 n_uncoded_bits = 0; 599 n_packets_rx = 0; 600 last_i = 0; 601 602 % might as well include RS232 framing bits in uncoded error count 603 604 uw_errs = luw*ones(1, nbits); 605 for i=luw:nbits 606 uw_errs(i) = sum(xor(rx_bit_stream(i-luw+1:i), uw_rs232)); 607 end 608 609 frame_start = find(uw_errs < 2)+1; 610 nframes = length(frame_start) 611 for i=1:nframes 612 613 % double check UW OK 614 615 st_uw = frame_start(i) - luw; en_uw = frame_start(i) - 1; 616 uw_err_check = sum(xor(rx_bit_stream(st_uw:en_uw), uw_rs232)); 617 %printf("uw_err_check: %d\n", uw_err_check); 618 619 % strip off rs232 start/stop bits 620 621 nbits_rs232 = (256+2+65)*10; 622 nbits = (256+2+65)*8; 623 nbits_byte = 10; 624 rx_codeword = zeros(1,nbits); 625 pdb = 1; 626 627 for k=1:nbits_byte:nbits_rs232 628 for l=1:8 629 rx_codeword(pdb) = rx_bit_stream(frame_start(i)-1+k+l); 630 pdb++; 631 end 632 end 633 assert(pdb == (nbits+1)); 634 635 data_bits = rx_codeword(1:256*8); 636 checksum_bits = rx_codeword(256*8+1:258*8); 637 parity_bits = rx_codeword(258*8+1:258*8+516); 638 padding_bits = rx_codeword(258*8+516+1:258*8+516+1); 639 640 % stopped here as we found bug lol! 641 end 642 643 figure(1); clf; 644 plot(uw_errs); 645 title('Unique Word Hamming Distance') 646 figure(2); clf; 647 lframe_start = length(frame_start); 648 plot(frame_start(2:lframe_start)-frame_start(1:lframe_start-1)); 649 %title('Unique Word Hamming Distance') 650 651endfunction 652 653 654% Start simulation -------------------------------------------------------- 655 656more off; 657currentdir = pwd; 658thiscomp = computer; 659 660if strcmpi(thiscomp, 'MACI64')==1 661 if exist('CMLSimulate')==0 662 cd '/Users/bill/Current/Projects/DLR_FSO/Visit2013_FSO_GEO/cml' 663 addpath '../' % assume the source files stored here 664 CmlStartup % note that this is not in the cml path! 665 disp('added MACI64 path and run CmlStartup') 666 end 667end 668 669if strfind(thiscomp, 'pc-linux-gnu')==8 670 if exist('LdpcEncode')==0, 671 cd '~/tmp/cml' 672 CmlStartup 673 disp('CmlStartup has been run') 674 % rmpath '/home/bill/cml/mexhelp' % why is this needed? 675 % maybe different path order in octave cf matlab ? 676 end 677end 678 679cd(currentdir) 680 681ldpc_fsk_lib; 682randn('state',1); 683rand('state',1); 684 685% ------------------ select which demo/test to run here --------------- 686 687demo = 12; 688 689if demo == 1 690 printf("simple_ut....\n"); 691 data = simple_ut; 692end 693 694if demo == 2 695 printf("generate C header file....\n"); 696 data = simple_ut("../src/H2064_516_sparse.h"); 697end 698 699if demo == 3 700 printf("decode_from_file ......\n"); 701 data = simple_ut; 702 detected_data = decode_from_file("codeword.bin"); 703 error_positions = xor( detected_data(1:length(data)), data ); 704 Nerrs = sum(error_positions); 705 printf(" Nerrs = %d\n", Nerrs); 706end 707 708if demo == 4 709 printf("plot a curve....\n"); 710 plot_curve; 711end 712 713if demo == 5 714 715 % generate test data and save to disk 716 717 [data code_param] = simple_ut; 718 f = fopen("dat_in2064.bin","wb"); fwrite(f, data, "uint8"); fclose(f); 719 720 % Outboard C encoder 721 722 system("../src/ldpc_enc dat_in2064.bin dat_op2064.bin"); 723 724 % Test with Octave decoder 725 726 detected_data = decode_from_file("dat_op2064.bin"); 727 error_positions = xor(detected_data(1:length(data)), data); 728 Nerrs = sum(error_positions); 729 printf("Nerrs = %d\n", Nerrs); 730end 731 732if demo == 6 733 test_c_encoder; 734end 735 736if demo == 7 737 test_c_decoder; 738end 739 740% generates simulated demod soft decision symbols to drive C ldpc decoder with 741 742if demo == 8 743 frames = 100; 744 EsNodB = 3; 745 EsNo = 10^(EsNodB/10); 746 variance = 1/(2*EsNo); 747 748 frame_rs232 = []; 749 for i=1:frames 750 frame_rs232 = [frame_rs232 gen_sstv_frame]; 751 end 752 753 % write hard decn version to disk file, useful for fsk_mod input 754 755 f = fopen("sstv.bin","wb"); fwrite(f, frame_rs232, "char"); fclose(f); 756 757 % soft decision version (with noise) 758 759 s = 1 - 2*frame_rs232; 760 noise = sqrt(variance)*randn(1,length(frame_rs232)); 761 r = s + noise; 762 f = fopen("sstv_sd.bin","wb"); fwrite(f, r, "float32"); fclose(f); 763end 764 765 766if demo == 9 767 frames = 100; 768 EbNodB = 11; 769 770 frame_rs232 = []; 771 for i=1:frames 772 frame_rs232 = [frame_rs232 gen_sstv_frame]; 773 end 774 775 % Use C FSK modulator to generate modulated signal 776 777 f = fopen("sstv.bin","wb"); fwrite(f, frame_rs232, "char"); fclose(f); 778 system("../build_linux/src/fsk_mod 2 9600 1200 1200 2400 sstv.bin fsk_mod.raw"); 779 780 % Add some channel noise here in Octave 781 782 f = fopen("fsk_mod.raw","rb"); tx = fread(f, "short")'; fclose(f); tx_pwr = var(tx); 783 Fs = 9600; Rs=1200; EbNolin = 10 ^ (EbNodB/10); 784 variance = (tx_pwr/2)*states.Fs/(states.Rs*EbNolin*states.bitspersymbol); 785 noise = sqrt(variance)*randn(1,length(tx)); 786 SNRdB = 10*log10(var(tx)/var(noise)); 787 rx = tx + noise; 788 f = fopen("fsk_demod.raw","wb"); tx = fwrite(f, rx, "short"); fclose(f); 789 790 % Demodulate using C modem and C de-framer/LDPC decoder 791 792 system("../build_linux/src/fsk_demod 2XS 8 9600 1200 fsk_demod.raw - | ../src/drs232_ldpc - dummy_out.bin"); 793end 794 795 796% Plots uncoded BER curves for two different SSTV simulations. Used 797% to compare results with different processing steps as we build up C 798% command line. BER curves are powerful ways to confirm system is 799% operating as expected, several bugs were found using this system. 800 801if demo == 10 802 sim_in.frames = 10; 803 EbNodBvec = 7:10; 804 805 sim_in.demod_type = 3; 806 ber_test1 = []; 807 for i = 1:length(EbNodBvec) 808 [n_uncoded_errs n_uncoded_bits] = run_sstv_sim(sim_in, EbNodBvec(i)); 809 ber_test1(i) = n_uncoded_errs/n_uncoded_bits; 810 end 811 812 sim_in.demod_type = 4; 813 ber_c = []; 814 for i = 1:length(EbNodBvec) 815 [n_uncoded_errs n_uncoded_bits] = run_sstv_sim(sim_in, EbNodBvec(i)); 816 ber_test2(i) = n_uncoded_errs/n_uncoded_bits; 817 end 818 819 figure(1); 820 clf; 821 semilogy(EbNodBvec, ber_test1, '+-;first test;') 822 grid; 823 xlabel('Eb/No (dB)') 824 ylabel('BER') 825 826 hold on; 827 semilogy(EbNodBvec, ber_test2, 'g+-;second test;') 828 legend("boxoff"); 829 hold off; 830 831end 832 833% Measure PER of complete coded and uncoded system 834 835if demo == 11 836 sim_in.frames = 10; 837 EbNodB = 9; 838 sim_in.demod_type = 7; 839 run_sstv_sim(sim_in, EbNodB); 840end 841 842 843% Compare parity bits from an off-air stream of demodulated bits 844% Use something like: 845% cat ~/Desktop/923096fs_wenet.iq | ../build_linux/src/fsk_demod 2X 8 9600 1200 - fsk_demod.bin CU8 846% (note not soft dec mode) 847if demo == 12 848 f = fopen("fsk_demod.bin","rb"); rx_bit_stream = fread(f, "uint8")'; fclose(f); 849 850 compare_parity_bits(rx_bit_stream); 851end 852