%
% 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;
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 ------------------------------------
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)
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;
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)));
% 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
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';
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)
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);
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')
-
-
-
-
-
projects / freetel-svn-tracking.git / commitdiff