From d588e8394258cf536fc0947ad3e4ceeaf33e216e Mon Sep 17 00:00:00 2001 From: drowe67 Date: Fri, 22 May 2015 02:14:04 +0000 Subject: [PATCH] timing offset simulation and correction by Octave demod. +/- 1000ppm works well for AWGN and fading channels git-svn-id: https://svn.code.sf.net/p/freetel/code@2145 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/cohpsk.m | 145 +++--------------------------------- codec2-dev/octave/tcohpsk.m | 94 +++++++++++++++-------- 2 files changed, 75 insertions(+), 164 deletions(-) diff --git a/codec2-dev/octave/cohpsk.m b/codec2-dev/octave/cohpsk.m index 868fa0d9..b222ba4d 100644 --- a/codec2-dev/octave/cohpsk.m +++ b/codec2-dev/octave/cohpsk.m @@ -64,7 +64,7 @@ function sim_in = symbol_rate_init(sim_in) sim_in.Nsymbrow = Nsymbrow = Nsymb/Nc; sim_in.Npilotsframe = Npilotsframe = 2; sim_in.Nsymbrowpilot = Nsymbrowpilot = Nsymbrow + Npilotsframe; - + if verbose == 2 printf("Each frame contains %d data bits or %d data symbols, transmitted as %d symbols by %d carriers.", framesize, Nsymb, Nsymbrow, Nc); printf(" There are %d pilot symbols in each carrier together at the start of each frame, then %d data symbols.", Npilotsframe, Ns); @@ -459,147 +459,20 @@ function test_bits_coh_file(test_bits_coh) endfunction -% Frequency offset estimation -------------------------------------------------- - -% This function was used in initial Nov 2014 experiments - -function [f_max s_max] = freq_off_est(rx_fdm, tx_pilot, offset, n) - - Fs = 8000; - nc = 1800; % portion we wish to correlate over (first 2 rows on pilots) - - % downconvert to complex baseband to remove images - - f = 1500; - 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 - - -% Set of functions to implement latest and greatest freq offset -% estimation, March 2015 ---------------------- - -% 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; - 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 - h = 0.5 - 0.5*cos(2*pi*(0:Ncm-1)/(Ncm-1)); - T = abs(fft(coarse_mem .* h, Ndft)).^2; - sc = Ndft/Fs; - - 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); - - 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 - - function [ch_symb rx_timing rx_filt rx_baseband afdmdv f_est] = rate_Fs_rx_processing(afdmdv, ch_fdm_frame, f_est, nsymb, nin, freq_track) M = afdmdv.M; - + rx_baseband = []; rx_filt = []; rx_timing = []; + ch_fdm_frame_index = 1; + for r=1:nsymb % shift signal to nominal baseband, this will put Nc/2 carriers either side of 0 Hz - [rx_fdm_frame_bb afdmdv.fbb_phase_rx] = freq_shift(ch_fdm_frame(1+(r-1)*M:r*M), -f_est, afdmdv.Fs, afdmdv.fbb_phase_rx); + [rx_fdm_frame_bb afdmdv.fbb_phase_rx] = freq_shift(ch_fdm_frame(ch_fdm_frame_index:ch_fdm_frame_index + nin - 1), -f_est, afdmdv.Fs, afdmdv.fbb_phase_rx); + ch_fdm_frame_index += nin; % downconvert each FDM carrier to Nc separate baseband signals @@ -613,6 +486,12 @@ function [ch_symb rx_timing rx_filt rx_baseband afdmdv f_est] = rate_Fs_rx_proce ch_symb(r,:) = rx_onesym; + % we only allow a timing shift on one symbol per frame + + if nin != M + nin = M; + end + % freq tracking, see test_ftrack.m for unit test. Placed in % this function as it needs to work on a symbol by symbol basis % rather than frame by frame. This means the control loop diff --git a/codec2-dev/octave/tcohpsk.m b/codec2-dev/octave/tcohpsk.m index 3abc785b..713c8ec7 100644 --- a/codec2-dev/octave/tcohpsk.m +++ b/codec2-dev/octave/tcohpsk.m @@ -1,23 +1,26 @@ % tcohpsk.m % David Rowe Oct 2014 % -% Octave script that: +% Octave coherent PSK modem script that hs two modes: % % 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. +% the reference versions of the same modem written in Octave. % -% or (ii) can optionally be used to run an Octave version of the cohpsk -% modem to tune and develop it. + +% (ii) Runs the Octave version of the cohpsk modem to tune and develop +% it, including extensive channel simulations such as AWGN noise, +% fading/HF, frequency offset, frequency drift, and tx/rx sample +% rate differences. % TODO: % -% [ ] Test +% [X] Test % [X] AWGN channel % [X] freq offset % [X] fading channel % [X] freq drift -% [ ] timing drift +% [X] timing drift % [X] tune perf/impl loss to get closer to ideal % [X] linear interp of phase for better fading perf % [X] freq offset/drift feedback loop @@ -28,10 +31,13 @@ % [X] ability to "unsync" when signal dissapears % [ ] some calibrated tests against FreeDV 1600 % + compare sound quality at various Es/Nos +% [ ] sync +% + set some req & implement % [ ] way to handle eom w/o nasties % + like mute ouput when signal has gone or v low snr % + instantaneous snr % [ ] nasty rig filter passband +% [ ] pilot based EsNo estimation graphics_toolkit ("gnuplot"); more off; @@ -64,18 +70,20 @@ if strcmp(test, 'compare to c') fading_en = 0; hf_delay_ms = 2; compare_with_c = 1; + sample_rate_ppm = 0; end % should be BER around 0.015 to 0.02 if strcmp(test, 'awgn') frames = 100; - foff = 0; - dfoff = 0; + foff = 58.8; + dfoff = -0.5/Fs; EsNodB = 8; fading_en = 0; hf_delay_ms = 2; compare_with_c = 0; + sample_rate_ppm = 1000; end % Similar to AWGN - should be BER around 0.015 to 0.02 @@ -83,11 +91,12 @@ end if strcmp(test, 'fading'); frames = 100; foff = -53.1; - dfoff = 0.5/Fs; + dfoff = 0.0/Fs; EsNodB = 12; fading_en = 1; hf_delay_ms = 2; compare_with_c = 0; + sample_rate_ppm = -100; end EsNo = 10^(EsNodB/10); @@ -225,8 +234,9 @@ ch_fdm_delay = zeros(1, acohpsk.Nsymbrowpilot*M + nhfdelay); % main loop -------------------------------------------------------------------- +% run mod and channel as aseparate loop so we can resample to simulate sample rate differences + for f=1:frames - acohpsk.frame = f; 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) @@ -256,9 +266,6 @@ for f=1:frames % Channel -------------------------------------------------------------------- % - %[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); @@ -292,6 +299,35 @@ for f=1:frames ch_fdm_frame += noise; ch_fdm_frame_log = [ch_fdm_frame_log ch_fdm_frame]; +end + +% simulate difference in sample clocks + +ch_fdm_frame_log = resample(ch_fdm_frame_log, (1E6 + sample_rate_ppm), 1E6); + +% Now run demod ---------------------------------------------------------------- + +ch_fdm_frame_log_index = 1; +f = 0; + +while (ch_fdm_frame_log_index + acohpsk.Nsymbrowpilot*M+M/P) < length(ch_fdm_frame_log) + acohpsk.frame = f++; + + if sync == 0 + nin = M; + else + nin = M; + if rx_timing(1) > M/P + nin = M + M/P; + end + if rx_timing(1) < -M/P + nin = M - M/P; + end + end + + nin_frame = (acohpsk.Nsymbrowpilot-1)*M + nin; + ch_fdm_frame = ch_fdm_frame_log(ch_fdm_frame_log_index:ch_fdm_frame_log_index + nin_frame - 1); + ch_fdm_frame_log_index += nin_frame; % % Demod ---------------------------------------------------------------------- @@ -299,11 +335,10 @@ for f=1:frames % store two frames of received samples so we can rewind if we get a good candidate - ch_fdm_frame_buf(1:(Nsw-1)*acohpsk.Nsymbrowpilot*M) = ch_fdm_frame_buf(acohpsk.Nsymbrowpilot*M+1:Nsw*acohpsk.Nsymbrowpilot*M); - ch_fdm_frame_buf((Nsw-1)*acohpsk.Nsymbrowpilot*M+1:Nsw*acohpsk.Nsymbrowpilot*M) = ch_fdm_frame; + ch_fdm_frame_buf(1:Nsw*acohpsk.Nsymbrowpilot*M-nin_frame) = ch_fdm_frame_buf(nin_frame+1:Nsw*acohpsk.Nsymbrowpilot*M); + ch_fdm_frame_buf(Nsw*acohpsk.Nsymbrowpilot*M-nin_frame+1:Nsw*acohpsk.Nsymbrowpilot*M) = ch_fdm_frame; next_sync = sync; - nin = M; % if out of sync do Initial Freq offset estimation over NSW frames to flush out memories @@ -407,8 +442,7 @@ for f=1:frames % BER stats if f > 2 - error_positions = xor(prev_tx_bits2, rx_bits); - %error_positions = xor(prev_tx_bits, rx_bits); + error_positions = xor(tx_bits_log((f-3)*framesize+1:(f-2)*framesize), rx_bits); Nerrs += sum(error_positions); nerr_log = [nerr_log sum(error_positions)]; Tbits += length(error_positions); @@ -416,22 +450,14 @@ for f=1:frames printf("\r [%d]", f); end - - %rx_fdm_frame_bb_log = [rx_fdm_frame_bb_log rx_fdm_frame_bb]; - %rx_baseband_log = [rx_baseband_log rx_baseband]; - %rx_filt_log = [rx_filt_log rx_filt]; - %rx_timing_log = [rx_timing_log rx_timing]; - %ch_symb_log = [ch_symb_log; ch_symb]; - % ct_symb_ff_log = [ct_symb_ff_log; acohpsk.ct_symb_ff_buf(1:acohpsk.Nsymbrowpilot,:)]; + % reset BER stats if we lose sync - - if sync == 0 - Nerrs = 0; - Tbits = 0; - nerr_log = []; + if sync == 1 + %Nerrs = 0; + %Tbits = 0; + %nerr_log = []; end - % printf("f: %d sync: %d next_sync: %d\n", f, sync, next_sync); [sync acohpsk] = sync_state_machine(acohpsk, sync, next_sync); prev_tx_bits2 = prev_tx_bits; @@ -563,6 +589,12 @@ else plot(foff_log(1:length(f_est_log)) - f_est_log + Fcentre) title('freq offset estimation error'); + figure(8) + clf + subplot(211) + stem(real(rx_filt_log(1,52:4:200))) + subplot(212) + stem(imag(rx_filt_log(1,52:4:200))) end -- 2.25.1