From ac3a26fce0fde7507e2340777b1e5406d7e5c8ca Mon Sep 17 00:00:00 2001 From: drowe67 Date: Mon, 4 May 2015 04:40:37 +0000 Subject: [PATCH] working on freq estimation, about to try two path model git-svn-id: https://svn.code.sf.net/p/freetel/code@2122 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/cohpsk.m | 143 +++++++++++--- codec2-dev/octave/tcohpsk.m | 191 +++++++++++------- codec2-dev/octave/test_cohpsk.m | 75 +------ codec2-dev/octave/test_foff.m | 338 ++++++++++++++++++++++++++++++++ 4 files changed, 580 insertions(+), 167 deletions(-) create mode 100644 codec2-dev/octave/test_foff.m diff --git a/codec2-dev/octave/cohpsk.m b/codec2-dev/octave/cohpsk.m index c93cca3d..42329032 100644 --- a/codec2-dev/octave/cohpsk.m +++ b/codec2-dev/octave/cohpsk.m @@ -534,33 +534,58 @@ endfunction % returns an estimate of frequency offset, advances to next sync state function [next_sync cohpsk] = coarse_freq_offset_est(cohpsk, fdmdv, ch_fdm_frame, sync, next_sync) - Fcentre = fdmdv.Fcentre; - Nc = fdmdv.Nc; - Fsep = fdmdv.Fsep; - M = fdmdv.M; - Fs = fdmdv.Fs; - Ndft = cohpsk.Ndft; + Fcentre = fdmdv.Fcentre; + Nc = fdmdv.Nc; + Fsep = fdmdv.Fsep; + M = fdmdv.M; + Fs = fdmdv.Fs; + Ndft = cohpsk.Ndft; + coarse_mem = cohpsk.coarse_mem; + Ncm = cohpsk.Ncm; + + % ll + % |--------|-----| + + ll = length(ch_fdm_frame); + sz_mem = Ncm-ll; + for i=1:sz_mem + coarse_mem(i) = coarse_mem(i+ll); + end + coarse_mem(Ncm-ll+1:Ncm) = ch_fdm_frame; if sync == 0 - f_start = Fcentre - ((Nc/2)+2)*Fsep; - f_stop = Fcentre + ((Nc/2)+2)*Fsep; - ll = length(ch_fdm_frame); - h = 0.5 - 0.5*cos(2*pi*(0:ll-1)/(ll-1)); - T = abs(fft(ch_fdm_frame .* h, Ndft)).^2; + h = 0.5 - 0.5*cos(2*pi*(0:Ncm-1)/(Ncm-1)); + T = abs(fft(coarse_mem .* h, Ndft)).^2; sc = Ndft/Fs; - bin_start = floor(f_start*sc+0.5)+1; - bin_stop = floor(f_stop*sc+0.5)+1; + for i=1:5 + f_start = Fcentre - ((Nc/2)+2)*Fsep; + f_stop = Fcentre + ((Nc/2)+2)*Fsep; + bin_start = floor(f_start*sc+0.5)+1; + bin_stop = floor(f_stop*sc+0.5)+1; + x = bin_start-1:bin_stop-1; + bin_est = x*T(bin_start:bin_stop)'/sum(T(bin_start:bin_stop)); + f_est = floor(bin_est/sc+0.5); + Fcentre = f_est; + end %printf("f_start: %f f_stop: %f sc: %f bin_start: %d bin_stop: %d\n", f_start, f_stop, sc, bin_start, bin_stop); - x = bin_start-1:bin_stop-1; - bin_est = x*T(bin_start:bin_stop)'/sum(T(bin_start:bin_stop)); - cohpsk.f_est = floor(bin_est/sc+0.5); + cohpsk.f_est = f_est; + printf(" coarse freq est: %f\n", cohpsk.f_est); next_sync = 1; - + figure(5) + clf + subplot(211) + plot(T) + hold on + plot([bin_est bin_est],[0 max(T)],'g') + hold off + axis([bin_start bin_stop 0 max(T)]) + end + cohpsk.coarse_mem = coarse_mem; endfunction @@ -596,7 +621,7 @@ function [next_sync cohpsk] = frame_sync_fine_freq_est(cohpsk, ch_symb, sync, ne % sample correlation over 2D grid of time and fine freq points max_corr = 0; - for f_fine=-20:1:20 + for f_fine=-20:0.25:20 f_fine_rect = exp(-j*f_fine*2*pi*sampling_points/Rs)'; for t=0:cohpsk.Nsymbrowpilot-1 corr = 0; mag = 0; @@ -619,14 +644,52 @@ function [next_sync cohpsk] = frame_sync_fine_freq_est(cohpsk, ch_symb, sync, ne end printf(" fine freq f: %f max_corr: %f max_mag: %f ct: %d\n", cohpsk.f_fine_est, abs(max_corr), max_mag, cohpsk.ct); - if max_corr/max_mag > 0.9 - printf(" in sync!\n"); + if abs(max_corr/max_mag) > 0.7 + printf(" [%d] in sync!\n", cohpsk.frame); + cohpsk.sync_timer = 0; + %cohpsk.f_est -= cohpsk.f_fine_est; + %cohpsk.f_fine_est = 0; + %cohpsk.ff_rect = 1; + printf(" .... adjusting to %f\n", cohpsk.f_est); next_sync = 4; else next_sync = 0; printf(" back to coarse freq offset est...\n"); end + cohpsk.ratio = abs(max_corr/max_mag); + end + + + if sync == 4 + + % we are in sync so just sample correlation over 1D grid of fine freq points + + max_corr = 0; + st = cohpsk.f_fine_est - 1; + en = cohpsk.f_fine_est + 1; + for f_fine = st:0.25*en + f_fine_rect = exp(-j*f_fine*2*pi*sampling_points/Rs)'; + corr = 0; mag = 0; + for c=1:Nc*Nd + f_corr_vec = f_fine_rect .* ct_symb_buf(cohpsk.ct+sampling_points,c); + for p=1:length(sampling_points) + corr += pilot2(p,c-Nc*floor((c-1)/Nc)) * f_corr_vec(p); + mag += abs(f_corr_vec(p)); + end + end + if corr > max_corr + max_corr = corr; + max_mag = mag; + f_fine_est = f_fine; + end + end + + %cohpsk.f_est -= 0.5*f_fine_est; + %printf(" coarse: %f fine: %f\n", cohpsk.f_est, f_fine_est); + cohpsk.ratio = abs(max_corr/max_mag); end + + endfunction @@ -641,9 +704,8 @@ function acohpsk = fine_freq_correct(acohpsk, sync, next_sync); % havent fine freq corrected ct_symb_buf if next_sync == 4 correct % all 8 if sync == 2 correct latest 6 - if (next_sync == 4) || (sync == 4) - if (next_sync == 4) && (sync == 2) + if (next_sync == 4) && (sync == 2) % first frame, we've just gotten sync so fine freq correct all Nsymbrowpilot+2 samples @@ -652,9 +714,9 @@ function acohpsk = fine_freq_correct(acohpsk, sync, next_sync); acohpsk.ff_phase *= acohpsk.ff_rect'; ct_symb_ff_buf(r,:) *= acohpsk.ff_phase; end + end - else - + if sync == 4 % second and subsequent frames, just fine freq correct the latest Nsymbrowpilot ct_symb_ff_buf(1:2,:) = ct_symb_ff_buf(acohpsk.Nsymbrowpilot+1:acohpsk.Nsymbrowpilot+2,:); @@ -663,23 +725,44 @@ function acohpsk = fine_freq_correct(acohpsk, sync, next_sync); acohpsk.ff_phase *= acohpsk.ff_rect'; ct_symb_ff_buf(r,:) *= acohpsk.ff_phase; end + end - end + mag = abs(acohpsk.ff_phase); + acohpsk.ff_phase /= mag; - mag = abs(acohpsk.ff_phase); - acohpsk.ff_phase /= mag; + acohpsk.ct_symb_ff_buf = ct_symb_ff_buf; - acohpsk.ct_symb_ff_buf = ct_symb_ff_buf; - end endfunction % misc sync state machine code, just wanted it in a function -function sync = sync_state_machine(sync, next_sync) +function [sync cohpsk] = sync_state_machine(cohpsk, sync, next_sync) + if sync == 1 next_sync = 2; end + if sync == 5 + next_sync = 4; + end + + if sync == 4 + + % check that sync is still good, fall out of sync on consecutive bad frames */ + + if cohpsk.ratio < 0.5 + cohpsk.sync_timer++; + else + cohpsk.sync_timer = 0; + end + %printf(" ratio: %f sync timer: %d\n", cohpsk.ratio, cohpsk.sync_timer); + + if cohpsk.sync_timer == 5 + printf(" lost sync ....\n"); + next_sync = 0; + end + end + sync = next_sync; endfunction diff --git a/codec2-dev/octave/tcohpsk.m b/codec2-dev/octave/tcohpsk.m index 72aff531..d401766c 100644 --- a/codec2-dev/octave/tcohpsk.m +++ b/codec2-dev/octave/tcohpsk.m @@ -1,14 +1,17 @@ % tcohpsk.m % David Rowe Oct 2014 % -% Octave script that tests the C port of the coherent PSK modem. This -% script loads the output of unittest/tcohpsk.c and compares it to the -% output of the reference versions of the same functions written in -% Octave. +% Octave script that: % - +% i) tests the C port of the coherent PSK modem. This script loads +% the output of unittest/tcohpsk.c and compares it to the output of +% the reference versions of the same functions written in Octave. +% +% or (ii) can optionally be used to run an Octave version of the cohpsk +% modem to tune and develop it. graphics_toolkit ("gnuplot"); +more off; cohpsk; fdmdv; @@ -19,10 +22,13 @@ randn('state',1); % test parameters ---------------------------------------------------------- -n = 840; -frames = 35; -foff = 10.5; -EsNodB = 8; +frames = 100; +foff = 100; +dfoff = 0; +EsNodB = 10; +fading_en = 0; +hf_delay_ms = 2; +compare_with_c = 0; EsNo = 10^(EsNodB/10); @@ -108,8 +114,10 @@ rx_fdm_filter_log = []; rx_baseband_log = []; rx_fdm_frame_log = []; ct_symb_ff_log = []; +rx_timing_log = []; +ratio_log = []; -% BER measurement ------------------------------------------------------------- +% Channel modeling and BER measurement ---------------------------------------- rand('state',1); tx_bits_coh = round(rand(1,framesize*10)); @@ -121,13 +129,15 @@ prev_tx_bits = []; phase_ch = 1; sync = 0; -% Output vectors from C port --------------------------------------------------- - -load ../build_linux/unittest/tcohpsk_out.txt +[spread spread_2ms hf_gain] = init_hf_model(Fs, Fs, frames*acohpsk.Nsymbrowpilot*afdmdv.M); +hf_n = 1; +nhfdelay = floor(hf_delay_ms*Fs/1000); +ch_fdm_delay = zeros(1, acohpsk.Nsymbrowpilot*M + nhfdelay); % main loop -------------------------------------------------------------------- for i=1:frames + acohpsk.frame = i; tx_bits = tx_bits_coh(ptx_bits_coh:ptx_bits_coh+framesize-1); ptx_bits_coh += framesize; if ptx_bits_coh > length(tx_bits_coh) @@ -153,8 +163,29 @@ for i=1:frames % Channel -------------------------------------------------------------------- % - [ch_fdm_frame phase_ch] = freq_shift(tx_fdm_frame, foff, Fs, phase_ch); - + %[ch_fdm_frame phase_ch] = freq_shift(tx_fdm_frame, foff, Fs, phase_ch); + %foff += dfoff*acohpsk.Nsymbrowpilot*M; + + ch_fdm_frame = zeros(1,acohpsk.Nsymbrowpilot*M); + for i=1:acohpsk.Nsymbrowpilot*M + foff_rect = exp(j*2*pi*foff/Fs); + foff += dfoff; + phase_ch *= foff_rect; + ch_fdm_frame(i) = tx_fdm_frame(i) * phase_ch; + end + phase_ch /= abs(phase_ch); + + if fading_en + ch_fdm_delay(1:nhfdelay) = ch_fdm_delay(acohpsk.Nsymbrowpilot*M+1:nhfdelay+acohpsk.Nsymbrowpilot*M); + ch_fdm_delay(nhfdelay+1:nhfdelay+acohpsk.Nsymbrowpilot*M) = ch_fdm_frame; + + for i=1:acohpsk.Nsymbrowpilot*M + ahf_model = hf_gain*(spread(hf_n)*ch_fdm_frame(i) + spread_2ms(hf_n)*ch_fdm_delay(i)); + ch_fdm_frame(i) = ahf_model; + hf_n++; + end + end + % each carrier has power = 2, total power 2Nc, total symbol rate NcRs, noise BW B=Fs % Es/No = (C/Rs)/(N/B), N = var = 2NcFs/NcRs(Es/No) = 2Fs/Rs(Es/No) @@ -176,7 +207,6 @@ for i=1:frames next_sync = sync; [next_sync acohpsk] = coarse_freq_offset_est(acohpsk, afdmdv, ch_fdm_frame, sync, next_sync); - nin = M; % shift entire FDM signal to 0 Hz @@ -199,6 +229,8 @@ for i=1:frames rx_filt_log = [rx_filt_log rx_filt]; [rx_onesym rx_timing env afdmdv] = rx_est_timing(afdmdv, rx_filt, nin); + rx_timing_log = [rx_timing_log rx_timing]; + ch_symb(r,:) = rx_onesym; end @@ -217,6 +249,7 @@ for i=1:frames rx_phi_log = [rx_phi_log; phi_]; rx_bits_log = [rx_bits_log rx_bits]; tx_bits_prev_log = [tx_bits_prev_log prev_tx_bits2]; + ratio_log = [ratio_log acohpsk.ratio]; % BER stats @@ -226,68 +259,90 @@ for i=1:frames Tbits += length(error_positions); end - %printf("f: %d sync: %d next_sync: %d\n", i, sync, next_sync); - sync = sync_state_machine(sync, next_sync); + if sync == 0 + Nerrs = 0; + Tbits = 0; + nerr_log = []; + end + + % printf("f: %d sync: %d next_sync: %d\n", i, sync, next_sync); + [sync acohpsk] = sync_state_machine(acohpsk, sync, next_sync); prev_tx_bits2 = prev_tx_bits; prev_tx_bits = tx_bits; end -stem_sig_and_error(1, 111, tx_bits_log_c(1:n), tx_bits_log(1:n) - tx_bits_log_c(1:n), 'tx bits', [1 n -1.5 1.5]) -stem_sig_and_error(2, 211, real(tx_symb_log_c(1:n)), real(tx_symb_log(1:n) - tx_symb_log_c(1:n)), 'tx symb re', [1 n -1.5 1.5]) -stem_sig_and_error(2, 212, imag(tx_symb_log_c(1:n)), imag(tx_symb_log(1:n) - tx_symb_log_c(1:n)), 'tx symb im', [1 n -1.5 1.5]) - -stem_sig_and_error(3, 211, real(tx_fdm_frame_log_c(1:n)), real(tx_fdm_frame_log(1:n) - tx_fdm_frame_log_c(1:n)), 'tx fdm frame re', [1 n -10 10]) -stem_sig_and_error(3, 212, imag(tx_fdm_frame_log_c(1:n)), imag(tx_fdm_frame_log(1:n) - tx_fdm_frame_log_c(1:n)), 'tx fdm frame im', [1 n -10 10]) -stem_sig_and_error(4, 211, real(ch_fdm_frame_log_c(1:n)), real(ch_fdm_frame_log(1:n) - ch_fdm_frame_log_c(1:n)), 'ch fdm frame re', [1 n -10 10]) -stem_sig_and_error(4, 212, imag(ch_fdm_frame_log_c(1:n)), imag(ch_fdm_frame_log(1:n) - ch_fdm_frame_log_c(1:n)), 'ch fdm frame im', [1 n -10 10]) -stem_sig_and_error(5, 211, real(rx_fdm_frame_bb_log_c(1:n)), real(rx_fdm_frame_bb_log(1:n) - rx_fdm_frame_bb_log_c(1:n)), 'rx fdm frame bb re', [1 n -10 10]) -stem_sig_and_error(5, 212, imag(rx_fdm_frame_bb_log_c(1:n)), imag(rx_fdm_frame_bb_log(1:n) - rx_fdm_frame_bb_log_c(1:n)), 'rx fdm frame bb im', [1 n -10 10]) - -[n m] = size(ch_symb_log); -stem_sig_and_error(6, 211, real(ch_symb_log_c), real(ch_symb_log - ch_symb_log_c), 'ch symb re', [1 n -1.5 1.5]) -stem_sig_and_error(6, 212, imag(ch_symb_log_c), imag(ch_symb_log - ch_symb_log_c), 'ch symb im', [1 n -1.5 1.5]) -stem_sig_and_error(7, 211, real(ct_symb_ff_log_c), real(ct_symb_ff_log - ct_symb_ff_log_c), 'ct symb ff re', [1 n -1.5 1.5]) -stem_sig_and_error(7, 212, imag(ct_symb_ff_log_c), imag(ct_symb_ff_log - ct_symb_ff_log_c), 'ct symb ff im', [1 n -1.5 1.5]) -stem_sig_and_error(8, 211, rx_amp_log_c, rx_amp_log - rx_amp_log_c, 'Amp Est', [1 n -1.5 1.5]) -stem_sig_and_error(8, 212, rx_phi_log_c, rx_phi_log - rx_phi_log_c, 'Phase Est', [1 n -4 4]) -stem_sig_and_error(9, 211, real(rx_symb_log_c), real(rx_symb_log - rx_symb_log_c), 'rx symb re', [1 n -1.5 1.5]) -stem_sig_and_error(9, 212, imag(rx_symb_log_c), imag(rx_symb_log - rx_symb_log_c), 'rx symb im', [1 n -1.5 1.5]) -stem_sig_and_error(10, 111, rx_bits_log_c, rx_bits_log - rx_bits_log_c, 'rx bits', [1 n -1.5 1.5]) - -check(tx_bits_log, tx_bits_log_c, 'tx_bits'); -check(tx_symb_log, tx_symb_log_c, 'tx_symb'); -check(tx_fdm_frame_log, tx_fdm_frame_log_c, 'tx_fdm_frame'); -check(ch_fdm_frame_log, ch_fdm_frame_log_c, 'ch_fdm_frame'); -check(rx_fdm_frame_bb_log, rx_fdm_frame_bb_log_c, 'rx_fdm_frame_bb', 0.01); -check(ch_symb_log, ch_symb_log_c, 'ch_symb',0.01); -check(ct_symb_ff_log, ct_symb_ff_log_c, 'ct_symb_ff',0.01); -check(rx_amp_log, rx_amp_log_c, 'rx_amp_log',0.01); -check(rx_phi_log, rx_phi_log_c, 'rx_phi_log',0.01); -check(rx_symb_log, rx_symb_log_c, 'rx_symb',0.01); -check(rx_bits_log, rx_bits_log_c, 'rx_bits'); - -% Determine bit error rate - -sz = length(tx_bits_log_c); -Nerrs_c = sum(xor(tx_bits_prev_log, rx_bits_log_c)); -Tbits_c = length(tx_bits_prev_log); -ber_c = Nerrs_c/Tbits_c; ber = Nerrs/Tbits; +printf("EsNodB: %4.1f ber..: %4.3f Nerrs..: %d Tbits..: %d\n", EsNodB, ber, Nerrs, Tbits); + +if compare_with_c + + % Output vectors from C port --------------------------------------------------- + + load ../build_linux/unittest/tcohpsk_out.txt + + stem_sig_and_error(1, 111, tx_bits_log_c(1:n), tx_bits_log(1:n) - tx_bits_log_c(1:n), 'tx bits', [1 n -1.5 1.5]) + stem_sig_and_error(2, 211, real(tx_symb_log_c(1:n)), real(tx_symb_log(1:n) - tx_symb_log_c(1:n)), 'tx symb re', [1 n -1.5 1.5]) + stem_sig_and_error(2, 212, imag(tx_symb_log_c(1:n)), imag(tx_symb_log(1:n) - tx_symb_log_c(1:n)), 'tx symb im', [1 n -1.5 1.5]) + + stem_sig_and_error(3, 211, real(tx_fdm_frame_log_c(1:n)), real(tx_fdm_frame_log(1:n) - tx_fdm_frame_log_c(1:n)), 'tx fdm frame re', [1 n -10 10]) + stem_sig_and_error(3, 212, imag(tx_fdm_frame_log_c(1:n)), imag(tx_fdm_frame_log(1:n) - tx_fdm_frame_log_c(1:n)), 'tx fdm frame im', [1 n -10 10]) + stem_sig_and_error(4, 211, real(ch_fdm_frame_log_c(1:n)), real(ch_fdm_frame_log(1:n) - ch_fdm_frame_log_c(1:n)), 'ch fdm frame re', [1 n -10 10]) + stem_sig_and_error(4, 212, imag(ch_fdm_frame_log_c(1:n)), imag(ch_fdm_frame_log(1:n) - ch_fdm_frame_log_c(1:n)), 'ch fdm frame im', [1 n -10 10]) + stem_sig_and_error(5, 211, real(rx_fdm_frame_bb_log_c(1:n)), real(rx_fdm_frame_bb_log(1:n) - rx_fdm_frame_bb_log_c(1:n)), 'rx fdm frame bb re', [1 n -10 10]) + stem_sig_and_error(5, 212, imag(rx_fdm_frame_bb_log_c(1:n)), imag(rx_fdm_frame_bb_log(1:n) - rx_fdm_frame_bb_log_c(1:n)), 'rx fdm frame bb im', [1 n -10 10]) + + [n m] = size(ch_symb_log); + stem_sig_and_error(6, 211, real(ch_symb_log_c), real(ch_symb_log - ch_symb_log_c), 'ch symb re', [1 n -1.5 1.5]) + stem_sig_and_error(6, 212, imag(ch_symb_log_c), imag(ch_symb_log - ch_symb_log_c), 'ch symb im', [1 n -1.5 1.5]) + stem_sig_and_error(7, 211, real(ct_symb_ff_log_c), real(ct_symb_ff_log - ct_symb_ff_log_c), 'ct symb ff re', [1 n -1.5 1.5]) + stem_sig_and_error(7, 212, imag(ct_symb_ff_log_c), imag(ct_symb_ff_log - ct_symb_ff_log_c), 'ct symb ff im', [1 n -1.5 1.5]) + stem_sig_and_error(8, 211, rx_amp_log_c, rx_amp_log - rx_amp_log_c, 'Amp Est', [1 n -1.5 1.5]) + stem_sig_and_error(8, 212, rx_phi_log_c, rx_phi_log - rx_phi_log_c, 'Phase Est', [1 n -4 4]) + stem_sig_and_error(9, 211, real(rx_symb_log_c), real(rx_symb_log - rx_symb_log_c), 'rx symb re', [1 n -1.5 1.5]) + stem_sig_and_error(9, 212, imag(rx_symb_log_c), imag(rx_symb_log - rx_symb_log_c), 'rx symb im', [1 n -1.5 1.5]) + stem_sig_and_error(10, 111, rx_bits_log_c, rx_bits_log - rx_bits_log_c, 'rx bits', [1 n -1.5 1.5]) + + check(tx_bits_log, tx_bits_log_c, 'tx_bits'); + check(tx_symb_log, tx_symb_log_c, 'tx_symb'); + check(tx_fdm_frame_log, tx_fdm_frame_log_c, 'tx_fdm_frame'); + check(ch_fdm_frame_log, ch_fdm_frame_log_c, 'ch_fdm_frame'); + check(rx_fdm_frame_bb_log, rx_fdm_frame_bb_log_c, 'rx_fdm_frame_bb', 0.01); + + check(ch_symb_log, ch_symb_log_c, 'ch_symb',0.01); + check(ct_symb_ff_log, ct_symb_ff_log_c, 'ct_symb_ff',0.01); + check(rx_amp_log, rx_amp_log_c, 'rx_amp_log',0.01); + check(rx_phi_log, rx_phi_log_c, 'rx_phi_log',0.01); + check(rx_symb_log, rx_symb_log_c, 'rx_symb',0.01); + check(rx_bits_log, rx_bits_log_c, 'rx_bits'); + + % Determine bit error rate + + sz = length(tx_bits_log_c); + Nerrs_c = sum(xor(tx_bits_prev_log, rx_bits_log_c)); + Tbits_c = length(tx_bits_prev_log); + ber_c = Nerrs_c/Tbits_c; + printf("EsNodB: %4.1f ber_c: %4.3f Nerrs_c: %d Tbits_c: %d\n", EsNodB, ber_c, Nerrs_c, Tbits_c); + +else + + % some other useful plots + + figure(1) + plot(rx_symb_log*exp(j*pi/4),'+') + figure(2) + plot(rx_timing_log) + figure(3) + stem(nerr_log) + figure(4) + stem(ratio_log) -printf("EsNodB: %4.1f ber_c: %3.2f Nerrs_c: %d Tbits_c: %d\n", EsNodB, ber_c, Nerrs_c, Tbits_c); - -% some other useful plots - -if 0 -figure -plot(rx_symb_log,'+') -figure -plot(xor(tx_bits_prev_log, rx_bits_log_c),'+') end -% C header file of noise samples so C version gives extacly the same results + +% function to write C header file of noise samples so C version gives +% extactly the same results function write_noise_file(uvnoise_log) diff --git a/codec2-dev/octave/test_cohpsk.m b/codec2-dev/octave/test_cohpsk.m index c680f267..3b0abed7 100644 --- a/codec2-dev/octave/test_cohpsk.m +++ b/codec2-dev/octave/test_cohpsk.m @@ -147,7 +147,7 @@ function test_single sim_in.ldpc_code = 0; sim_in.Ntrials = 400; - sim_in.Esvec = 13; + sim_in.Esvec = 12; sim_in.hf_sim = 1; sim_in.hf_mag_only = 0; sim_in.modulation = 'qpsk'; @@ -531,79 +531,14 @@ function rate_Fs_rx(rx_filename) endfunction -% ideas: cld estimate timing with freq offset and decimate to save cpu load -% fft to do cross correlation -function [f_max s_max] = test_freq_off_est(rx_filename, offset, n) - fpilot = fopen("tx_zero.raw","rb"); tx_pilot = fread(fpilot, "short"); fclose(fpilot); - frx=fopen(rx_filename,"rb"); rx_fdm = fread(frx, "short"); fclose(frx); - Fs = 8000; - nc = 1800; % portion we wish to correlate over (first 2 rows on pilots) - - % downconvert to complex baseband to remove images - - f = 1000; - foff_rect = exp(j*2*pi*f*(1:2*n)/Fs); - tx_pilot_bb = tx_pilot(1:n) .* foff_rect(1:n)'; - rx_fdm_bb = rx_fdm(offset:offset+2*n-1) .* foff_rect'; - - % remove -2000 Hz image - - b = fir1(50, 1000/Fs); - tx_pilot_bb_lpf = filter(b,1,tx_pilot_bb); - rx_fdm_bb_lpf = filter(b,1,rx_fdm_bb); - - % decimate by M - - M = 4; - tx_pilot_bb_lpf = tx_pilot_bb_lpf(1:M:length(tx_pilot_bb_lpf)); - rx_fdm_bb_lpf = rx_fdm_bb_lpf(1:M:length(rx_fdm_bb_lpf)); - n /= M; - nc /= M; - - % correlate over a range of frequency offsets and delays - - c_max = 0; - f_n = 1; - f_range = -75:2.5:75; - c_log=zeros(n, length(f_range)); - - for f=f_range - foff_rect = exp(j*2*pi*(f*M)*(1:nc)/Fs); - for s=1:n - - c = abs(tx_pilot_bb_lpf(1:nc)' * (rx_fdm_bb_lpf(s:s+nc-1) .* foff_rect')); - c_log(s,f_n) = c; - if c > c_max - c_max = c; - f_max = f; - s_max = s; - end - end - f_n++; - %printf("f: %f c_max: %f f_max: %f s_max: %d\n", f, c_max, f_max, s_max); - end - - figure(1); - y = f_range; - x = max(s_max-25,1):min(s_max+25, n); - mesh(y,x, c_log(x,:)); - grid - - s_max *= M; - s_max -= floor(s_max/6400)*6400; - printf("f_max: %f s_max: %d\n", f_max, s_max); - - % decimated position at sample rate. need to relate this to symbol - % rate position. - -endfunction function snr = EsNo_to_SNR(EsNo) snr = interp1([20 11.8 9.0 6.7 4.9 3.3 -10], [ 3 3 0 -3 -6 -9 -9], EsNo); endfunction + function gen_test_bits() sim_in = standard_init(); framesize = 32*10; @@ -611,16 +546,18 @@ function gen_test_bits() test_bits_coh_file(tx_bits); endfunction + + + % Start simulations --------------------------------------- more off; %close all; %test_curves(); -test_single(); +%test_single(); %rate_Fs_tx("tx_zero.raw"); %rate_Fs_tx("tx.raw"); %rate_Fs_rx("tx_-4dB.wav") %rate_Fs_rx("tx.raw") -%test_freq_off_est("tx.raw",40,6400) %gen_test_bits(); diff --git a/codec2-dev/octave/test_foff.m b/codec2-dev/octave/test_foff.m new file mode 100644 index 00000000..4745fdee --- /dev/null +++ b/codec2-dev/octave/test_foff.m @@ -0,0 +1,338 @@ +% test_foff.m +% David Rowe April 2015 +% +% Octave script for testing the cohpsk freq offset estimator + +graphics_toolkit ("gnuplot"); +more off; + +cohpsk; +fdmdv; +autotest; + +rand('state',1); +randn('state',1); + +% Core function for testing frequency offset estimator. Performs one test + +function sim_out = freq_off_est_test(sim_in) + global Nfilter; + global M; + + Rs = 50; + Nc = 4; + Nd = 2; + framesize = 32; + Fs = 8000; + Fcentre = 1500; + + afdmdv.Nsym = 2; + afdmdv.Nt = 3; + + afdmdv.Fs = 8000; + afdmdv.Nc = Nd*Nc-1; + afdmdv.Rs = Rs; + afdmdv.M = afdmdv.Fs/afdmdv.Rs; + afdmdv.tx_filter_memory = zeros(afdmdv.Nc+1, Nfilter); + afdmdv.Nfilter = Nfilter; + afdmdv.gt_alpha5_root = gen_rn_coeffs(0.5, 1/Fs, Rs, afdmdv.Nsym, afdmdv.M); + afdmdv.Fsep = 75; + afdmdv.phase_tx = ones(afdmdv.Nc+1,1); + freq_hz = afdmdv.Fsep*( -Nc*Nd/2 - 0.5 + (1:Nc*Nd) ); + afdmdv.freq_pol = 2*pi*freq_hz/Fs; + afdmdv.freq = exp(j*afdmdv.freq_pol); + afdmdv.Fcentre = 1500; + + afdmdv.fbb_rect = exp(j*2*pi*Fcentre/Fs); + afdmdv.fbb_phase_tx = 1; + afdmdv.fbb_phase_rx = 1; + afdmdv.phase_rx = ones(afdmdv.Nc+1,1); + + nin = M; + + P = afdmdv.P = 4; + afdmdv.Nfilter = afdmdv.Nsym*afdmdv.M; + afdmdv.rx_filter_mem_timing = zeros(afdmdv.Nc+1, afdmdv.Nt*afdmdv.P); + afdmdv.Nfiltertiming = afdmdv.M + afdmdv.Nfilter + afdmdv.M; + afdmdv.rx_filter_memory = zeros(afdmdv.Nc+1, afdmdv.Nfilter); + + acohpsk = standard_init(); + acohpsk.framesize = framesize; + acohpsk.ldpc_code = 0; + acohpsk.ldpc_code_rate = 1; + acohpsk.Nc = Nc; + acohpsk.Rs = Rs; + acohpsk.Ns = 4; + acohpsk.Nd = Nd; + acohpsk.modulation = 'qpsk'; + acohpsk.do_write_pilot_file = 0; + acohpsk = symbol_rate_init(acohpsk); + acohpsk.Ncm = 10*acohpsk.Nsymbrowpilot*M; + acohpsk.coarse_mem = zeros(1,acohpsk.Ncm); + acohpsk.Ndft = 2^(ceil(log2(acohpsk.Ncm))); + + frames = sim_in.frames; + EsNodB = sim_in.EsNodB; + foff = sim_in.foff; + dfoff = sim_in.dfoff; + fading_en = sim_in.fading_en; + + EsNo = 10^(EsNodB/10); + hf_delay_ms = 2; + phase_ch = 1; + + rand('state',1); + tx_bits_coh = round(rand(1,framesize*10)); + ptx_bits_coh = 1; + [spread spread_2ms hf_gain] = init_hf_model(Fs, Fs, frames*acohpsk.Nsymbrowpilot*afdmdv.M); + + hf_n = 1; + nhfdelay = floor(hf_delay_ms*Fs/1000); + ch_fdm_delay = zeros(1, acohpsk.Nsymbrowpilot*M + nhfdelay); + + sync = 0; next_sync = 1; + sync_start = 1; + freq_offset_log = []; + sync_time_log = []; + ch_fdm_frame_log = []; + ch_symb_log = []; + tx_fdm_frame_log = []; + + for f=1:frames + + acohpsk.frame = f; + + % + % Mod -------------------------------------------------------------------- + % + + tx_bits = tx_bits_coh(ptx_bits_coh:ptx_bits_coh+framesize-1); + ptx_bits_coh += framesize; + if ptx_bits_coh > length(tx_bits_coh) + ptx_bits_coh = 1; + end + + [tx_symb tx_bits] = bits_to_qpsk_symbols(acohpsk, tx_bits, [], []); + + tx_fdm_frame = []; + for r=1:acohpsk.Nsymbrowpilot + tx_onesymb = tx_symb(r,:); + [tx_baseband afdmdv] = tx_filter(afdmdv, tx_onesymb); + [tx_fdm afdmdv] = fdm_upconvert(afdmdv, tx_baseband); + tx_fdm_frame = [tx_fdm_frame tx_fdm]; + end + tx_fdm_frame_log = [tx_fdm_frame_log tx_fdm_frame]; + + % + % Channel -------------------------------------------------------------------- + % + + ch_fdm_frame = zeros(1,acohpsk.Nsymbrowpilot*M); + for i=1:acohpsk.Nsymbrowpilot*M + foff_rect = exp(j*2*pi*foff/Fs); + foff += dfoff; + phase_ch *= foff_rect; + ch_fdm_frame(i) = tx_fdm_frame(i) * phase_ch; + end + phase_ch /= abs(phase_ch); + + if fading_en + ch_fdm_delay(1:nhfdelay) = ch_fdm_delay(acohpsk.Nsymbrowpilot*M+1:nhfdelay+acohpsk.Nsymbrowpilot*M); + ch_fdm_delay(nhfdelay+1:nhfdelay+acohpsk.Nsymbrowpilot*M) = ch_fdm_frame; + + for i=1:acohpsk.Nsymbrowpilot*M + ahf_model = hf_gain*(spread(hf_n)*ch_fdm_frame(i) + spread_2ms(hf_n)*ch_fdm_delay(i)); + ch_fdm_frame(i) = ahf_model; + hf_n++; + end + end + + % each carrier has power = 2, total power 2Nc, total symbol rate NcRs, noise BW B=Fs + % Es/No = (C/Rs)/(N/B), N = var = 2NcFs/NcRs(Es/No) = 2Fs/Rs(Es/No) + + variance = 2*Fs/(acohpsk.Rs*EsNo); + uvnoise = sqrt(0.5)*(randn(1,acohpsk.Nsymbrowpilot*M) + j*randn(1,acohpsk.Nsymbrowpilot*M)); + noise = sqrt(variance)*uvnoise; + + ch_fdm_frame += noise; + ch_fdm_frame_log = [ch_fdm_frame_log ch_fdm_frame]; + + % + % Try to achieve sync -------------------------------------------------------------------- + % + + next_sync = sync; + + if sync == 0 + next_sync = 2; + acohpsk.f_est = Fcentre; + end + + [rx_fdm_frame_bb afdmdv.fbb_phase_rx] = freq_shift(ch_fdm_frame, -acohpsk.f_est, Fs, afdmdv.fbb_phase_rx); + + for r=1:acohpsk.Nsymbrowpilot + + % downconvert each FDM carrier to Nc separate baseband signals + + [rx_baseband afdmdv] = fdm_downconvert(afdmdv, rx_fdm_frame_bb(1+(r-1)*M:r*M), nin); + [rx_filt afdmdv] = rx_filter(afdmdv, rx_baseband, nin); + [rx_onesym rx_timing env afdmdv] = rx_est_timing(afdmdv, rx_filt, nin); + + ch_symb(r,:) = rx_onesym; + end + ch_symb_log = [ch_symb_log; ch_symb]; + + % coarse timing (frame sync) and initial fine freq est --------------------------------------------- + + [next_sync acohpsk] = frame_sync_fine_freq_est(acohpsk, ch_symb, sync, next_sync); + + % if we've acheived sync gather stats + + if (next_sync == 4) + freq_offset_log = [freq_offset_log acohpsk.f_fine_est+foff]; + sync_time_log = [sync_time_log f-sync_start]; + sync = 0; next_sync = 2; sync_start = f; + end + + %printf("f: %d sync: %d next_sync: %d\n", f, sync, next_sync); + [sync acohpsk] = sync_state_machine(acohpsk, sync, next_sync); + + end + + % ftx=fopen("coarse_tx.raw","wb"); fwrite(ftx, 1000*ch_fdm_frame_log, "short"); fclose(ftx); + + sim_out.freq_offset_log = freq_offset_log; + sim_out.sync_time_log = sync_time_log; + sim_out.ch_fdm_frame_log = ch_fdm_frame_log; + sim_out.ch_symb_log = ch_symb_log; + sim_out.tx_fdm_frame_log = tx_fdm_frame_log; +endfunction + + +function freq_off_est_test_single + sim_in.frames = 100; + sim_in.EsNodB = 20; + sim_in.foff = -15; + sim_in.dfoff = 0; + sim_in.fading_en = 1; + + sim_out = freq_off_est_test(sim_in); + + figure(1) + subplot(211) + plot(sim_out.freq_offset_log) + subplot(212) + hist(sim_out.freq_offset_log) + + figure(2) + subplot(211) + plot(sim_out.sync_time_log) + subplot(212) + hist(sim_out.sync_time_log) + + figure(3) + subplot(211) + plot(real(sim_out.tx_fdm_frame_log(1:2*960))) + subplot(212) + plot(real(sim_out.ch_symb_log(1:24,:)),'+') +endfunction + + +function [freq_off_log EsNodBSet] = freq_off_est_test_curves + EsNodBSet = [20 12 8]; + + sim_in.frames = 100; + sim_in.foff = -20; + sim_in.dfoff = 0; + sim_in.fading_en = 1; + freq_off_log = 1E6*ones(sim_in.frames, length(EsNodBSet) ); + sync_time_log = 1E6*ones(sim_in.frames, length(EsNodBSet) ); + + for i=1:length(EsNodBSet) + sim_in.EsNodB = EsNodBSet(i); + printf("%f\n", sim_in.EsNodB); + + sim_out = freq_off_est_test(sim_in); + freq_off_log(1:length(sim_out.freq_offset_log),i) = sim_out.freq_offset_log; + sync_time_log(1:length(sim_out.sync_time_log),i) = sim_out.sync_time_log; + end + + figure(1) + clf + hold on; + for i=1:length(EsNodBSet) + data = freq_off_log(find(freq_off_log(:,i) < 1E6),i); + s = std(data); + m = mean(data); + stdbar = [m-s; m+s]; + plot(EsNodBSet(i), data, '+') + plot([EsNodBSet(i) EsNodBSet(i)]+0.5, stdbar,'+-') + end + hold off + + axis([6 22 -25 25]) + if sim_in.fading_en + title_str = sprintf('foff = %d Hz Fading', sim_in.foff); + else + title_str = sprintf('foff = %d Hz AWGN', sim_in.foff); + end + title(title_str); + xlabel('Es/No (dB)') + ylabel('freq offset error (Hz)'); + + figure(2) + clf + hold on; + for i=1:length(EsNodBSet) + leg = sprintf("%d;%d dB;", i, EsNodBSet(i)); + plot(freq_off_log(find(freq_off_log(:,i) < 1E6),i),leg) + end + hold off; + title(title_str); + xlabel('test'); + ylabel('freq offset error (Hz)'); + legend('boxoff'); + + figure(3) + clf + hold on; + for i=1:length(EsNodBSet) + data = sync_time_log(find(sync_time_log(:,i) < 1E6),i); + if length(data) + s = std(data); + m = mean(data); + stdbar = [m-s; m+s]; + plot(EsNodBSet(i), data, '+') + plot([EsNodBSet(i) EsNodBSet(i)]+0.5, stdbar,'+-') + end + end + hold off; + axis([6 22 0 10]) + ylabel('sync time (frames)') + xlabel('Es/No (dB)'); + title(title_str); + + figure(4) + clf + hold on; + for i=1:length(EsNodBSet) + leg = sprintf("%d;%d dB;", i, EsNodBSet(i)); + plot(sync_time_log(find(sync_time_log(:,i) < 1E6),i),leg) + end + hold off; + title(title_str); + xlabel('test'); + ylabel('sync time (frames)'); + legend('boxoff'); + +endfunction + + +freq_off_est_test_single; +%freq_off_est_test_curves; + +% 1. start with +/- 20Hz offset +% 2. Measure frames to sync. How to define sync? Foff to withn 1 Hz. Sync state +% Need to see if we get false sync +% 3. Try shortened filter +% 4. Extend to parallel demods at +/- -- 2.25.1