From bf399cdc1be3e78b1733b40298031725e9aed5d3 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Sun, 7 May 2017 20:23:40 +0000 Subject: [PATCH] modified to test Bill's short LDPC codes git-svn-id: https://svn.code.sf.net/p/freetel/code@3125 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/ldpcut.m | 183 +++++++++++++++++++++++++++---------- 1 file changed, 137 insertions(+), 46 deletions(-) diff --git a/codec2-dev/octave/ldpcut.m b/codec2-dev/octave/ldpcut.m index 0cb2cadf..b78eb1ea 100644 --- a/codec2-dev/octave/ldpcut.m +++ b/codec2-dev/octave/ldpcut.m @@ -20,18 +20,29 @@ qpsk; function sim_out = run_simulation(sim_in) - rate = sim_in.rate; - framesize = sim_in.framesize; + % Note this is effective Eb/No of payload data bits, sorta thing we + % plot on BER versus Eb/No graphs of decoded data. So if we have a + % rate 1/2 code, each codeword bit will have Eb/No - 3dB. + + EbNodBvec = sim_in.EbNodBvec; + + % for wimax code frame size specifies code + + if isfield(sim_in, "framesize") + framesize = sim_in.framesize; + rate = sim_in.rate; + end + Ntrials = sim_in.Ntrials; - EsNodBvec = sim_in.EsNodBvec; verbose = sim_in.verbose; if isfield(sim_in, "hf_en") hf_en = sim_in.hf_en; else hf_en = 0; end + wimax_en = sim_in.wimax_en; - % Start simulation + % Init LDPC code ------------------------------------ mod_order = 4; bps = 2; modulation = 'QPSK'; @@ -41,15 +52,22 @@ function sim_out = run_simulation(sim_in) decoder_type = 0; max_iterations = 100; - code_param = ldpc_init(rate, framesize, modulation, mod_order, mapping); + if wimax_en + code_param = ldpc_init_wimax(rate, framesize, modulation, mod_order, mapping); + else + load HRA_112_112.txt + [code_param framesize rate] = ldpc_init_user(HRA_112_112, modulation, mod_order, mapping); + end + + % set up option HF (multipath) model ------------------------------------ - % set up HF model + fading = ones(1,Ntrials*code_param.code_bits_per_frame/bps); if hf_en - % some typical values, or replace with user supplied - % Note we need to specify a symbol rate Rs + % We assume symbols spread acroos Nc OFDM carriers + Nc = 14; Rs = 50; dopplerSpreadHz = 1.0; path_delay = 1E-3*Rs; if isfield(sim_in, "dopplerSpreadHz") @@ -64,27 +82,43 @@ function sim_out = run_simulation(sim_in) randn('seed',1); + Ns = Ntrials*code_param.code_bits_per_frame/bps; + Nr = ceil(Ns/Nc); + hf_model = zeros(Nr,Nc); + % note we ask for 10% more samples than we need, as % doppler_spread() function can sometimes return slightly less % than we need due to round off - spread1 = doppler_spread(dopplerSpreadHz, Rs, Ntrials*framesize*1.1); - spread2 = doppler_spread(dopplerSpreadHz, Rs, Ntrials*framesize*1.1); - spread1 = spread1(1:Ntrials*framesize/bps); - spread2 = spread2(1:Ntrials*framesize/bps); + spread1 = doppler_spread(dopplerSpreadHz, Rs, Nr*1.1); + spread2 = doppler_spread(dopplerSpreadHz, Rs, Nr*1.1); + spread1 = spread1(1:Nr); + spread2 = spread2(1:Nr); hf_gain = 1.0/sqrt(var(spread1)+var(spread2)); end - % --------------------------------- - % run simulation at each EsNo point - % --------------------------------- + % ---------------------------------- + % run simulation at each EB/No point + % ---------------------------------- - for ne = 1:length(EsNodBvec) - EsNodB = EsNodBvec(ne); + for ne = 1:length(EbNodBvec) + randn('seed',1); + rand('seed',1); + + % Given Eb/No of payload data bits, work out Es/No we need to + % apply to each channel symbol: + % + % i) Each codeword bit gets noise: Eb/No - 3 (for a rate 1/2 code) + % ii) QPSK means two bits/symbol.: Es/No = Eb/No + 3 + % + % -> which neatly cancel out ...... (at least for rate 1/2) + + EsNodB = EbNodBvec(ne) + 10*log10(rate) + 10*log10(bps); EsNo = 10^(EsNodB/10); variance = 1/EsNo; - - Tbits = Terrs = Ferrs = Terrs_raw = 0; + hf_r = 1; + + Tbits = Terrs = Ferrs = Terrs_raw = Tbits_raw = 0; tx_bits = []; tx_symbols = []; @@ -101,22 +135,34 @@ function sim_out = run_simulation(sim_in) rx_symbols = tx_symbols; - hf_model = ones(1,length(tx_symbols)); - if hf_en + % Optional HF (multipath) channel model - % Simplified rate Rs simulation model of single carrier. Note - % if freq of carrier is 0, model can be simplified further, as - % path delay term collapses. + if hf_en - w = 0; - hf_model = hf_gain*(spread1 + exp(j*w*path_delay)*spread2); - rx_symbols = tx_symbols .* abs(hf_model); + % Simplified rate Rs OFDM simulation model that doesn't + % include ISI, just freq filtering. We assume perfect phase + % estimation so it's just amplitude distortion. We distribute + % symbols across Nc carriers + + Ns = length(tx_symbols); + w = (1:Nc)*2*pi; + rx_symbols = [rx_symbols zeros(1,Nr*Nc-Ns)]; % pad out to integer number of rows + + for r=1:Nr + for c=1:Nc + hf_model(hf_r,c) = hf_gain*(spread1(hf_r) + exp(-j*w(c)*path_delay)*spread2(hf_r)); + rx_symbols(Nc*(r-1)+c) *= abs(hf_model(hf_r,c)); + fading(Nc*(r-1)+c) = abs(hf_model(hf_r,c)); + end + hf_r++; + end + rx_symbols = rx_symbols(1:Ns); % remove padding end - % Add AWGN noise, 0.5 factor ensures var(noise) == variance , i.e. splits power between Re & Im + % Add AWGN noise, 0.5 factor splits power evenly between Re & Im noise = sqrt(variance*0.5)*(randn(1,length(tx_symbols)) + j*randn(1,length(tx_symbols))); - rx_symbols += + noise; + rx_symbols += noise; % Decode a bunch of frames @@ -126,7 +172,7 @@ function sim_out = run_simulation(sim_in) % coded - arx_codeword = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, rx_symbols(st:en), EsNo, abs(hf_model(st:en))); + arx_codeword = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, rx_symbols(st:en), EsNo, fading(st:en)); st = (nn-1)*code_param.data_bits_per_frame + 1; en = (nn)*code_param.data_bits_per_frame; error_positions = xor(arx_codeword(1:code_param.data_bits_per_frame), tx_bits(st:en)); @@ -139,6 +185,7 @@ function sim_out = run_simulation(sim_in) rx_bits = [rx_bits qpsk_demod(rx_symbols(st+s-1))]; end Nerrs_raw = sum(xor(rx_bits, tx_bits(st:en))); + Nbits_raw = code_param.data_bits_per_frame; if verbose == 2 % print "." if frame decoded without errors, 'x' if we can't decode @@ -149,20 +196,48 @@ function sim_out = run_simulation(sim_in) if Nerrs > 0, Ferrs = Ferrs + 1; end Terrs += Nerrs; - Terrs_raw += Nerrs_raw; Tbits += code_param.data_bits_per_frame; + Terrs_raw += Nerrs_raw; + Tbits_raw += Nbits_raw; end if verbose - printf("\nEsNodB: %3.2f BER: %f Tbits: %d Terrs: %d Ferrs: %d Terrs_raw: %d Raw BER: %3.2f\n", - EsNodB, Terrs/Tbits, Tbits, Terrs, Ferrs, Terrs_raw, Terrs_raw/Tbits) + printf("\nCoded EbNodB: %3.2f BER: %4.3f Tbits: %6d Terrs: %6d FER: %4.3f Tframes: %d Ferrs: %d\n", + EsNodB, Terrs/Tbits, Tbits, Terrs, Ferrs/Ntrials, Ntrials, Ferrs); + EbNodB_raw = EbNodBvec(ne) + 10*log10(rate); + printf("Raw EbNodB..: %3.2f BER: %4.3f Tbits: %6d Terrs: %6d\n", + EbNodB_raw, Terrs_raw/Tbits_raw, Tbits_raw, Terrs_raw); end + sim_out.rate = rate; sim_out.Tbits(ne) = Tbits; sim_out.Terrs(ne) = Terrs; sim_out.Ferrs(ne) = Ferrs; sim_out.BER(ne) = Terrs/Tbits; sim_out.FER(ne) = Ferrs/Ntrials; + + if hf_en && (verbose > 1) + figure(2); clf; + plot(real(rx_symbols(Ns/2:Ns)), imag(rx_symbols(Ns/2:Ns)), '+'); + axis([-2 2 -2 2]); + title('Scatter') + + figure(3); clf; + + % limit mesh plot to Np points to plot quickly + + Np = 500; + step = ceil(hf_r/Np); + mesh(1:Nc, (1:step:hf_r-1)/Rs, abs(hf_model(1:step:hf_r-1,:))) + title('HF channel amplitude'); + xlabel('Carrier'); + ylabel('Time (s)'); + + figure(4) + subplot(211); plot(abs(spread1)); + subplot(212); plot(abs(spread2)); + title('HF spreading function amplitudes') + end end endfunction @@ -173,6 +248,8 @@ endfunction more off; format +wimax_en = 1; + % --------------------------------------------------------------------------------- % 1/ Simplest possible one frame simulation % --------------------------------------------------------------------------------- @@ -182,14 +259,21 @@ printf("Test 1\n------\n"); mod_order = 4; modulation = 'QPSK'; mapping = 'gray'; -framesize = 576; % CML library has a bunch of different framesizes available -rate = 1/2; demod_type = 0; decoder_type = 0; max_iterations = 100; + +if wimax_en + framesize = 576*2; % CML library has a bunch of different framesizes available + rate = 1/2; + code_param = ldpc_init_wimax(rate, framesize, modulation, mod_order, mapping); +else + load HRA_112_112.txt + [code_param framesize rate] = ldpc_init_user(HRA_112_112, modulation, mod_order, mapping); +end + EsNo = 10; % decoder needs an estimated channel EsNo (linear ratio, not dB) -code_param = ldpc_init(rate, framesize, modulation, mod_order, mapping); tx_bits = round(rand(1, code_param.data_bits_per_frame)); [tx_codeword, qpsk_symbols] = ldpc_enc(tx_bits, code_param); rx_codeword = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, qpsk_symbols, EsNo); @@ -205,12 +289,16 @@ printf("Nerr: %d\n\n", Nerr); printf("Test 2\n------\n"); -sim_in.rate = 0.5; -sim_in.framesize = 2*576; +sim_in.wimax_en = wimax_en; +if sim_in.wimax_en + % these are inputs for Wimax mode, e.g. framesize defines code used + sim_in.rate = 0.5; + sim_in.framesize = 576*2; % long codes smooth over fades but increase latency +end sim_in.verbose = 2; -sim_in.Ntrials = 100; -sim_in.EsNodBvec = 2; -sim_in.hf_en = 0; +sim_in.Ntrials = 1000/4; +sim_in.EbNodBvec = 9; +sim_in.hf_en = 1; run_simulation(sim_in); @@ -220,13 +308,13 @@ run_simulation(sim_in); printf("\n\nTest 3\n------\n"); -sim_in.EsNodBvec = -2:10; +sim_in.EbNodBvec = -2:10; sim_in.hf_en = 0; sim_out = run_simulation(sim_in); sim_in.hf_en = 1; sim_out_hf = run_simulation(sim_in); -EbNodB = sim_in.EsNodBvec - 10*log10(2); % QPSK has two bits/symbol +EbNodB = sim_in.EbNodBvec; uncoded_theory = 0.5*erfc(sqrt(10.^(EbNodB/10))); EbNoLin = 10.^(EbNodB/10); @@ -237,14 +325,17 @@ clf semilogy(EbNodB, uncoded_theory,'r-+;AWGN;') hold on; semilogy(EbNodB, uncoded_hf_theory,'r-o;HF;'); -semilogy(EbNodB-10*log10(sim_in.rate), sim_out.BER+1E-10,'g-+;AWGN LDPC;'); -semilogy(EbNodB-10*log10(sim_in.rate), sim_out_hf.BER+1E-10,'g-o;HF LDPC;'); +semilogy(EbNodB, sim_out.BER+1E-10,'g-+;AWGN LDPC;'); +semilogy(EbNodB, sim_out_hf.BER+1E-10,'g-o;HF LDPC;'); hold off; grid('minor') xlabel('Eb/No (dB)') ylabel('BER') -axis([min(EbNodB) max(EbNodB) min(uncoded_BER_theory) 1]) +axis([min(EbNodB) max(EbNodB) min(uncoded_theory) 1]) legend('boxoff') + + + -- 2.25.1