From 14826b498001677e44ea11a0de893ec2561ffc28 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Sat, 28 Mar 2015 22:43:19 +0000 Subject: [PATCH] octave modem with freq offset est working OK, good BER, time to re-port C code git-svn-id: https://svn.code.sf.net/p/freetel/code@2094 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/cohpsk.m | 149 ++++++++---------------------------- codec2-dev/octave/tcohpsk.m | 69 +++++++++++++---- 2 files changed, 86 insertions(+), 132 deletions(-) diff --git a/codec2-dev/octave/cohpsk.m b/codec2-dev/octave/cohpsk.m index a8e990c9..f7978dd5 100644 --- a/codec2-dev/octave/cohpsk.m +++ b/codec2-dev/octave/cohpsk.m @@ -93,6 +93,8 @@ function sim_in = symbol_rate_init(sim_in) sim_in.Nct_sym_buf = 2*Nsymbrowpilot + 2; sim_in.ct_symb_buf = zeros(sim_in.Nct_sym_buf, Nc); + sim_in.ff_phase = 1; + % Init LDPC -------------------------------------------------------------------- if ldpc_code @@ -189,121 +191,42 @@ end % Symbol rate processing for rx side (demodulator) ------------------------------------------------------- -function [rx_symb rx_bits rx_symb_linear amp_linear amp_ phi_ EsNo_ prev_sym_rx sim_in] = qpsk_symbols_to_bits(sim_in, s_ch, prev_sym_rx) - framesize = sim_in.framesize; - Nsymb = sim_in.Nsymb; - Nsymbrow = sim_in.Nsymbrow; - Nsymbrowpilot = sim_in.Nsymbrowpilot; - Nc = sim_in.Nc; - Npilotsframe = sim_in.Npilotsframe; - Ns = sim_in.Ns; - Np = sim_in.Np; - Nchip = sim_in.Nchip; - modulation = sim_in.modulation; - pilot = sim_in.pilot; - rx_symb_buf = sim_in.rx_symb_buf; - rx_pilot_buf = sim_in.rx_pilot_buf; - tx_pilot_buf = sim_in.tx_pilot_buf; - verbose = sim_in.verbose; - - % demodulate stage 1 - - for r=1:Nsymbrowpilot - for c=1:Nc*Nchip - rx_symb(r,c) = s_ch(r, c); - if strcmp(modulation,'dqpsk') - tmp = rx_symb(r,c); - rx_symb(r,c) *= conj(prev_sym_rx(c)/abs(prev_sym_rx(c))); - prev_sym_rx(c) = tmp; - end - end - end - - % strip out pilots - - rx_symb_pilot = rx_symb; - rx_symb = zeros(Nsymbrow, Nc*Nchip); - rx_pilot = zeros(Npilotsframe, Nc*Nchip); +function [rx_symb rx_bits amp_ phi_ EsNo_ cohpsk] = qpsk_symbols_to_bits(cohpsk, ct_symb_buf) + framesize = cohpsk.framesize; + Nsymb = cohpsk.Nsymb; + Nsymbrow = cohpsk.Nsymbrow; + Nsymbrowpilot = cohpsk.Nsymbrowpilot; + Nc = cohpsk.Nc; + Npilotsframe = cohpsk.Npilotsframe; + pilot = cohpsk.pilot; + verbose = cohpsk.verbose; + + % average pilots to get phase and amplitude estimates + % we assume there are two samples at the start of each frame and two at the end + + sampling_points = [1 2 7 8]; + pilot2 = [cohpsk.pilot(1,:); cohpsk.pilot(2,:); cohpsk.pilot(1,:); cohpsk.pilot(2,:);]; + phi_ = zeros(Nsymbrow, Nc); + amp_ = zeros(Nsymbrow, Nc); + + for c=1:Nc + corr = pilot2(:,c)' * ct_symb_buf(sampling_points,c); + mag = sum(abs(ct_symb_buf(sampling_points,c))); - for p=1:Npilotsframe - % printf("%d %d %d %d %d\n", (p-1)*Ns+1, p*Ns, (p-1)*(Ns+1)+2, p*(Ns+1), (p-1)*(Ns+1)+1); - rx_symb((p-1)*Ns+1:p*Ns,:) = rx_symb_pilot((p-1)*(Ns+1)+2:p*(Ns+1),:); - rx_pilot(p,:) = rx_symb_pilot((p-1)*(Ns+1)+1,:); + phi_(:, c) = angle(corr); + amp_(:, c) = mag/length(sampling_points); end - % buffer three frames of symbols (and pilots) for phase recovery - - rx_symb_buf(1:2*Nsymbrow,:) = rx_symb_buf(Nsymbrow+1:3*Nsymbrow,:); - rx_symb_buf(2*Nsymbrow+1:3*Nsymbrow,:) = rx_symb; - rx_pilot_buf(1:2*Npilotsframe,:) = rx_pilot_buf(Npilotsframe+1:3*Npilotsframe,:); - rx_pilot_buf(2*Npilotsframe+1:3*Npilotsframe,:) = rx_pilot; - sim_in.rx_symb_buf = rx_symb_buf; - sim_in.rx_pilot_buf = rx_pilot_buf; - - % pilot assisted phase estimation and correction of middle frame in rx symb buffer + % now correct phase of data symbols and make decn on bits - rx_symb = rx_symb_buf(Nsymbrow+1:2*Nsymbrow,:); - - phi_ = zeros(Nsymbrow, Nc*Nchip); - amp_ = ones(Nsymbrow, Nc*Nchip); - - for c=1:Nc*Nchip - - if verbose > 2 - printf("phi_ : "); - end - - for r=1:Nsymbrow - st = Npilotsframe+1+floor((r-1)/Ns) - floor(Np/2) + 1; - en = st + Np - 1; - ch_est = tx_pilot_buf(st:en,c)'*rx_pilot_buf(st:en,c)/Np; - phi_(r,c) = angle(ch_est); - amp_(r,c) = abs(ch_est); - %amp_(r,c) = abs(rx_symb(r,c)); - if verbose > 2 - printf("% 4.3f ", phi_(r,c)) - end - rx_symb(r,c) *= exp(-j*phi_(r,c)); - end - - if verbose > 2 - printf("\nrx_symb: "); - for r=1:Nsymbrow - printf("% 4.3f ", angle(rx_symb(r,c))) - end - printf("\nindexes: "); - for r=1:Nsymbrow - st = Npilotsframe+1+floor((r-1)/Ns) - floor(Np/2) + 1; - en = st + Np - 1; - printf("%2d,%2d ", st,en) - end - printf("\npilots : "); - for p=1:3*Npilotsframe - printf("% 4.3f ", angle(rx_pilot_buf(p,c))); - end - printf("\n\n"); - end - end - - % de-spread - - for r=1:Nsymbrow - for c=Nc+1:Nc:Nchip*Nc - rx_symb(r,1:Nc) = rx_symb(r,1:Nc) + rx_symb(r,c:c+Nc-1); - amp_(r,1:Nc) = amp_(r,1:Nc) + amp_(r,c:c+Nc-1); - end - end - - % demodulate stage 2 - - rx_symb_linear = zeros(1,Nsymb); - amp_linear = zeros(1,Nsymb); + rx_symb = zeros(Nsymbrow, Nc); + rx_symb_linear = zeros(1, Nsymbrow*Nc); rx_bits = zeros(1, framesize); for c=1:Nc for r=1:Nsymbrow i = (c-1)*Nsymbrow + r; - rx_symb_linear(i) = rx_symb(r,c); - amp_linear(i) = amp_(r,c); + rx_symb(i) = ct_symb_buf(2+r,c)*exp(-j*phi_(c)); + rx_symb_linear(i) = rx_symb(i); rx_bits((2*(i-1)+1):(2*i)) = qpsk_demod(rx_symb(r,c)); end end @@ -337,18 +260,11 @@ function [rx_symb rx_bits rx_symb_linear amp_linear amp_ phi_ EsNo_ prev_sym_rx % Estimate signal power - Es_ = mean(amp_linear .^ 2); + Es_ = mean(amp_ .^ 2); EsNo_ = Es_/No_; %printf("Es_: %f No_: %f Es/No: %f Es/No dB: %f\n", Es_, No_, Es_/No_, 10*log10(EsNo_)); - - % LDPC decoder requires some amplitude normalisation - % (AGC), was found to break ow. So we adjust the symbol - % amplitudes so that they are an averge of 1 - - rx_symb_linear /= mean(amp_linear); - amp_linear /= mean(amp_linear); - + endfunction @@ -625,6 +541,7 @@ function [next_sync cohpsk] = frame_sync_fine_timing_est(cohpsk, ch_symb, sync, max_mag = mag; cohpsk.ct = t; cohpsk.f_fine_est = f_fine; + cohpsk.ff_rect = exp(-j*f_fine*2*pi/Rs); end end end diff --git a/codec2-dev/octave/tcohpsk.m b/codec2-dev/octave/tcohpsk.m index 5f1e26c6..593dbfb1 100644 --- a/codec2-dev/octave/tcohpsk.m +++ b/codec2-dev/octave/tcohpsk.m @@ -25,7 +25,7 @@ randn('state',1); n = 2000; frames = 35*4; framesize = 32; -foff = -80; +foff = 0; EsNodB = 8; EsNo = 10^(EsNodB/10); @@ -233,36 +233,71 @@ for i=1:frames % coarse timing (frame sync) and initial fine freq est --------------------------------------------- [next_sync acohpsk] = frame_sync_fine_timing_est(acohpsk, ch_symb, sync, next_sync); + %acohpsk.ff_rect = exp(j*2*pi*(-1.0)/Rs); if sync == 1 next_sync = 2; end - printf("i: %d sync: %d next_sync: %d\n", i, sync, next_sync); + %printf("i: %d sync: %d next_sync: %d\n", i, sync, next_sync); sync = next_sync; - if (i==10) + if (i==1000) xx end - if 0 - [rx_symb rx_bits rx_symb_linear amp_linear amp_ phi_ EsNo_ prev_sym_rx acohpsk] = qpsk_symbols_to_bits(acohpsk, ct_symb_buf(ct+1:ct+acohpsk.Nsymbrowpilot,:), []); - rx_symb_log = [rx_symb_log; rx_symb]; - rx_amp_log = [rx_amp_log; amp_]; - rx_phi_log = [rx_phi_log; phi_]; - rx_bits_log = [rx_bits_log rx_bits]; + % We can decode first that we get sync on. Need to fine freq correct all of it's symbols, + % including pilots. From then on, just correct new symbols into frame. + % make copy, so if we lose sync we havent fine freq corrected ct_symb_buf + % if next_sync == 4 correct all 8 + % if sync == 2 correct latest 6 - % BER stats + if (next_sync == 4) || (sync == 4) - if i > 3 - error_positions = xor(prev_tx_bits2, rx_bits); - Nerrs += sum(error_positions); - nerr_log = [nerr_log sum(error_positions)]; - Tbits += length(error_positions); + if next_sync == 4 + + % first frame, we've just gotten sync so fine freq correct all Nsymbrowpilot+2 samples + + ct_symb_ff_buf = acohpsk.ct_symb_buf(acohpsk.ct+1:acohpsk.ct+acohpsk.Nsymbrowpilot+2,:); + for r=1:acohpsk.Nsymbrowpilot+2 + acohpsk.ff_phase *= acohpsk.ff_rect'; + ct_symb_ff_buf(r,:) *= acohpsk.ff_phase; + end + + else + + % second and subsequent frames, just fine freq correct the latest Nsymbrowpilot + + ct_symb_ff_buf(1:2,:) = ct_symb_ff_buf(Nsymbrowpilot+1:Nsymbrowpilot+2,:); + ct_symb_ff_buf(3:Nsymbrowpilot+2,:) = acohpsk.ct_symb_buf(acohpsk.ct+3:acohpsk.ct+acohpsk.Nsymbrowpilot+2,:); + for r=3:acohpsk.Nsymbrowpilot+2 + afdmdv.ff_phase *= afdmdv.ff_rect'; + ct_symb_ff_buf(r,:) *= afdmdv.ff_phase; + end + + end + + mag = abs(acohpsk.ff_phase); + acohpsk.ff_phase /= mag; + + [rx_symb rx_bits amp_ phi_ EsNo_] = qpsk_symbols_to_bits(acohpsk, ct_symb_ff_buf); + rx_symb_log = [rx_symb_log; rx_symb]; + rx_amp_log = [rx_amp_log; amp_]; + rx_phi_log = [rx_phi_log; phi_]; + rx_bits_log = [rx_bits_log rx_bits]; + + % BER stats + + if i > 2 + error_positions = xor(prev_tx_bits2, rx_bits); + Nerrs += sum(error_positions); + nerr_log = [nerr_log sum(error_positions)]; + Tbits += length(error_positions); + end end + prev_tx_bits2 = prev_tx_bits; prev_tx_bits = tx_bits; - end end @@ -281,6 +316,7 @@ stem_sig_and_error(6, 211, real(rx_symb_log_c(1:n)), real(rx_symb_log(1:n) - rx_ stem_sig_and_error(6, 212, imag(rx_symb_log_c(1:n)), imag(rx_symb_log(1:n) - rx_symb_log_c(1:n)), 'rx symb im', [1 n -1.5 1.5]) stem_sig_and_error(7, 111, rx_bits_log_c(1:n), rx_bits_log(1:n) - rx_bits_log_c(1:n), 'rx bits', [1 n -1.5 1.5]) +if 0 check(tx_bits_log, tx_bits_log_c, 'tx_bits'); check(tx_symb_log, tx_symb_log_c, 'tx_symb'); check(tx_fdm_log, tx_fdm_log_c, 'tx_fdm'); @@ -292,6 +328,7 @@ check(rx_amp_log, rx_amp_log_c, 'rx_amp_log',0.01); check(rx_phi_log, rx_phi_log_c, 'rx_phi_log'); check(rx_symb_log, rx_symb_log_c, 'rx_symb',0.01); check(rx_bits_log, rx_bits_log_c, 'rx_bits'); +end % Determine bit error rate -- 2.25.1