1% fsk_lib_ldpc_demo.m 2% 3% LDPC coded 4FSK modem demo, demonstrating soft dec using CML library functions 4 5fsk_lib; 6ldpc; 7 8% set up waveform 9function [states M] = modem_init(Rs,Fs,df) 10 M = 4; 11 states = fsk_init(Fs,Rs,M,P=8,nsym=100); 12 states.tx_real = 0; 13 states.tx_tone_separation = 250; 14 states.ftx = -2.5*states.tx_tone_separation + states.tx_tone_separation*(1:M); 15 states.fest_fmin = -Fs/2; 16 states.fest_fmax = +Fs/2; 17 states.fest_min_spacing = Rs/2; 18 states.df = df; 19 20 states.ber_valid_thresh = 0.1; 21 states.ber_invalid_thresh = 0.2; 22end 23 24% Run a complete modem (freq and timing estimators running) at a 25% single Eb/No point. At low Eb/No the estimators occasionally fall 26% over so we get complete junk, we consider that case a packet error 27% and exclude it from the BER estimation. 28 29function [states uber cber cper] = modem_run_test(HRA, EbNodB = 10, num_frames=10, Fs=8000, Rs=100, df=0, plots=0) 30 rand('seed',1); randn('seed',1); 31 [states M] = modem_init(Rs, Fs, df); 32 N = states.N; 33 if plots; states.verbose = 0x4; end 34 35 % set up LDPC code 36 Hsize=size(HRA); 37 Krate = (Hsize(2)-Hsize(1))/Hsize(2); states.rate = Krate; 38 code_param = ldpc_init_user(HRA, modulation='FSK', mod_order=states.M, mapping='gray'); 39 states.coden = code_param.coded_bits_per_frame; 40 states.codek = code_param.data_bits_per_frame; 41 42 % set up AWGN noise 43 EcNodB = EbNodB + 10*log10(Krate); 44 EcNo = 10^(EcNodB/10); 45 variance = states.Fs/(states.Rs*EcNo*states.bitspersymbol); 46 47 data_bits = round(rand(1,code_param.data_bits_per_frame)); tx_bits = []; 48 for f=1:num_frames 49 codeword_bits = LdpcEncode(data_bits, code_param.H_rows, code_param.P_matrix); 50 tx_bits = [tx_bits codeword_bits]; 51 end 52 53 % modulator and AWGN channel 54 tx = fsk_mod(states, tx_bits); 55 noise = sqrt(variance/2)*randn(length(tx),1) + j*sqrt(variance/2)*randn(length(tx),1); 56 rx = tx + noise; 57 58 % freq estimator and demod 59 run_frames = floor(length(rx)/N)-1; 60 st = 1; f_log = []; rx_bits = []; rx_filt = []; 61 for f=1:run_frames 62 63 % extract nin samples from input stream 64 nin = states.nin; 65 en = st + states.nin - 1; 66 67 % due to nin variations it's possible to overrun buffer 68 if en < length(rx) 69 sf = rx(st:en); 70 states = est_freq(states, sf, states.M); states.f = states.f2; 71 [arx_bits states] = fsk_demod(states, sf); 72 rx_bits = [rx_bits arx_bits]; 73 rx_filt = [rx_filt abs(states.f_int_resample)]; 74 f_log = [f_log; states.f]; 75 st += nin; 76 end 77 end 78 79 % count bit errors in test frames 80 81 num_frames=floor(length(rx_bits)/code_param.coded_bits_per_frame); 82 log_nerrs = []; num_frames_rx = 0; Tbits = Terrs = Tperr = Tpackets = 0; 83 uber = cber = 0.5; cper = 1; 84 for f=1:num_frames-1 85 st = (f-1)*code_param.coded_bits_per_frame + 1; en = (f+1)*code_param.coded_bits_per_frame; 86 states = ber_counter(states, codeword_bits, rx_bits(st:en)); 87 log_nerrs = [log_nerrs states.nerr]; 88 if states.ber_state num_frames_rx++; end 89 90 % Using sync provided by ber_counter() state machine for LDPC frame alignment 91 if states.ber_state 92 st_bit = (f-1)*code_param.coded_bits_per_frame + states.coarse_offset; 93 st_symbol = (st_bit-1)/states.bitspersymbol + 1; 94 en_symbol = st_symbol + code_param.coded_bits_per_frame/states.bitspersymbol - 1; 95 %printf("coded_bits: %d bps: %d st_bit: %d st_symbol: %d en_symbol: %d\n", 96 %code_param.coded_bits_per_frame, states.bitspersymbol, st_bit, st_symbol, en_symbol); 97 98 % map FSK filter ouputs to LLRs, then LDPC decode (see also fsk_cml_sam.m) 99 symL = DemodFSK(1/states.v_est*rx_filt(:,st_symbol:en_symbol), states.SNRest, 1); 100 llr = -Somap(symL); 101 [x_hat, PCcnt] = MpDecode(llr, code_param.H_rows, code_param.H_cols, max_iterations=100, decoder_type=0, 1, 1); 102 Niters = sum(PCcnt~=0); 103 detected_data = x_hat(Niters,:); 104 Nerrs = sum(xor(data_bits, detected_data(1:code_param.data_bits_per_frame))); 105 Terrs += Nerrs; 106 Tbits += code_param.data_bits_per_frame; 107 if Nerrs Tperr++; end 108 Tpackets++; 109 end 110 end 111 112 if states.Terrs 113 printf("Fs: %d Rs: %d df % 3.2f EbNo: %4.2f ftx: %3d frx: %3d\n",Fs, Rs, df, EbNodB, num_frames, num_frames_rx); 114 uber = states.Terrs/states.Tbits; cber = Terrs/Tbits; cper = Tperr/Tpackets; 115 printf(" Uncoded: nbits: %6d nerrs: %6d ber: %4.3f\n", states.Tbits, states.Terrs, uber); 116 printf(" Coded..: nbits: %6d nerrs: %6d ber: %4.3f\n", Tbits, Terrs, cber); 117 printf(" Coded..: npckt: %6d perrs: %6d per: %4.3f\n", Tpackets, Tperr, cper); 118 end 119 120 if plots 121 figure(1); clf; 122 ideal=ones(length(f_log),1)*states.ftx; 123 plot((1:length(f_log)),ideal(:,1),'bk;ideal;') 124 hold on; plot((1:length(f_log)),ideal(:,2:states.M),'bk'); hold off; 125 hold on; 126 plot(f_log(:,1), 'linewidth', 2, 'b;peak;'); 127 plot(f_log(:,2:states.M), 'linewidth', 2, 'b'); 128 hold off; 129 xlabel('Time (frames)'); ylabel('Frequency (Hz)'); 130 figure(2); clf; plot(log_nerrs); title('Errors per frame'); 131 end 132 133end 134 135 136function freq_run_curve_peak_mask(HRA, num_frames=100) 137 138 EbNodB = 4:10; 139 m4fsk_ber_theory = [0.23 0.18 0.14 0.09772 0.06156 0.03395 0.01579 0.00591 0.00168 3.39E-4]; 140 uber_log = []; cber_log = []; cper_log = []; 141 for ne = 1:length(EbNodB) 142 [states uber cber cper] = modem_run_test(HRA, EbNodB(ne), num_frames); 143 uber_log = [uber_log uber]; cber_log = [cber_log cber]; cper_log = [cper_log cper]; 144 end 145 146 figure(1); clf; 147 EbNodB_raw = EbNodB+10*log10(states.rate) 148 semilogy(EbNodB_raw, m4fsk_ber_theory(round(EbNodB_raw+1)), 'linewidth', 2, 'bk+-;uber theory;'); 149 grid; hold on; 150 semilogy(EbNodB_raw, uber_log+1E-12, 'linewidth', 2, '+-;uber;'); 151 semilogy(EbNodB, cber_log+1E-12, 'linewidth', 2, 'r+-;cber;'); 152 semilogy(EbNodB, cper_log+1E-12, 'linewidth', 2, 'c+-;cper;'); hold off; 153 xlabel('Eb/No (info bits, dB)'); ylabel('BER/PER'); axis([min(EbNodB_raw) max(EbNodB) 1E-4 1]); 154 title(sprintf("%dFSK rate %3.1f (%d,%d) Ncodewords=%d NCodewordBits=%d Fs=%d Rs=%d", 155 states.M, states.rate, states.coden, states.codek, num_frames, states.Tbits, states.Fs, states.Rs)); 156 print("fsk_lib_ldpc.png", "-dpng") 157end 158 159% Choose simulation here --------------------------------------------------- 160 161init_cml(); 162load H_256_512_4.mat; HRA=H; 163more off; 164 165% single point 166[states uber cber cper] = modem_run_test(HRA, EbNodB=8); 167if cber == 0 168 printf("PASS\n"); 169end 170 171% curve 172%freq_run_curve_peak_mask(HRA, 200) 173