From: drowe67 Date: Tue, 9 May 2017 20:53:07 +0000 (+0000) Subject: brought ldpcut back to a simplified demo for blog post X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=6327a80195bd833921c9afc0601dc11f9389c361;p=freetel-svn-tracking.git brought ldpcut back to a simplified demo for blog post git-svn-id: https://svn.code.sf.net/p/freetel/code@3129 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/octave/ldpcut.m b/codec2-dev/octave/ldpcut.m index 18c7ca9b..cf45bd2d 100644 --- a/codec2-dev/octave/ldpcut.m +++ b/codec2-dev/octave/ldpcut.m @@ -2,18 +2,10 @@ % % David Rowe 18 Dec 2013 % -% Octave LDPC unit test script, based on simulation by Bill Cowley VK5DSP -% - -% Start CML library (see CML set up instructions in ldpc.m) - -currentdir = pwd; -addpath '/home/david/Desktop/cml/mat' % assume the source files stored here -cd /home/david/Desktop/cml -CmlStartup % note that this is not in the cml path! -cd(currentdir) +% Octave LDPC unit test script using CML library, based on simulation +% by Bill Cowley VK5DSP -% Our LDPC library +% Libraries we need ldpc; qpsk; @@ -26,21 +18,10 @@ function sim_out = run_simulation(sim_in) 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 - + framesize = sim_in.framesize; + rate = sim_in.rate; Ntrials = sim_in.Ntrials; 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; % Init LDPC code ------------------------------------ @@ -52,53 +33,10 @@ function sim_out = run_simulation(sim_in) decoder_type = 0; max_iterations = 100; - 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 optional HF (multipath) model ------------------------------------ - - fading = ones(1,Ntrials*code_param.code_bits_per_frame/bps); - - if hf_en - - % 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") - dopplerSpreadHz = sim_in.dopplerSpreadHz; - end - if isfield(sim_in, "path_delay") - path_delay = sim_in.path_delay; - end - printf("Doppler Spread: %3.2f Hz Path Delay: %3.2f symbols\n", dopplerSpreadHz, path_delay); - - % reset seed so we get same fading channel on every simulation - - 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, 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 + code_param = ldpc_init_wimax(rate, framesize, modulation, mod_order, mapping); % ---------------------------------- - % run simulation at each EB/No point + % run simulation at each Eb/No point % ---------------------------------- for ne = 1:length(EbNodBvec) @@ -116,7 +54,6 @@ function sim_out = run_simulation(sim_in) EsNodB = EbNodBvec(ne) + 10*log10(rate) + 10*log10(bps); EsNo = 10^(EsNodB/10); variance = 1/EsNo; - hf_r = 1; Tbits = Terrs = Ferrs = Terrs_raw = Tbits_raw = 0; @@ -135,30 +72,6 @@ function sim_out = run_simulation(sim_in) rx_symbols = tx_symbols; - % Optional HF (multipath) channel model - - if hf_en - - % 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 splits power evenly between Re & Im noise = sqrt(variance*0.5)*(randn(1,length(tx_symbols)) + j*randn(1,length(tx_symbols))); @@ -166,25 +79,28 @@ function sim_out = run_simulation(sim_in) % Decode a bunch of frames + rx_bits_log = []; + for nn = 1: Ntrials st = (nn-1)*code_param.symbols_per_frame + 1; en = (nn)*code_param.symbols_per_frame; % coded - arx_codeword = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, rx_symbols(st:en), EsNo, fading(st:en)); + arx_codeword = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, rx_symbols(st:en), EsNo); 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)); Nerrs = sum( error_positions); + rx_bits_log = [rx_bits_log arx_codeword(1:code_param.data_bits_per_frame)]; % uncoded - to est raw BER compare first half or received frame to tx_bits as code is systematic - rx_bits = []; + raw_rx_bits = []; for s=1:code_param.symbols_per_frame*rate - rx_bits = [rx_bits qpsk_demod(rx_symbols(st+s-1))]; + raw_rx_bits = [raw_rx_bits qpsk_demod(rx_symbols(st+s-1))]; end - Nerrs_raw = sum(xor(rx_bits, tx_bits(st:en))); + Nerrs_raw = sum(xor(raw_rx_bits, tx_bits(st:en))); Nbits_raw = code_param.data_bits_per_frame; if verbose == 2 @@ -210,53 +126,35 @@ function sim_out = run_simulation(sim_in) 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 + sim_out.PER(ne) = Ferrs/Ntrials; end endfunction - -% START SIMULATIONS --------------------------------------------------------------- +% -------------------------------------------------------------------------------- +% START SIMULATIONS +% -------------------------------------------------------------------------------- more off; -format +format; -% set to 1 to use wimax codes built into CML library +% --------------------------------------------------------------------------------- +% Start CML library (see CML set up instructions in ldpc.m) +% --------------------------------------------------------------------------------- + +currentdir = pwd; +addpath '/home/david/Desktop/cml/mat' % assume the source files stored here +cd /home/david/Desktop/cml +CmlStartup % note that this is not in the cml path! +cd(currentdir) -wimax_en = 1; % --------------------------------------------------------------------------------- % 1/ Simplest possible one frame simulation % --------------------------------------------------------------------------------- -printf("Test 1\n------\n"); +printf("\nTest 1\n------\n"); mod_order = 4; modulation = 'QPSK'; @@ -265,14 +163,9 @@ 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 +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); EsNo = 10; % decoder needs an estimated channel EsNo (linear ratio, not dB) @@ -291,17 +184,14 @@ printf("Nerr: %d\n\n", Nerr); printf("Test 2\n------\n"); -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 - % e.g. try *4 and note difference in HF perf -end +% these are inputs for Wimax mode, e.g. framesize defines code used + +sim_in.rate = 0.5; +sim_in.framesize = 576*4; % long codes smooth over fades but increase latency + sim_in.verbose = 2; sim_in.Ntrials = 100; sim_in.EbNodBvec = 9; -sim_in.hf_en = 0; run_simulation(sim_in); @@ -312,33 +202,20 @@ run_simulation(sim_in); printf("\n\nTest 3\n------\n"); 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.EbNodBvec; -uncoded_theory = 0.5*erfc(sqrt(10.^(EbNodB/10))); - -EbNoLin = 10.^(EbNodB/10); -uncoded_hf_theory = 0.5.*(1-sqrt(EbNoLin./(EbNoLin+1))); +uncoded_awgn_ber_theory = 0.5*erfc(sqrt(10.^(EbNodB/10))); figure(1) clf -semilogy(EbNodB, uncoded_theory,'r-+;AWGN;') +semilogy(EbNodB, uncoded_awgn_ber_theory,'r-+;AWGN;') hold on; -semilogy(EbNodB, uncoded_hf_theory,'r-o;HF;'); 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_theory) 1]) +axis([min(EbNodB) max(EbNodB) 1E-3 1]) legend('boxoff') - - - - -