From b69f1336ff0f5c40d35306c779a36955378680e7 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Thu, 23 Mar 2017 06:47:58 +0000 Subject: [PATCH] working on handling phase est ambiguity neat +/- pi/2 git-svn-id: https://svn.code.sf.net/p/freetel/code@3081 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/bpsk.m | 180 ++++++++++++++++++++++++++++++--------- 1 file changed, 138 insertions(+), 42 deletions(-) diff --git a/codec2-dev/octave/bpsk.m b/codec2-dev/octave/bpsk.m index 2a43d225..91dfbc11 100644 --- a/codec2-dev/octave/bpsk.m +++ b/codec2-dev/octave/bpsk.m @@ -5,6 +5,54 @@ 1; +% Differential BPSK detector (m=2) +% ML detection which gains about 0.5dB (m=3, m=4) +% +% Based on JPL publication 89-38 "Multiple Symbol Differential +% Detection of Uncoded and Trellis Coded MPSK" by Divsalar, Simon, +% Shahshahani. Thanks Johhn Gibbs NN7F for advice. + +function rx_symb = dbpsk_demod(m, r) + tx_set = [1 -1]; + + if m == 2 + + % regular DBPSK detection + + rx_symb(i) = r(1) * r(2)'/(abs(r(2))); + + else + + % ML DBPSK detection + + max_eta = 0; rx_symb(i) = tx_set(1); + + for k=1:2 + for k_1=1:2 + for k_2=1:2 + + if m == 3 + eta = abs(r(3) + r(1)*tx_set(k)'*tx_set(k_1)' + r(2)*tx_set(k_1)')^2; + end + + if m == 4 + eta = abs(r(4) + r(1)*tx_set(k)'*tx_set(k_1)'*tx_set(k_2)' + r(2)*tx_set(k_1)'*tx_set(k_2)' + r(3)*tx_set(k_2)')^2; + end + + %printf(" %d %d %f \n", k_1, k, eta); + + if eta > max_eta + max_eta = eta; + rx_symb(i) = tx_set(k); + end + + end + end + end + end +endfunction + + function sim_out = run_sim(sim_in) Rs = 50; @@ -14,21 +62,31 @@ function sim_out = run_sim(sim_in) verbose = sim_in.verbose; m = sim_in.m; phase_est_mem = sim_in.phase_est_mem; - phase_offset = sim_in.phase_offset; + assert(mod(phase_est_mem,2) == 1, "phase_est_mem must be odd"); + phase_est_delay = floor(phase_est_mem/2); + printf("phase_est_mem: %d phase_est_delay: %d\n", phase_est_mem, phase_est_delay); + woffset = 2*pi*sim_in.freq_offset_hz/Rs; - tx = zeros(1,Nbits); - rx = zeros(1,Nbits); for nn=1:length(EbNodB) EbNo = 10 .^ (EbNodB(nn)/10); variance = 1/(EbNo/2); noise = sqrt(variance)*(0.5*randn(1,Nbits) + j*0.5*randn(1,Nbits)); + tx = zeros(1,Nbits); + rx = zeros(1,Nbits); tx_bits = rand(1,Nbits) > 0.5; prev_tx = 1; prev_rx = 1; r = ones(1,max(4,phase_est_mem)); + phase_amb = 0; + phase_offset = sim_in.phase_offset + + phase_offset_log = zeros(1,Nbits); + phase_amb_log = zeros(1,Nbits); + phase_est_stripped = zeros(1,Nbits); + phase_est = zeros(1,Nbits); Nerrs = 0; for i=1:Nbits @@ -40,61 +98,99 @@ function sim_out = run_sim(sim_in) end rx(i) = tx(i)*exp(j*phase_offset) + noise(i); + phase_offset_log(i) = phase_offset; phase_offset += woffset; + if (phase_offset > pi) + phase_offset -= 2*pi; + end r(2:phase_est_mem) = r(1:phase_est_mem-1); r(1) = rx(i); if dbpsk - if m == 2 - %rx_symb(i) = rx(i) * conj(prev_rx)/(abs(prev_rx)); - rx_symb(i) = r(1) * r(2)'/(abs(prev_rx)); - prev_rx = rx(i); - else - tx_set = [1 -1]; - max_eta = 0; rx_symb(i) = tx_set(1); - for k=1:2 - for k_1=1:2 - for k_2=1:2 - if m == 3 - eta = abs(r(3) + r(1)*tx_set(k)'*tx_set(k_1)' + r(2)*tx_set(k_1)')^2; - end - if m == 4 - eta = abs(r(4) + r(1)*tx_set(k)'*tx_set(k_1)'*tx_set(k_2)' + r(2)*tx_set(k_1)'*tx_set(k_2)' + r(3)*tx_set(k_2)')^2; - end - %printf(" %d %d %f \n", k_1, k, eta); - if eta > max_eta - max_eta = eta; - rx_symb(i) = tx_set(k); - end - end - end - end - end + rx_symb(i) = dbpsk_demod(m, r); else - rx_symb(i) = rx(i); if phase_est_mem - % modulation strip - stripped(i) = angle(sum(r(1:phase_est_mem) .^ 2)); - rx_symb(i) *= exp(-j*stripped(i)/2); + + if i >= phase_est_mem + + % demod symbol at centre of phase est window + + centre = i - phase_est_delay; + + % modulation strip + + stripped = r(1:phase_est_mem) .^ 2; + + phase_est_stripped(centre) = angle(sum(stripped))/2; + + % determine if phase has jumped from - -> + + + if (phase_est_stripped(centre-1) < -pi/4) && (phase_est_stripped(centre) > pi/4) + %printf("- -> +\n"); + phase_amb -= pi; + if (phase_amb < -pi) + phase_amb += 2*pi; + end + end + + % determine if phase has jumped from + -> - + + if (phase_est_stripped(centre-1) > pi/4) && (phase_est_stripped(centre) < -pi/4) + %printf("+ -> -\n"); + phase_amb += pi; + if (phase_amb > pi) + phase_amb -= 2*pi; + end + end + + phase_amb_log(centre) = phase_amb; + phase_est(centre) = phase_est_stripped(centre) + phase_amb; + + % keep phase est in range of -pi .. pi to aid plotting + + if (phase_est(centre) < -pi) + phase_est(centre) += 2*pi; + end + if (phase_est(centre) > pi) + phase_est(centre) -= 2*pi; + end + + rx_symb(centre) *= exp(-j*phase_est(centre)); + end end end end + % if using block based phase est, strip off first few and last few + % symbols where phase ests are invalid + + st = phase_est_delay+1; + %en = Nbits-phase_est_delay; + en = 50; rx_bits = real(rx_symb) < 0; errors = xor(tx_bits, rx_bits); - Nerrs = sum(errors); - printf("EbNodB: %3.2f BER: %4.3f Nbits: %d Nerrs: %d\n", EbNodB(nn), Nerrs/Nbits, Nbits, Nerrs); + Nerrs = sum(errors(st:en)); + printf("EbNodB: %3.2f BER: %4.3f Nbits: %d Nerrs: %d\n", EbNodB(nn), Nerrs/(en-st+1), (en-st+1), Nerrs); if verbose figure(1); clf; - plot(rx_symb,'+'); + plot(rx_symb(st:en),'+'); axis([-2 2 -2 2]); figure(2); clf; - plot(stripped,'+'); - %axis([-2 2 -2 2]); + plot(phase_offset_log(st:en),'b+;phase offset;', 'markersize', 10); + hold on; + plot(phase_amb_log(st:en),'c-;phase amb;'); + plot(phase_est_stripped(st:en),'ro;phase est stripped;'); + plot(phase_est(st:en),'g*;phase est;'); + hold off; + axis([1 en-st+1 -pi pi]) + legend('boxoff'); + figure(3); clf; + plot(errors(st:en),'+'); + axis([1 en-st+1 -0.5 1.5]); end sim_out.ber(nn) = Nerrs/Nbits; @@ -163,14 +259,14 @@ end function run_single - sim_in.Nbits = 100; - sim_in.EbNodB = 40; + sim_in.Nbits = 1000; + sim_in.EbNodB = 4; sim_in.dbpsk = 0; sim_in.m = 2; - sim_in.phase_est_mem = 10; + sim_in.phase_est_mem = 11; sim_in.verbose = 1; - sim_in.phase_offset = pi/4; - sim_in.freq_offset_hz = .01; + sim_in.phase_offset = 0; + sim_in.freq_offset_hz = 1; run_sim(sim_in); end -- 2.25.1