--- /dev/null
+% ldpc_qpsk.m
+%
+% David Rowe 18 Dec 2013
+%
+% Similar to ldpc_short.m, but derived from ldpcut.m and uses QPSK and
+% CML 2D functunctions and QPSK. Probably should combine this and
+% ldpc_short.m some day.
+
+% Our LDPC library
+
+ldpc;
+qpsk;
+gp_interleaver;
+
+
+function sim_out = run_simulation(sim_in)
+
+ % 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;
+ 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;
+ interleave_en = sim_in.interleave_en;
+
+ % Init LDPC code ------------------------------------
+
+ mod_order = 4; bps = 2;
+ modulation = 'QPSK';
+ mapping = 'gray';
+
+ demod_type = 0;
+ 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);
+ Nc = 14; Rs = 50; Tp = (framesize/Nc)/Rs; Tp_codec2 = 0.04;
+
+ if hf_en
+
+ % We assume symbols spread acroos Nc OFDM carriers
+
+ 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
+
+ % ----------------------------------
+ % run simulation at each Eb/No point
+ % ----------------------------------
+
+ 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;
+ hf_r = 1;
+
+ Tbits = Terrs = Ferrs = Terrs_raw = Tbits_raw = 0;
+
+ tx_bits = [];
+ tx_symbols = [];
+ rx_symbols = [];
+
+ % Encode a bunch of frames
+
+ for nn=1:Ntrials
+ atx_bits = round(rand( 1, code_param.data_bits_per_frame));
+ tx_bits = [tx_bits atx_bits];
+ [tx_codeword atx_symbols] = ldpc_enc(atx_bits, code_param);
+ if interleave_en
+ atx_symbols = gp_interleave(atx_symbols);
+ end
+ tx_symbols = [tx_symbols atx_symbols];
+ end
+
+ 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)));
+ rx_symbols += noise;
+
+ % Decode a bunch of frames
+
+ rx_bits_log = []; error_positions_log = [];
+ Nerrs_raw_log = [];
+
+ for nn = 1: Ntrials
+ st = (nn-1)*code_param.symbols_per_frame + 1;
+ en = (nn)*code_param.symbols_per_frame;
+
+ % coded
+
+ arx_symbols = rx_symbols(st:en);
+ afading = fading(st:en);
+ if interleave_en
+ arx_symbols = gp_deinterleave(arx_symbols);
+ afading = gp_deinterleave(afading);
+ end
+
+ arx_codeword = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, arx_symbols, EsNo, afading);
+ st = (nn-1)*code_param.data_bits_per_frame + 1;
+ en = (nn)*code_param.data_bits_per_frame;
+ arx_bits = arx_codeword(1:code_param.data_bits_per_frame);
+ error_positions = xor(arx_bits, tx_bits(st:en));
+ error_positions_log = [error_positions_log error_positions];
+ Nerrs = sum( error_positions);
+ rx_bits_log = [rx_bits_log arx_bits];
+
+ % uncoded - to est raw BER compare first half or received frame to tx_bits as code is systematic
+
+ raw_rx_bits = [];
+ for s=1:code_param.symbols_per_frame*rate
+ raw_rx_bits = [raw_rx_bits qpsk_demod(arx_symbols(s))];
+ end
+ Nerrs_raw = sum(xor(raw_rx_bits, tx_bits(st:en)));
+ Nbits_raw = code_param.data_bits_per_frame;
+ Nerrs_raw_log = [Nerrs_raw_log Nerrs_raw];
+
+ if verbose == 2
+ % print "." if frame decoded without errors, 'x' if we can't decode
+
+ if Nerrs > 0, printf('x'), else printf('.'), end
+ if (rem(nn, 50)==0), printf('\n'), end
+ end
+
+ if Nerrs > 0, Ferrs = Ferrs + 1; end
+ Terrs += Nerrs;
+ Tbits += code_param.data_bits_per_frame;
+ Terrs_raw += Nerrs_raw;
+ Tbits_raw += Nbits_raw;
+ end
+
+ % Alternative Codec 2 packet rate measurement indep of framesize
+
+ Nerrs_codec2_log = []; Ncodecpacketsize = 28;
+ Perrs = 0; Npackets = floor(length(tx_bits)/Ncodecpacketsize);
+ for p=1:Ncodecpacketsize:Npackets*Ncodecpacketsize
+ Nerrs = sum(xor(tx_bits(p:p+Ncodecpacketsize-1), rx_bits_log(p:p+Ncodecpacketsize-1)));
+ if Nerrs
+ Perrs++;
+ end
+ Nerrs_codec2_log = [Nerrs_codec2_log Nerrs];
+ end
+
+ if verbose
+ 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);
+ printf("Codec 2 PER: %5.4f Npackets: %d Perrs: %d\n", Perrs/Npackets, Npackets, Perrs);
+ end
+
+ sim_out.rate = rate;
+ sim_out.BER(ne) = Terrs/Tbits;
+ sim_out.PER(ne) = Perrs/Npackets;
+ sim_out.error_positions = error_positions_log;
+
+ 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;
+ subplot(211);
+ stem((1:Ntrials)*Tp, Nerrs_raw_log);
+ subplot(212);
+ stem((1:Npackets)*Tp_codec2, Nerrs_codec2_log);
+
+ figure(4); 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(5)
+ subplot(211); plot((1:hf_r-1)/Rs, abs(spread1(1:hf_r-1)));
+ subplot(212); plot((1:hf_r-1)/Rs, abs(spread2(1:hf_r-1)));
+ title('HF spreading function amplitudes')
+ end
+ end
+
+endfunction
+
+
+% ---------------------------------------------------------------------------------
+% Run a bunch of trials at just one EbNo point
+% ---------------------------------------------------------------------------------
+
+function run_single(Nbits=700*10, EbNodB=9, hf_en=0, wimax_en=1, framesize=576, interleave_en=0, error_pattern_filename)
+ sim_in.wimax_en = wimax_en;
+ if sim_in.wimax_en
+ % CML wimax codes
+ sim_in.rate = 0.5;
+ sim_in.framesize = framesize;
+ else
+ % Our HRA short LDPC code
+ sim_in.rate=0.5;
+ sim_in.framesize=224;
+ end
+ sim_in.verbose = 2;
+ sim_in.Ntrials = floor(Nbits/(sim_in.framesize*sim_in.rate));
+ sim_in.EbNodBvec = EbNodB;
+ sim_in.hf_en = hf_en;
+ sim_in.interleave_en = interleave_en;
+
+ sim_out = run_simulation(sim_in);
+
+ if nargin == 7
+ fep = fopen(error_pattern_filename, "wb");
+ fwrite(fep, sim_out.error_positions, "short");
+ fclose(fep);
+ end
+
+end
+
+% ---------------------------------------------------------------------------------
+% Lets draw some Eb/No versus BER curves
+% ---------------------------------------------------------------------------------
+
+function plot_curves(Nbits=700*60)
+ sim_in.EbNodBvec = -2:10;
+ sim_in.verbose = 2;
+ sim_in.interleave_en = 2;
+
+ % Wimax codes
+
+ sim_in.wimax_en = 1;
+ sim_in.rate = 0.5;
+ sim_in.framesize = 576*4;
+ sim_in.Ntrials = floor(Nbits/(sim_in.framesize*sim_in.rate));
+
+ sim_in.hf_en = 0; sim_out_awgn_wimax = run_simulation(sim_in);
+ sim_in.hf_en = 1; sim_out_hf_wimax = run_simulation(sim_in);
+
+ % Our short HRA codes
+
+ sim_in.wimax_en = 0;
+ sim_in.rate = 0.5;
+ sim_in.framesize = 224;
+ sim_in.Ntrials = floor(Nbits/(sim_in.framesize*sim_in.rate));
+
+ sim_in.hf_en = 0; sim_out_awgn_short = run_simulation(sim_in);
+ sim_in.hf_en = 1; sim_out_hf_short = run_simulation(sim_in);
+
+ EbNodB = sim_in.EbNodBvec;
+ uncoded_awgn_ber_theory = 0.5*erfc(sqrt(10.^(EbNodB/10)));
+
+ EbNoLin = 10.^(EbNodB/10);
+ uncoded_hf_ber_theory = 0.5.*(1-sqrt(EbNoLin./(EbNoLin+1)));
+
+ figure(1)
+ clf
+ semilogy(EbNodB, uncoded_awgn_ber_theory,'r-+;AWGN Uncoded;','markersize', 10, 'linewidth', 2)
+ hold on;
+ semilogy(EbNodB, uncoded_hf_ber_theory,'r-o;HF Uncoded;','markersize', 10, 'linewidth', 2);
+ semilogy(EbNodB, sim_out_awgn_wimax.BER+1E-10,'g-+;AWGN LDPC (2304,1152);','markersize', 10, 'linewidth', 2);
+ semilogy(EbNodB, sim_out_hf_wimax.BER+1E-10,'g-o;HF LDPC (2304,1152);','markersize', 10, 'linewidth', 2);
+ semilogy(EbNodB, sim_out_awgn_short.BER+1E-10,'b-+;AWGN LDPC (224,112);','markersize', 10, 'linewidth', 2);
+ semilogy(EbNodB, sim_out_hf_short.BER+1E-10,'b-o;HF LDPC (224,112);','markersize', 10, 'linewidth', 2);
+ hold off;
+ grid('minor')
+ xlabel('Eb/No (dB)')
+ ylabel('BER')
+ axis([min(EbNodB) max(EbNodB) 1E-3 5E-1])
+ legend('boxoff')
+
+ uncoded_awgn_per_theory = 1 - (1-uncoded_awgn_ber_theory).^28;
+ uncoded_hf_per_theory = 1 - (1-uncoded_hf_ber_theory).^28;
+
+ figure(2)
+ clf
+ semilogy(EbNodB, uncoded_awgn_per_theory,'r-+;AWGN Uncoded;','markersize', 10, 'linewidth', 2)
+ hold on;
+ semilogy(EbNodB, uncoded_hf_per_theory,'r-o;HF Uncoded;','markersize', 10, 'linewidth', 2);
+ semilogy(EbNodB, sim_out_awgn_wimax.PER+1E-10,'g-+;AWGN LDPC (2304,1152);','markersize', 10, 'linewidth', 2);
+ semilogy(EbNodB, sim_out_hf_wimax.PER+1E-10,'g-o;HF LDPC (2304,1152);','markersize', 10, 'linewidth', 2);
+ semilogy(EbNodB, sim_out_awgn_short.PER+1E-10,'b-+;AWGN LDPC (224,112);','markersize', 10, 'linewidth', 2);
+ semilogy(EbNodB, sim_out_hf_short.PER+1E-10,'b-o;HF LDPC (224,112);','markersize', 10, 'linewidth', 2);
+ hold off;
+ grid('minor')
+ xlabel('Eb/No (dB)')
+ ylabel('PER')
+ axis([min(EbNodB) max(EbNodB) 1E-2 1])
+ legend('boxoff')
+ legend("location", "southwest");
+end
+
+
+% ---------------------------------------------------------------------------------
+% Start simulations here
+% ---------------------------------------------------------------------------------
+
+more off;
+format;
+
+% 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)
+
+run_single(Nbits=700*60, EbNo=6, hf_en=1, wimax_en=1, framesize=4*576, 1)
+%plot_curves
+
+
% Based on Bill Cowley's LDPC simulations\r
%\r
% Octave simulation of BPSK with short LDPC codes developed by Bill. First step\r
-% in use of LDPC codes with FreeDV and Codec 2 700C.\r
+% in use of LDPC codes with FreeDV and Codec 2 700C. BPSK makes life a bits easier\r
+% for simulation\r
%\r
% NOTE: You will need to set the CML path in the call to init_cml() below\r
% for you CML install. See lpdc.m for instructions on how to install \r
\r
1;\r
\r
+gp_interleaver;\r
\r
function init_cml(path_to_cml)\r
currentdir = pwd;\r
\r
function sim_out = run_sim(sim_in, HRA, Ntrials)\r
\r
+ rand('seed',1);\r
+ randn('seed',1);\r
+\r
% Note this is effective Eb/No of payload data bits, sorta thing we\r
% plot on BER versus Eb/No graphs of decoded data. So if we have a\r
% rate 1/2 code, each codeword bit will have Eb/No - 3dB.\r
EbNovec = sim_in.EbNovec;\r
\r
genie_Es = sim_in.genie_Es;\r
- packet_size = sim_in.packet_size;\r
code = sim_in.code;\r
hf_en = sim_in.hf_en;\r
verbose = sim_in.verbose;\r
\r
- if strcmp(code, 'golay')\r
- rate = 0.5;\r
- code_param.data_bits_per_frame = 12;\r
- framesize = 24;\r
+ if isfield(sim_in, "interleave_en") \r
+ interleave_en = sim_in.interleave_en;\r
+ interleave_frames = sim_in.interleave_frames;\r
+ else\r
+ interleave_en = 0;\r
+ interleave_frames = 1;\r
end\r
+\r
+ % default number of carriers for HF model, we fiddle with this a bit\r
+ % for different FEC schemes\r
+\r
+ Nc = 28;\r
+\r
+ % set up for different FEC, get roughly the same frame size so for\r
+ % HF channels the simulation runs over roughly the same time\r
+\r
if strcmp(code, 'ldpc')\r
- [Nr Nc] = size(HRA); \r
- rate = (Nc-Nr)/Nc;\r
- framesize = Nc;\r
+ [numr numc] = size(HRA);\r
+ rate = (numc-numr)/numc;\r
modulation = 'BPSK';\r
demod_type = 0;\r
decoder_type = 0;\r
code_param.H_cols = H_cols;\r
code_param.P_matrix = [];\r
code_param.data_bits_per_frame = length(code_param.H_cols) - length( code_param.P_matrix ); \r
+ code_param.code_bits_per_frame = numc;\r
+ Ncodewordsperframe = interleave_frames*224/numc;\r
+ framesize = Ncodewordsperframe*numc;\r
+ end\r
+\r
+ if strcmp(code, 'golay')\r
+ rate = 0.5;\r
+ code_param.data_bits_per_frame = 12;\r
+ code_param.code_bits_per_frame = 24;\r
+\r
+ % one Golay (24,12) codeword per row\r
+\r
+ Nc = 24; \r
+ Ncodewordsperframe = 9;\r
+ framesize = Nc*Ncodewordsperframe;\r
end\r
+\r
if strcmp(code, 'diversity')\r
rate = 0.5;\r
- code_param.data_bits_per_frame = 7;\r
- framesize = 14;\r
+ code_param.data_bits_per_frame = Nc/2;\r
+ code_param.code_bits_per_frame = Nc;\r
+ Ncodewordsperframe = 8;\r
+ framesize = Nc*Ncodewordsperframe;\r
end\r
\r
+ Rs = 50; Tp = (framesize/Nc)/Rs; Tp_codec2 = 0.04;\r
+\r
+ % we are using BPSK here\r
+\r
mod_order = 2; \r
bps = code_param.bits_per_symbol = log2(mod_order);\r
\r
\r
if hf_en\r
\r
- % some typical values, assume symbols spread acroos Nc OFDM carriers\r
+ % some typical values, assume symbols spread across Nc OFDM carriers\r
\r
- Nc = 28;\r
- Rs = 50; dopplerSpreadHz = 0.5; path_delay = 2E-3*Rs;\r
+ dopplerSpreadHz = 1; path_delay = 1E-3*Rs;\r
\r
- nsymb = Ntrials*framesize;\r
- Tp = (framesize/Nc)/Rs;\r
+ nsymb = Ntrials*framesize*bps;\r
spread1 = doppler_spread(dopplerSpreadHz, Rs, 1.1*ceil(nsymb/Nc));\r
spread2 = doppler_spread(dopplerSpreadHz, Rs, 1.1*ceil(nsymb/Nc));\r
hf_gain = 1.0/sqrt(var(spread1)+var(spread2));\r
% printf("nsymb: %d lspread1: %d\n", nsymb, length(spread1));\r
hf_model = zeros(ceil(nsymb/Nc),Nc);\r
- end\r
+ end\r
+\r
+ % --------------------------------------------------------------------\r
+ % Sun a simulation for every EbNovec point\r
+ %---------------------------------------------------------------------\r
\r
- % loop around each EbNovec point\r
- \r
for ne = 1:length(EbNovec)\r
\r
% Given Eb/No of payload data bits, work out Es/No we need to\r
EbNodB = EbNovec(ne);\r
EsNodB = EbNodB + 10*log10(code_param.bits_per_symbol*rate);\r
EsNo = 10^(EsNodB/10);\r
+\r
+ % reset Hf channel index here so each run gets exactly the same HF channel\r
+\r
hf_r = 1;\r
\r
+ % bunch of counters and logs\r
+\r
Terrs = 0; Tbits = 0; Ferrs = 0; Terrs_raw = Tbits_raw = 0;\r
- tx_bits = rx_bits = [];\r
r_log = []; Nerrs_log = []; Nerrs_raw_log = [];\r
+ tx_bits_log = rx_bits_log = [];\r
\r
for nn = 1: Ntrials \r
- data = round( rand( 1, code_param.data_bits_per_frame ) );\r
- tx_bits = [tx_bits data];\r
- if strcmp(code, 'golay')\r
- codeword = egolayenc(data);\r
- end\r
+\r
+ % Each trial is one complete frame - OK first set up frame to tx\r
+\r
+ codeword = tx_bits = [];\r
+\r
if strcmp(code, 'ldpc')\r
- codeword = LdpcEncode( data, code_param.H_rows, code_param.P_matrix );\r
+ for f=1:Ncodewordsperframe\r
+ data = round( rand( 1, code_param.data_bits_per_frame ) );\r
+ tx_bits = [tx_bits data];\r
+ codeword = [codeword LdpcEncode(data, code_param.H_rows, code_param.P_matrix) ];\r
+ end\r
+ if interleave_en\r
+ codeword = gp_interleave(codeword);\r
+ end\r
end\r
+\r
+ if strcmp(code, 'golay')\r
+ for f=1:Ncodewordsperframe\r
+ data = round( rand( 1, code_param.data_bits_per_frame ) );\r
+ tx_bits = [tx_bits data];\r
+ codeword = [codeword egolayenc(data)];\r
+ end\r
+ end\r
+\r
if strcmp(code, 'diversity')\r
- codeword = [data data];\r
+ for f=1:Ncodewordsperframe\r
+ data = round( rand( 1, code_param.data_bits_per_frame ) );\r
+ tx_bits = [tx_bits data];\r
+ codeword = [codeword data data];\r
+ end\r
end\r
\r
- code_param.code_bits_per_frame = length( codeword );\r
+ tx_bits_log = [tx_bits_log tx_bits];\r
+\r
Nsymb = code_param.code_bits_per_frame/bps; \r
\r
% modulate\r
Tbits_raw += length(codeword);\r
Nerrs_raw_log = [Nerrs_raw_log Nerrs_raw];\r
\r
- % other decoders\r
+ % FEC decoder\r
\r
- if strcmp(code, 'golay')\r
- detected_data = egolaydec(real(r) < 0);\r
- detected_data = detected_data(code_param.data_bits_per_frame+1:framesize);\r
- end\r
+ rx_bits = [];\r
\r
if strcmp(code, 'ldpc')\r
+ if interleave_en\r
+ r = gp_deinterleave(r);\r
+ end\r
\r
% in the binary case the LLRs are just a scaled version of the rx samples ...\r
% if the EsNo is known -- use the following \r
input_decoder_c = 4 * estEsN0 * r;\r
end\r
\r
- [x_hat, PCcnt] = MpDecode( input_decoder_c, code_param.H_rows, code_param.H_cols, ...\r
+ for f=1:Ncodewordsperframe\r
+ st = (f-1)*code_param.code_bits_per_frame + 1;\r
+ en = st + code_param.code_bits_per_frame - 1;\r
+ [x_hat, PCcnt] = MpDecode(input_decoder_c(st:en), code_param.H_rows, code_param.H_cols, ...\r
max_iterations, decoder_type, 1, 1);\r
- Niters = sum(PCcnt!=0);\r
- detected_data = x_hat(Niters,:);\r
+ Niters = sum(PCcnt!=0);\r
+ arx_bits = x_hat(Niters,:);\r
+ rx_bits = [rx_bits arx_bits(1:code_param.data_bits_per_frame)];\r
\r
- if isfield(sim_in, "c_include_file") \r
+ #{\r
+ if isfield(sim_in, "c_include_file") \r
\r
- % optionally dump code and unit test data to a C header file\r
+ % optionally dump code and unit test data to a C header file\r
\r
- code_param.c_include_file = sim_in.c_include_file;\r
- ldpc_gen_h_file(code_param, max_iterations, decoder_type, input_decoder_c, x_hat, detected_data);\r
+ code_param.c_include_file = sim_in.c_include_file;\r
+ ldpc_gen_h_file(code_param, max_iterations, decoder_type, input_decoder_c, x_hat, detected_data);\r
+ #}\r
end\r
\r
detected_data = detected_data(1:code_param.data_bits_per_frame);\r
end\r
\r
+ if strcmp(code, 'golay')\r
+ for f=1:Ncodewordsperframe\r
+ st = (f-1)*code_param.code_bits_per_frame+1; en=st+code_param.code_bits_per_frame-1;\r
+ arx_codeword = egolaydec(real(r(st:en)) < 0);\r
+ rx_bits = [rx_bits arx_codeword(code_param.data_bits_per_frame+1:code_param.code_bits_per_frame)];\r
+ end\r
+ end\r
+\r
if strcmp(code, 'diversity')\r
- detected_data = real(r(1:7) + r(8:14)) < 0;\r
+ for f=1:Ncodewordsperframe\r
+ st = (f-1)*Nc+1;\r
+ r_combined = r(st:st+Nc/2-1) + r(st+Nc/2:st+Nc-1);\r
+ arx_data = real(r_combined) < 0;\r
+ rx_bits = [rx_bits arx_data];\r
+ end\r
+\r
+ #{\r
+ % simulate low rate code to mop up errors\r
+\r
+ error_positions = xor(rx_bits, tx_bits);\r
+ Nerrs = sum(error_positions);\r
+ if Nerrs < 6\r
+ rx_bits = tx_bits;\r
+ end \r
+ #}\r
end\r
\r
- rx_bits = [rx_bits detected_data];\r
+ rx_bits_log = [rx_bits_log rx_bits];\r
\r
- error_positions = xor( detected_data, data );\r
+ error_positions = xor(rx_bits, tx_bits);\r
Nerrs = sum(error_positions);\r
Nerrs_log = [Nerrs_log Nerrs];\r
\r
- if Nerrs>0, Ferrs = Ferrs +1; end\r
+ if Nerrs>0, Ferrs += 1; end\r
Terrs = Terrs + Nerrs;\r
- Tbits = Tbits + code_param.data_bits_per_frame;\r
+ Tbits = Tbits + code_param.data_bits_per_frame*Ncodewordsperframe;\r
end\r
\r
- % count packet errors using supplied packet size. Operating one\r
- % one big string of bits lets us account for cases were FEC framesize\r
- % is less than packet size\r
-\r
- error_positions = xor(tx_bits, rx_bits);\r
- Perrs = 0; Tpackets = 0;\r
- for i=1:packet_size:length(tx_bits)-packet_size\r
- Nerrs = sum(error_positions(i:i+packet_size-1));\r
- if Nerrs>0, Perrs = Perrs +1; end\r
- Tpackets++;\r
+ % Alternative Codec 2 packet rate measurement indep of framesize\r
+\r
+ Nerrs_codec2_log = [];\r
+ Ncodecpacketsize = 28;\r
+ Perrs = 0; Npackets = floor(length(tx_bits_log)/Ncodecpacketsize);\r
+ for p=1:Ncodecpacketsize:Npackets*Ncodecpacketsize\r
+ Nerrs = sum(xor(tx_bits_log(p:p+Ncodecpacketsize-1), rx_bits_log(p:p+Ncodecpacketsize-1)));\r
+ if Nerrs\r
+ Perrs++;\r
+ end\r
+ Nerrs_codec2_log = [Nerrs_codec2_log Nerrs];\r
end\r
\r
printf("Coded EbNo: %3.1f dB BER: %5.4f PER: %5.4f Nbits: %4d Nerrs: %4d Tpackets: %4d Perr: %4d\n", \r
- EbNodB, Terrs/Tbits, Perrs/ Tpackets, Tbits, Terrs, Tpackets, Perrs);\r
+ EbNodB, Terrs/Tbits, Ferrs/ Ntrials, Tbits, Terrs, Ntrials, Ferrs);\r
EbNodB_raw = EsNodB - 10*log10(code_param.bits_per_symbol);\r
printf("Raw EbNo..: %3.1f dB BER: %5.4f Nbits: %4d Nerrs: %4d\n", EbNodB_raw, \r
Terrs_raw/Tbits_raw, Tbits_raw, Terrs_raw);\r
-\r
- TERvec(ne) = Terrs;\r
- FERvec(ne) = Ferrs;\r
- BERvec(ne) = Terrs/ Tbits;\r
- PERvec(ne) = Perrs/ Tpackets;\r
+ printf("Codec 2 PER: %5.4f Npackets: %d Perrs: %d\n", Perrs/Npackets, Npackets, Perrs);\r
+ \r
+ BERvec(ne) = Terrs/Tbits;\r
+ PERvec(ne) = Perrs/Npackets;\r
\r
sim_out.BERvec = BERvec;\r
sim_out.PERvec = PERvec;\r
- sim_out.FERvec = FERvec;\r
- sim_out.TERvec = TERvec;\r
- sim_out.framesize = framesize;\r
+\r
+ error_positions = xor(tx_bits_log, rx_bits_log);\r
sim_out.error_positions = error_positions;\r
\r
if verbose\r
axis([-2 2 -2 2])\r
title('Scatter');\r
\r
- figure(4)\r
-\r
+ figure(4); clf;\r
subplot(211);\r
stem((1:Ntrials)*Tp, Nerrs_raw_log);\r
subplot(212);\r
- stem((1:Ntrials)*Tp, Nerrs_log);\r
- %axis([1 Ntrials*Tp -0.2 1.2])\r
+ stem((1:Npackets)*Tp_codec2, Nerrs_codec2_log);\r
\r
if hf_en\r
- figure(5); clf;\r
+ figure(6); clf;\r
\r
% limit mesh plot to Np points to plot quickly\r
\r
endfunction\r
\r
\r
-function plot_curves(hf_en)\r
+function plot_curves(Ntrials=100, hf_en=0)\r
\r
if hf_en\r
epslabel = 'hf';\r
epslabel = 'awgn';\r
end\r
\r
- Ntrials = 500;\r
-\r
sim_in.genie_Es = 1;\r
- sim_in.packet_size = 28;\r
sim_in.code = 'ldpc';\r
sim_in.hf_en = hf_en;\r
sim_in.verbose = 0;\r
\r
printf("HRA_112_112-------------\n");\r
sim_out1 = run_sim(sim_in, HRA_112_112, Ntrials);\r
+\r
+#{\r
printf("HRA_112_56-------------\n");\r
sim_out2 = run_sim(sim_in, HRA_112_56 , Ntrials);\r
-#{\r
+#}\r
+\r
printf("HRA_56_56-------------\n");\r
sim_out3 = run_sim(sim_in, HRA_56_56 , Ntrials*2);\r
+#{\r
printf("HRA_56_28-------------\n");\r
sim_out4 = run_sim(sim_in, HRA_56_28 , Ntrials*2);\r
- %printf("Golay -------------\n");\r
- %sim_in.code = 'golay';\r
- %sim_out5 = run_sim(sim_in, [], Ntrials*10);\r
- %printf("Diversity -------------\n");\r
- %sim_in.code = 'diversity';\r
- %sim_out6 = run_sim(sim_in, [], Ntrials*10);\r
+\r
+ printf("Golay -------------\n");\r
+ sim_in.code = 'golay';\r
+ sim_out5 = run_sim(sim_in, [], Ntrials);\r
+ \r
#}\r
+ printf("Diversity -------------\n");\r
+ sim_in.code = 'diversity';\r
+ sim_out6 = run_sim(sim_in, [], Ntrials);\r
\r
if hf_en\r
Ebvec_theory = 2:0.5:12;\r
% no packet error if all bits ok (1-p(0))*(1-p(1))\r
% P(packet error) = p(0)+p(1)+....\r
\r
- uncoded_PER_theory = 1 - (1-uncoded_BER_theory).^sim_in.packet_size;\r
+ uncoded_PER_theory = 1 - (1-uncoded_BER_theory).^28;\r
\r
figure(1); clf;\r
semilogy(Ebvec_theory, uncoded_BER_theory, 'b+-;BPSK theory;','markersize', 10, 'linewidth', 2)\r
hold on;\r
semilogy(EbNovec, sim_out1.BERvec, 'g+-;rate 1/2 HRA 112 112;','markersize', 10, 'linewidth', 2)\r
- semilogy(EbNovec, sim_out2.BERvec, 'r+-;rate 2/3 HRA 112 56;','markersize', 10, 'linewidth', 2)\r
- %semilogy(EbNovec, sim_out3.BERvec, 'c+-;rate 1/2 HRA 56 56;','markersize', 10, 'linewidth', 2)\r
+ %semilogy(EbNovec, sim_out2.BERvec, 'r+-;rate 2/3 HRA 112 56;','markersize', 10, 'linewidth', 2)\r
+ semilogy(EbNovec, sim_out3.BERvec, 'c+-;rate 1/2 HRA 56 56;','markersize', 10, 'linewidth', 2)\r
%semilogy(EbNovec, sim_out4.BERvec, 'k+-;rate 2/3 HRA 56 28;','markersize', 10, 'linewidth', 2)\r
%semilogy(EbNovec, sim_out5.BERvec, 'm+-;rate 1/2 Golay (24,12);','markersize', 10, 'linewidth', 2)\r
- %semilogy(EbNovec, sim_out6.BERvec, 'go-;rate 1/2 Diversity;','markersize', 10, 'linewidth', 2)\r
+ semilogy(EbNovec, sim_out6.BERvec, 'go-;rate 1/2 Diversity;','markersize', 10, 'linewidth', 2)\r
hold off;\r
xlabel('Eb/No')\r
ylabel('BER')\r
semilogy(Ebvec_theory, uncoded_PER_theory, 'b+-;BPSK theory;','markersize', 10, 'linewidth', 2)\r
hold on;\r
semilogy(EbNovec, sim_out1.PERvec, 'g+-;rate 1/2 HRA 112 112;','markersize', 10, 'linewidth', 2)\r
- semilogy(EbNovec, sim_out2.PERvec, 'r+-;rate 2/3 HRA 112 56;','markersize', 10, 'linewidth', 2)\r
- %semilogy(EbNovec, sim_out3.PERvec, 'c+-;rate 1/2 HRA 56 56;','markersize', 10, 'linewidth', 2)\r
+ %semilogy(EbNovec, sim_out2.PERvec, 'r+-;rate 2/3 HRA 112 56;','markersize', 10, 'linewidth', 2)\r
+ semilogy(EbNovec, sim_out3.PERvec, 'c+-;rate 1/2 HRA 56 56;','markersize', 10, 'linewidth', 2)\r
%semilogy(EbNovec, sim_out4.PERvec, 'k+-;rate 2/3 HRA 56 28;','markersize', 10, 'linewidth', 2)\r
%semilogy(EbNovec, sim_out5.PERvec, 'm+-;rate 1/2 Golay (24,12);','markersize', 10, 'linewidth', 2)\r
- %semilogy(EbNovec, sim_out6.PERvec, 'go-;rate 1/2 Diversity;','markersize', 10, 'linewidth', 2)\r
+ semilogy(EbNovec, sim_out6.PERvec, 'go-;rate 1/2 Diversity;','markersize', 10, 'linewidth', 2)\r
hold off;\r
xlabel('Eb/No')\r
ylabel('PER')\r
endfunction\r
\r
\r
-function run_single(bits, code = 'ldpc', channel = 'awgn', EbNodB, error_pattern_filename)\r
+function run_single(bits, code = 'ldpc', channel = 'awgn', EbNodB, interleave=0, error_pattern_filename)\r
\r
sim_in.code = code;\r
load HRA_112_112.txt\r
+ HRA = HRA_112_112;\r
\r
- if strcmp(code, 'ldpc')\r
- data_bits_per_frame = 112;\r
- rate = 0.5;\r
- end\r
-\r
- if strcmp(code, 'diversity')\r
- data_bits_per_frame = 7;\r
- rate = 0.5;\r
- end\r
-\r
- Ntrials = bits/data_bits_per_frame;\r
- sim_in.genie_Es = 1;\r
- sim_in.packet_size = 28;\r
+ sim_in.genie_Es = 1;\r
\r
+ sim_in.EbNovec = EbNodB;\r
if strcmp(channel, 'awgn')\r
sim_in.hf_en = 0;\r
- sim_in.EbNovec = EbNodB;\r
else\r
sim_in.hf_en = 1;\r
- sim_in.EbNovec = EbNodB;\r
end\r
\r
+ if interleave\r
+ sim_in.interleave_en = 1;\r
+ sim_in.interleave_frames = interleave;\r
+ else\r
+ sim_in.interleave_frames = 1;\r
+ end\r
+ Ntrials = floor(bits/(112*sim_in.interleave_frames));\r
sim_in.verbose = 1;\r
- sim_out = run_sim(sim_in, HRA_112_112, Ntrials);\r
+ sim_out = run_sim(sim_in, HRA, Ntrials);\r
\r
- if nargin == 5\r
+ if nargin == 6\r
fep = fopen(error_pattern_filename, "wb");\r
fwrite(fep, sim_out.error_positions, "short");\r
fclose(fep);\r
bits = data_bits_per_frame;\r
Ntrials = bits/data_bits_per_frame;\r
sim_in.genie_Es = 1;\r
- sim_in.packet_size = 28;\r
sim_in.hf_en = 0;\r
sim_in.Esvec = 2;\r
sim_in.c_include_file = "../src/HRA_112_112.h";\r
\r
init_cml('/home/david/Desktop/cml');\r
\r
-rand('seed',1);\r
-randn('seed',1);\r
more off;\r
format;\r
\r
% simple single point test\r
\r
-run_single(700*120, 'ldpc', 'hf', 12)\r
+run_single(700*60, 'ldpc', 'hf', 6, 32)\r
\r
% plotting curves (may take a while)\r
\r
%plot_curves(0);\r
-%plot_curves(1);\r
+%plot_curves(500,1);\r
\r
% generating error files \r
\r
-#{\r
-run_single(700*10, 'ldpc', 'awgn', 2, 'awgn_2dB_ldpc.err')\r
-run_single(700*10, 'diversity', 'awgn', 2, 'awgn_2dB_diversity.err')\r
-run_single(700*10, 'ldpc', 'hf', 6, 'hf_6dB_ldpc.err')\r
-run_single(700*10, 'diversity', 'hf', 6, 'hf_6dB_diversity.err')\r
-#}\r
+%run_single(700*10, 'ldpc', 'awgn', 3, 0, 'awgn_3dB_ldpc.err')\r
+%run_single(700*10, 'diversity', 'awgn', 3, 0, 'awgn_3dB_diversity.err')\r
+%run_single(700*10, 'ldpc', 'hf', 10, 0, 'hf_10dB_ldpc.err')\r
+%run_single(700*10, 'diversity', 'hf', 10, 0, 'hf_10dB_diversity.err')\r
\r
% generate C header for C port of code\r
\r
projects / freetel-svn-tracking.git / commitdiff