From: drowe67 Date: Fri, 11 Sep 2015 07:21:12 +0000 (+0000) Subject: building up symbol rate QPSK simulation, pilot insertion/removal, first pass phase... X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=98537294e9af9eef1b6c3771d9c78799fb61cf85;p=freetel-svn-tracking.git building up symbol rate QPSK simulation, pilot insertion/removal, first pass phase est, test framework for drawing curves git-svn-id: https://svn.code.sf.net/p/freetel/code@2306 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/octave/ldpc.m b/codec2-dev/octave/ldpc.m index 9a10d93e..00499625 100644 --- a/codec2-dev/octave/ldpc.m +++ b/codec2-dev/octave/ldpc.m @@ -22,134 +22,11 @@ function code_param = ldpc_init(rate, framesize, modulation, mod_order, mapping) code_param.data_bits_per_frame = length(code_param.H_cols) - length( code_param.P_matrix ); code_param.S_matrix = CreateConstellation( modulation, mod_order, mapping ); code_param.bits_per_symbol = log2(mod_order); + code_param.code_bits_per_frame = framesize; + code_param.symbols_per_frame = framesize/2; endfunction -% 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 a unique word into a frame of bits. The UW bits are spread -% throughout the input frame 2 bits at a time. - -function frameout = insert_uw(framein, uw) - - luw = length(uw); - lframein = length(framein); - spacing = 2*lframein/luw; - - frameout = []; - - pin = 1; pout = 1; puw = 1; - while (luw) - %printf("pin %d pout %d puw %d luw %d\n", pin, pout, puw, luw); - frameout(pout:pout+spacing-1) = framein(pin:pin+spacing-1); - pin += spacing; - pout += spacing; - frameout(pout:pout+1) = uw(puw:puw+1); - puw += 2; - pout += 2; - luw -= 2; - end -endfunction - -% removes a unique word from a frame of bits. The UW bits are spread -% throughout the input frame 2 bits at a time. - -function frameout = remove_uw(framein, lvd, luw) - - spacing = 2*lvd/luw; - - frameout = []; - - pin = 1; pout = 1; - while (luw) - %printf("pin %d pout %d luw %d ", pin, pout, luw); - %printf("pin+spacing-1 %d lvd %d lframein: %d\n", pin+spacing-1, lvd, length(framein)); - frameout(pout:pout+spacing-1) = framein(pin:pin+spacing-1); - pin += spacing + 2; - pout += spacing; - luw -= 2; - end - -endfunction - - -% removes a unique word from a frame of symbols. The UW symbols are spread -% throughout the input frame 1 symbol at a time. - -function framesymbolsout = remove_uw_symbols(framesymbolsin, ldatasymbols, luwsymbols) - - spacing = ldatasymbols/luwsymbols; - - framesymbolsout = []; - - pin = 1; pout = 1; - while (luwsymbols) - %printf("pin %d pout %d luw %d ", pin, pout, luwsymbols); - %printf("pin+spacing-1 %d ldatasymbols %d lframein: %d\n", pin+spacing-1, ldatasymbols, length(framesymbolsin)); - framesymbolsout(pout:pout+spacing-1) = framesymbolsin(pin:pin+spacing-1); - pin += spacing + 1; - pout += spacing; - luwsymbols--; - end - -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 - function [codeword s] = ldpc_enc(data, code_param) codeword = LdpcEncode( data, code_param.H_rows, code_param.P_matrix ); @@ -158,9 +35,14 @@ endfunction function detected_data = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, r, EsNo, fading) + if nargin == 6 + fading = ones(1, length(r)); + end + symbol_likelihood = Demod2D( r, code_param.S_matrix, EsNo, fading); % initialize the extrinsic decoder input + input_somap_c = zeros(1, code_param.code_bits_per_frame ); bit_likelihood = Somap( symbol_likelihood, demod_type, input_somap_c ); diff --git a/codec2-dev/octave/ldpcut.m b/codec2-dev/octave/ldpcut.m index f6d787ef..0b45ff79 100644 --- a/codec2-dev/octave/ldpcut.m +++ b/codec2-dev/octave/ldpcut.m @@ -16,6 +16,7 @@ cd(currentdir) % Our LDPC library ldpc; +qpsk_; function sim_out = run_simulation(sim_in) @@ -150,8 +151,9 @@ figure(1) clf semilogy(EbNodB, uncoded_BER_theory,'r;uncoded QPSK theory;') hold on; -semilogy(EbNodB-10*log10(sim_in.rate), sim_out.BER,'g;LDPC coded QPSK simulation;'); +semilogy(EbNodB-10*log10(sim_in.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]) diff --git a/codec2-dev/octave/qpsk.m b/codec2-dev/octave/qpsk.m new file mode 100644 index 00000000..b502e164 --- /dev/null +++ b/codec2-dev/octave/qpsk.m @@ -0,0 +1,140 @@ +% 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 + diff --git a/codec2-dev/octave/tqpsk.m b/codec2-dev/octave/tqpsk.m new file mode 100644 index 00000000..099a454e --- /dev/null +++ b/codec2-dev/octave/tqpsk.m @@ -0,0 +1,169 @@ +% 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 +