code_param.data_bits_per_frame = length(code_param.H_cols) - length( code_param.P_matrix );\r
code_param.S_matrix = CreateConstellation( modulation, mod_order, mapping );\r
code_param.bits_per_symbol = log2(mod_order);\r
+ code_param.code_bits_per_frame = framesize;\r
+ code_param.symbols_per_frame = framesize/2;\r
endfunction\r
\r
\r
-% Gray coded QPSK modulation function\r
-\r
-function symbol = qpsk_mod(two_bits)\r
- two_bits_decimal = sum(two_bits .* [2 1]); \r
- switch(two_bits_decimal)\r
- case (0) symbol = 1;\r
- case (1) symbol = j;\r
- case (2) symbol = -j;\r
- case (3) symbol = -1;\r
- endswitch\r
-endfunction\r
-\r
-\r
-% Gray coded QPSK demodulation function\r
-\r
-function two_bits = qpsk_demod(symbol)\r
- bit0 = real(symbol*exp(j*pi/4)) < 0;\r
- bit1 = imag(symbol*exp(j*pi/4)) < 0;\r
- two_bits = [bit1 bit0];\r
-endfunction\r
-\r
-\r
-% inserts a unique word into a frame of bits. The UW bits are spread\r
-% throughout the input frame 2 bits at a time.\r
-\r
-function frameout = insert_uw(framein, uw)\r
-\r
- luw = length(uw);\r
- lframein = length(framein);\r
- spacing = 2*lframein/luw;\r
-\r
- frameout = [];\r
-\r
- pin = 1; pout = 1; puw = 1;\r
- while (luw)\r
- %printf("pin %d pout %d puw %d luw %d\n", pin, pout, puw, luw);\r
- frameout(pout:pout+spacing-1) = framein(pin:pin+spacing-1);\r
- pin += spacing; \r
- pout += spacing;\r
- frameout(pout:pout+1) = uw(puw:puw+1);\r
- puw += 2;\r
- pout += 2;\r
- luw -= 2;\r
- end\r
-endfunction\r
-\r
-% removes a unique word from a frame of bits. The UW bits are spread\r
-% throughout the input frame 2 bits at a time.\r
-\r
-function frameout = remove_uw(framein, lvd, luw)\r
-\r
- spacing = 2*lvd/luw;\r
-\r
- frameout = [];\r
-\r
- pin = 1; pout = 1;\r
- while (luw)\r
- %printf("pin %d pout %d luw %d ", pin, pout, luw);\r
- %printf("pin+spacing-1 %d lvd %d lframein: %d\n", pin+spacing-1, lvd, length(framein));\r
- frameout(pout:pout+spacing-1) = framein(pin:pin+spacing-1);\r
- pin += spacing + 2; \r
- pout += spacing;\r
- luw -= 2;\r
- end\r
-\r
-endfunction\r
-\r
-\r
-% removes a unique word from a frame of symbols. The UW symbols are spread\r
-% throughout the input frame 1 symbol at a time.\r
-\r
-function framesymbolsout = remove_uw_symbols(framesymbolsin, ldatasymbols, luwsymbols)\r
-\r
- spacing = ldatasymbols/luwsymbols;\r
-\r
- framesymbolsout = [];\r
-\r
- pin = 1; pout = 1;\r
- while (luwsymbols)\r
- %printf("pin %d pout %d luw %d ", pin, pout, luwsymbols);\r
- %printf("pin+spacing-1 %d ldatasymbols %d lframein: %d\n", pin+spacing-1, ldatasymbols, length(framesymbolsin));\r
- framesymbolsout(pout:pout+spacing-1) = framesymbolsin(pin:pin+spacing-1);\r
- pin += spacing + 1; \r
- pout += spacing;\r
- luwsymbols--;\r
- end\r
-\r
-endfunction\r
-\r
-\r
-\r
-% builds up a sparse QPSK modulated version version of the UW for use\r
-% in UW sync at the rx\r
-\r
-function mod_uw = build_mod_uw(uw, spacing)\r
- luw = length(uw);\r
-\r
- mod_uw = [];\r
-\r
- pout = 1; puw = 1;\r
- while (luw)\r
- %printf("pin %d pout %d puw %d luw %d\n", pin, pout, puw, luw);\r
- pout += spacing/2;\r
- mod_uw(pout) = qpsk_mod(uw(puw:puw+1));\r
- puw += 2;\r
- pout += 1;\r
- luw -= 2;\r
- end\r
-endfunction\r
-\r
-\r
-% Uses the UW to determine when we have a full codeword ready for decoding\r
-\r
-function [found_uw corr] = look_for_uw(mem_rx_symbols, mod_uw)\r
- sparse_mem_rx_symbols = mem_rx_symbols(find(mod_uw));\r
-\r
- % correlate with ref UW\r
-\r
- num = (mem_rx_symbols * mod_uw') .^ 2;\r
- den = (sparse_mem_rx_symbols * sparse_mem_rx_symbols') * (mod_uw * mod_uw');\r
- \r
- corr = abs(num/(den+1E-6));\r
- found_uw = corr > 0.8;\r
-endfunction\r
-\r
\r
function [codeword s] = ldpc_enc(data, code_param)\r
codeword = LdpcEncode( data, code_param.H_rows, code_param.P_matrix );\r
\r
\r
function detected_data = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, r, EsNo, fading)\r
+ if nargin == 6\r
+ fading = ones(1, length(r));\r
+ end\r
+\r
symbol_likelihood = Demod2D( r, code_param.S_matrix, EsNo, fading);\r
\r
% initialize the extrinsic decoder input\r
+\r
input_somap_c = zeros(1, code_param.code_bits_per_frame );\r
bit_likelihood = Somap( symbol_likelihood, demod_type, input_somap_c );\r
\r
% Our LDPC library\r
\r
ldpc;\r
+qpsk_;\r
\r
function sim_out = run_simulation(sim_in)\r
\r
clf\r
semilogy(EbNodB, uncoded_BER_theory,'r;uncoded QPSK theory;')\r
hold on;\r
-semilogy(EbNodB-10*log10(sim_in.rate), sim_out.BER,'g;LDPC coded QPSK simulation;');\r
+semilogy(EbNodB-10*log10(sim_in.rate), sim_out.BER+1E-10,'g;LDPC coded QPSK simulation;');\r
hold off;\r
grid('minor')\r
xlabel('Eb/No (dB)')\r
ylabel('BER')\r
+axis([min(EbNodB) max(EbNodB) min(uncoded_BER_theory) 1])\r
--- /dev/null
+% qpsk.m
+%
+% David Rowe Sep 2015
+%
+% Octave functions to implement a QPSK modem
+
+1;
+
+% Gray coded QPSK modulation function
+
+function symbol = qpsk_mod(two_bits)
+ two_bits_decimal = sum(two_bits .* [2 1]);
+ switch(two_bits_decimal)
+ case (0) symbol = 1;
+ case (1) symbol = j;
+ case (2) symbol = -j;
+ case (3) symbol = -1;
+ endswitch
+endfunction
+
+
+% Gray coded QPSK demodulation function
+
+function two_bits = qpsk_demod(symbol)
+ bit0 = real(symbol*exp(j*pi/4)) < 0;
+ bit1 = imag(symbol*exp(j*pi/4)) < 0;
+ two_bits = [bit1 bit0];
+endfunction
+
+
+% Inserts pilot symbols a frame of symbols. The pilot symbols are
+% spread evenly throughout the input frame.
+
+function frameout = insert_pilots(framein, pilots, Npilotstep)
+
+ lpilots = length(pilots);
+ lframein = length(framein);
+ frameout = zeros(1, lframein + lpilots);
+
+ pin = 1; pout = 1; ppilots = 1;
+ while (lpilots)
+ %printf("pin %d pout %d ppilots %d lpilots %d\n", pin, pout, ppilots, lpilots);
+ frameout(pout:pout+Npilotstep-1) = framein(pin:pin+Npilotstep-1);
+ pin += Npilotstep;
+ pout += Npilotstep;
+ frameout(pout:pout) = pilots(ppilots);
+ ppilots++;
+ pout++;
+ lpilots--;
+ end
+endfunction
+
+
+% Removes the pilots symbols from a frame of symbols.
+
+function frameout = remove_pilots(framein, pilots, Npilotstep)
+
+ frameout = [];
+ lpilots = length(pilots);
+
+ pin = 1; pout = 1;
+ while (lpilots)
+ %printf("pin %d pout %d lpilots %d ", pin, pout, lpilots);
+ %printf("pin+spacing-1 %d lvd %d lframein: %d\n", pin+spacing-1, lvd, length(framein));
+ frameout(pout:pout+Npilotstep-1) = framein(pin:pin+Npilotstep-1);
+ pin += Npilotstep+1;
+ pout += Npilotstep;
+ lpilots--;
+ end
+
+endfunction
+
+
+% Estimate and correct phase offset using a window of Np pilots around
+% current symbol
+
+function symbpilot_rx = correct_phase_offset(aqpsk, symbpilot_rx)
+ rx_pilot_buf = aqpsk.rx_pilot_buf;
+ Npilotstep = aqpsk.Npilotstep;
+ Nsymb = aqpsk.Nsymb;
+
+ for ns=1:Npilotstep+1:Nsymb
+
+ % update buffer of recent pilots, note we need past ones
+
+ rx_pilot_buf(1) = rx_pilot_buf(2);
+ next_pilot_index = ceil(ns/(Npilotstep+1))*(Npilotstep+1);
+ rx_pilot_buf(2) = symbpilot_rx(next_pilot_index);
+
+ % average pilot symbols to get estimate of phase
+
+ phase_est = angle(sum(rx_pilot_buf));
+
+ %printf("next_pilot_index: %d phase_est: %f\n", next_pilot_index, phase_est);
+
+ % now correct the phase of each symbol
+
+ for s=ns:ns+Npilotstep
+ symbpilot_rx(s) *= exp(-j*phase_est);
+ end
+ end
+
+ aqpsk.rx_pilot_buf = rx_pilot_buf;
+endfunction
+
+
+% builds up a sparse QPSK modulated version version of the UW for use
+% in UW sync at the rx
+
+function mod_uw = build_mod_uw(uw, spacing)
+ luw = length(uw);
+
+ mod_uw = [];
+
+ pout = 1; puw = 1;
+ while (luw)
+ %printf("pin %d pout %d puw %d luw %d\n", pin, pout, puw, luw);
+ pout += spacing/2;
+ mod_uw(pout) = qpsk_mod(uw(puw:puw+1));
+ puw += 2;
+ pout += 1;
+ luw -= 2;
+ end
+endfunction
+
+
+% Uses the UW to determine when we have a full codeword ready for decoding
+
+function [found_uw corr] = look_for_uw(mem_rx_symbols, mod_uw)
+ sparse_mem_rx_symbols = mem_rx_symbols(find(mod_uw));
+
+ % correlate with ref UW
+
+ num = (mem_rx_symbols * mod_uw') .^ 2;
+ den = (sparse_mem_rx_symbols * sparse_mem_rx_symbols') * (mod_uw * mod_uw');
+
+ corr = abs(num/(den+1E-6));
+ found_uw = corr > 0.8;
+endfunction
+
--- /dev/null
+% tqpsk.m
+%
+% David Rowe Sep 2015
+%
+% QPSK modem Octave test sctipt
+
+% make sure we get same random numbers every time
+
+more off;
+qpsk;
+ldpc;
+
+function sim_out = run_simulation(sim_in)
+ Ntrials = sim_in.Ntrials;
+ EsNodBvec = sim_in.EsNodBvec;
+ verbose = sim_in.verbose;
+ phase_offset = sim_in.phase_offset;
+ phase_est_en = sim_in.phase_est_en;
+ verbose = sim_in.verbose;
+
+ rand('state',1);
+ randn('state',1);
+
+ % Init modem ------------------------------------------------------------------------------
+
+ aqpsk.Rbdata = 2400; % bit rate (Hz) of payload data
+ aqpsk.code_rate = 0.5;
+ aqpsk.framesize = 1152; % number of bits in LDPC encoded frame (data plus parity bits)
+ aqpsk.Ndatabits = aqpsk.framesize*aqpsk.code_rate;
+ aqpsk.Rbcoded = aqpsk.Rbdata/aqpsk.code_rate; % bit rate (Hz) of LPDC coded data
+ aqpsk.Rf = aqpsk.Rbcoded/aqpsk.framesize; % frame rate (Hz)
+ aqpsk.bits_per_symb = 2; % QPSK has two bits per symbol
+ aqpsk.Rscoded = aqpsk.Rbcoded/aqpsk.bits_per_symb; % symbol rate of LDPC coded data
+ aqpsk.Nsymb = aqpsk.framesize/aqpsk.bits_per_symb; % number of coded data symbols/frame
+ aqpsk.Npilotstep = 8; % number of data symbols between each pilot symbol
+
+ aqpsk.Npilotsframe = (aqpsk.framesize/aqpsk.bits_per_symb)/aqpsk.Npilotstep;
+ if aqpsk.Npilotsframe != floor(aqpsk.Npilotsframe)
+ error("Npilotsframe must be an integer");
+ end
+
+ aqpsk.Nsymbpilot = aqpsk.Nsymb + aqpsk.Npilotsframe; % total number of symbols per frame including pilots
+ aqpsk.Rs = aqpsk.Rf*aqpsk.Nsymbpilot; % total symbol rate (Hz)
+ aqpsk.pilots = ones(1, aqpsk.Npilotsframe);
+ aqpsk.Np = 2; % number of pilots used for phase estimation
+ aqpsk.rx_pilot_buf = zeros(1,aqpsk.Np);
+
+ Nsymb = aqpsk.Nsymb;
+ Nsymbpilot = aqpsk.Nsymbpilot;
+
+ % Init LDPC --------------------------------------------------------------------------
+
+ currentdir = pwd;
+ addpath '/home/david/tmp/cml/mat' % assume the source files stored here
+ cd /home/david/tmp/cml
+ CmlStartup % note that this is not in the cml path!
+ cd(currentdir)
+
+ mod_order = 4;
+ modulation = 'QPSK';
+ mapping = 'gray';
+ demod_type = 0;
+ decoder_type = 0;
+ max_iterations = 100;
+ EsNo = 10; % decoder needs an estimated channel EsNo (linear ratio, not dB)
+
+ code_param = ldpc_init(aqpsk.code_rate, aqpsk.framesize, modulation, mod_order, mapping);
+ assert(code_param.data_bits_per_frame == aqpsk.Ndatabits);
+
+ % Start simulation -------------------------------------------------------------------
+
+ for ne = 1:length(EsNodBvec)
+ EsNodB = EsNodBvec(ne);
+ EsNo = 10^(EsNodB/10);
+ variance = 1/EsNo;
+
+ Tbits = Terrs_uc = Terrs = Ferrs = 0;
+
+ % Encode a bunch of frames
+
+ for nn = 1: Ntrials
+ tx_bits = round(rand(1, aqpsk.Ndatabits));
+ [tx_codeword, tx_coded_symb] = ldpc_enc(tx_bits, code_param);
+ symbpilot_tx = insert_pilots(tx_coded_symb, aqpsk.pilots, aqpsk.Npilotstep);
+
+ noise = sqrt(variance*0.5)*(randn(1,aqpsk.Nsymbpilot) + j*randn(1,aqpsk.Nsymbpilot));
+ symbpilot_rx = symbpilot_tx * exp(j*phase_offset) + noise;
+
+ if phase_est_en
+ symbpilot_rx = correct_phase_offset(aqpsk, symbpilot_rx);
+ end
+ rx_coded_symb = remove_pilots(symbpilot_rx, aqpsk.pilots, aqpsk.Npilotstep);
+ rx_codeword = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, rx_coded_symb, EsNo);
+
+ % measure coded and uncoded errors
+
+ errors_positions = xor(tx_bits, rx_codeword(1:aqpsk.Ndatabits));
+ Nerrs = sum(errors_positions);
+
+ rx_bits_uncoded = zeros(1, aqpsk.Ndatabits);
+ for s=1:Nsymb*aqpsk.code_rate;
+ rx_bits_uncoded(2*(s-1)+1:2*s) = qpsk_demod(rx_coded_symb(s));
+ end
+ errors_positions = xor(tx_bits, rx_bits_uncoded);
+ Nerrs_uc = sum(errors_positions);
+
+ 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 = Terrs + Nerrs;
+ Terrs_uc = Terrs_uc + Nerrs_uc;
+ Tbits = Tbits + aqpsk.Ndatabits;
+ end
+
+ if verbose
+ printf("EsNodB: %3.2f BER: %f Tbits: %d Terrs: %d Ferrs: %d Terrs_uc: %d\n", EsNodB, Terrs/Tbits, Tbits, Terrs, Ferrs, Terrs_uc);
+ end
+
+ sim_out.Tbits(ne) = Tbits;
+ sim_out.Terrs(ne) = Terrs;
+ sim_out.Terrs_uc(ne) = Terrs_uc;
+ sim_out.Ferrs(ne) = Ferrs;
+ sim_out.BER(ne) = Terrs/Tbits;
+ sim_out.FER(ne) = Ferrs/Ntrials;
+ sim_out.BER_uc(ne) = Terrs_uc/Tbits;
+ end
+
+ sim_out.code_rate = aqpsk.code_rate;
+end
+
+
+% ----------------------------------------------------------------------
+
+sim_in.phase_est_en = 1;
+sim_in.phase_offset = 0;
+sim_in.verbose = 1;
+sim_in.EsNodBvec = -2:10;
+sim_in.Ntrials = 10;
+
+sim_out = run_simulation(sim_in);
+
+EbNodB = sim_in.EsNodBvec - 10*log10(2); % QPSK has two bits/symbol
+uncoded_BER_theory = 0.5*erfc(sqrt(10.^(EbNodB/10)));
+
+figure(1)
+clf
+semilogy(EbNodB, uncoded_BER_theory,'r;uncoded QPSK theory;')
+hold on;
+semilogy(EbNodB, sim_out.BER_uc,'r+;uncoded QPSK actual;')
+semilogy(EbNodB-10*log10(sim_out.code_rate), sim_out.BER+1E-10,'g;LDPC coded QPSK simulation;');
+hold off;
+grid('minor')
+xlabel('Eb/No (dB)')
+ylabel('BER')
+axis([min(EbNodB) max(EbNodB) min(uncoded_BER_theory) 1])
+
+% add noise
+% correct freq offset
+% output "demodulated" tx frame to a file of bits for 3rd party modulator
+% actually estimate Es/No
+% meas impl loss for pilot symbols
+% put in Es/No loop, measure coded and uncoded errors
+% separate into end and dec
+
projects / freetel-svn-tracking.git / commitdiff