From: drowe67 Date: Fri, 23 Mar 2018 21:03:07 +0000 (+0000) Subject: coarse timing with normalised max passes tofdm X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=fd09172e2928a7ba5d9e613a59422c9e5f479fb5;p=freetel-svn-tracking.git coarse timing with normalised max passes tofdm git-svn-id: https://svn.code.sf.net/p/freetel/code@3427 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/octave/ofdm_dev.m b/codec2-dev/octave/ofdm_dev.m index 520eb1fc..aaeba9e0 100644 --- a/codec2-dev/octave/ofdm_dev.m +++ b/codec2-dev/octave/ofdm_dev.m @@ -1011,4 +1011,4 @@ init_cml('/home/david/Desktop/cml/'); %run_curves %run_curves_estimators %acquisition_histograms(0, 0) -acquisition_test(10, -3, 5) +acquisition_test(10, 4, 5) diff --git a/codec2-dev/octave/ofdm_lib.m b/codec2-dev/octave/ofdm_lib.m index bd0303b9..72ebeea8 100644 --- a/codec2-dev/octave/ofdm_lib.m +++ b/codec2-dev/octave/ofdm_lib.m @@ -63,25 +63,32 @@ endfunction + more suitable for real time implementation #} -function [t_est foff_est] = coarse_sync(states, rx, rate_fs_pilot_samples) - Nsamperframe = states.Nsamperframe; Fs = states.Fs; +function [t_est foff_est timing_valid timing_mx] = coarse_sync(states, rx, rate_fs_pilot_samples) + ofdm_load_const; Npsam = length(rate_fs_pilot_samples); - verbose = states.verbose; Ncorr = length(rx) - (Nsamperframe+Npsam) + 1; assert(Ncorr > 0); - corr1 = corr2 = zeros(1,Ncorr); - av_level = Npsam*sqrt(rx*rx')/Ncorr + 1E-12; + corr = zeros(1,Ncorr); + + % normalise correlation so we can compare to a threshold across varying input levels + + av_level = 2*sqrt(states.timing_norm*(rx*rx')/length(rx)) + 1E-12; + + % correlate with pilots at start and end of frame to determine timing offset + for i=1:Ncorr rx1 = rx(i:i+Npsam-1); rx2 = rx(i+Nsamperframe:i+Nsamperframe+Npsam-1); - corr1(i) = rx1 * rate_fs_pilot_samples'; - corr2(i) = rx2 * rate_fs_pilot_samples'; + corr(i) = abs(rx1 * rate_fs_pilot_samples' + rx2 * rate_fs_pilot_samples')/av_level; end - corr = (abs(corr1) + abs(corr2))/av_level; - [mx t_est] = max(corr); - printf(" max: %f t_est: %d\n", mx, t_est); + [timing_mx t_est] = max(corr); + timing_valid = timing_mx > timing_mx_thresh; + if verbose > 1 + printf(" max: %f timing_est: %d timing_valid: %d\n", timing_mx, timing_est, timing_valid); + end + #{ % original freq offset est code that never made it into C. Have some concerns about CPU % load of performing FFT, althout a smaller one could have been used with interpolation @@ -113,7 +120,7 @@ function [t_est foff_est] = coarse_sync(states, rx, rate_fs_pilot_samples) %printf("t_est: %d\n", t_est); figure(7); clf; plot(abs(corr)) - axis([1 Ncorr 0 2]) + %axis([1 Ncorr 0 2]) #{ figure(8) plot(C) @@ -126,6 +133,7 @@ function [t_est foff_est] = coarse_sync(states, rx, rate_fs_pilot_samples) axis([-0.2 0.2 -0.2 0.2]) %hold on; plot(rx(t_est+Nsamperframe:t_est+Npsam+Nsamperframe-1) .* rate_fs_pilot_samples','g+'); hold off; end + endfunction @@ -199,6 +207,12 @@ function states = ofdm_init(bps, Rs, Tcp, Ns, Nc) rate_fs_pilot_samples = states.pilots * W/states.M; states.rate_fs_pilot_samples = [rate_fs_pilot_samples(states.M-states.Ncp+1:states.M) rate_fs_pilot_samples]; + + % pre-compute a constant used in coarse_sync() + + Npsam = length(states.rate_fs_pilot_samples); + states.timing_norm = Npsam*(states.rate_fs_pilot_samples * states.rate_fs_pilot_samples'); + % printf("timing_norm: %f\n", states.timing_norm) % LDPC code is optionally enabled @@ -316,18 +330,22 @@ function [rx_bits states aphase_est_pilot_log rx_np rx_amp] = ofdm_demod(states, st = M+Ncp + Nsamperframe + 1 - floor(ftwindow_width/2) + (timing_est-1); en = st + Nsamperframe-1 + M+Ncp + ftwindow_width-1; - [ft_est coarse_foff_est_hz] = coarse_sync(states, rxbuf(st:en) .* exp(-j*woff_est*(st:en)), rate_fs_pilot_samples); - timing_est = timing_est + ft_est - ceil(ftwindow_width/2); + [ft_est coarse_foff_est_hz timing_valid timing_mx] = coarse_sync(states, rxbuf(st:en) .* exp(-j*woff_est*(st:en)), rate_fs_pilot_samples); + + if timing_valid + timing_est = timing_est + ft_est - ceil(ftwindow_width/2); + + % Black magic to keep sample_point inside cyclic prefix. Or something like that. + delta_t = ft_est - ceil(ftwindow_width/2); + sample_point = max(timing_est+Ncp/4, sample_point); + sample_point = min(timing_est+Ncp, sample_point); + end + if verbose > 1 - printf(" ft_est: %2d timing_est: %2d sample_point: %2d\n", ft_est, timing_est, sample_point); + printf(" ft_est: %2d timing_est: %2d mx: %3.2f sample_point: %2d\n", ft_est, timing_est, timing_mx, sample_point); end - % Black magic to keep sample_point inside cyclic prefix. Or something like that. - - delta_t = ft_est - ceil(ftwindow_width/2); - sample_point = max(timing_est+Ncp/4, sample_point); - sample_point = min(timing_est+Ncp, sample_point); end % down convert at current timing instant---------------------------------- @@ -431,7 +449,7 @@ function [rx_bits states aphase_est_pilot_log rx_np rx_amp] = ofdm_demod(states, % Adjust nin to take care of sample clock offset nin = Nsamperframe; - if timing_en + if timing_en && timing_valid thresh = (M+Ncp)/8; tshift = (M+Ncp)/4; if timing_est > thresh @@ -449,6 +467,8 @@ function [rx_bits states aphase_est_pilot_log rx_np rx_amp] = ofdm_demod(states, states.rx_sym = rx_sym; states.rxbuf = rxbuf; states.nin = nin; + states.timing_valid = timing_valid; + states.timing_mx = timing_mx; states.timing_est = timing_est; states.sample_point = sample_point; states.delta_t = delta_t; diff --git a/codec2-dev/octave/ofdm_load_const.m b/codec2-dev/octave/ofdm_load_const.m index db66d6d9..70459ddb 100644 --- a/codec2-dev/octave/ofdm_load_const.m +++ b/codec2-dev/octave/ofdm_load_const.m @@ -12,6 +12,7 @@ bps = states.bps; Nbitsperframe = states.Nbitsperframe; Nrowsperframe = states.Nrowsperframe; Nsamperframe = states.Nsamperframe; +timing_mx_thresh = 0.3; W = states.W; w = states.w; diff --git a/codec2-dev/octave/tofdm.m b/codec2-dev/octave/tofdm.m index f7f15381..0e02baf0 100644 --- a/codec2-dev/octave/tofdm.m +++ b/codec2-dev/octave/tofdm.m @@ -19,6 +19,7 @@ autotest; Ts = 0.018; Tcp = 0.002; Rs = 1/Ts; bps = 2; Nc = 16; Ns = 8; states = ofdm_init(bps, Rs, Tcp, Ns, Nc); +states.verbose = 0; ofdm_load_const; rand('seed',1); @@ -48,7 +49,8 @@ states.rxbuf(Nrxbuf-nin+1:Nrxbuf) = rx_log(prx:nin); prx += nin; rxbuf_log = []; rxbuf_in_log = []; rx_sym_log = []; foff_hz_log = []; -timing_est_log = []; coarse_foff_est_hz_log = []; sample_point_log = []; +timing_est_log = timing_valid_log = timing_mx_log = []; +coarse_foff_est_hz_log = []; sample_point_log = []; phase_est_pilot_log = []; rx_amp_log = []; rx_np_log = []; rx_bits_log = []; @@ -87,6 +89,8 @@ for f=1:Nframes rx_amp_log = [rx_amp_log arx_amp]; foff_hz_log = [foff_hz_log; states.foff_est_hz]; timing_est_log = [timing_est_log; states.timing_est]; + timing_valid_log = [timing_valid_log; states.timing_valid]; + timing_mx_log = [timing_mx_log; states.timing_mx]; coarse_foff_est_hz_log = [coarse_foff_est_hz_log; states.coarse_foff_est_hz]; sample_point_log = [sample_point_log; states.sample_point]; rx_np_log = [rx_np_log arx_np]; @@ -109,6 +113,7 @@ system(path_to_tofdm); load tofdm_out.txt; fg = 1; + figure(fg++); clf; plot(rx_np_log,'+'); title('Octave Scatter Diagram'); axis([-1.5 1.5 -1.5 1.5]); figure(fg++); clf; plot(rx_np_log_c,'+'); title('C Scatter Diagram'); axis([-1.5 1.5 -1.5 1.5]); @@ -136,7 +141,9 @@ d = phase_est_pilot_log - phase_est_pilot_log_c; d = angle(exp(j*d)); stem_sig_and_error(fg, 211, phase_est_pilot_log_c, d, 'phase est pilot', [1 length(phase_est_pilot_log_c) -1.5 1.5]) stem_sig_and_error(fg++, 212, rx_amp_log_c, rx_amp_log - rx_amp_log_c, 'rx amp', [1 length(rx_amp_log_c) -1.5 1.5]) -stem_sig_and_error(fg++, 111, foff_hz_log_c, (foff_hz_log - foff_hz_log_c), 'foff hz', [1 length(foff_hz_log_c) -1.5 1.5]) +stem_sig_and_error(fg , 211, foff_hz_log_c, (foff_hz_log - foff_hz_log_c), 'foff hz', [1 length(foff_hz_log_c) -1.5 1.5]) + +stem_sig_and_error(fg++, 212, timing_mx_log_c, (timing_mx_log - timing_mx_log_c), 'timing mx', [1 length(timing_mx_log_c) 0 2]) stem_sig_and_error(fg, 211, timing_est_log_c, (timing_est_log - timing_est_log_c), 'timing est', [1 length(timing_est_log_c) -1.5 1.5]) stem_sig_and_error(fg++, 212, sample_point_log_c, (sample_point_log - sample_point_log_c), 'sample point', [1 length(sample_point_log_c) -1.5 1.5]) @@ -155,6 +162,8 @@ check(rx_sym_log, rx_sym_log_c, 'rx_sym', tol=5E-3); check(phase_est_pilot_log, phase_est_pilot_log_c, 'phase_est_pilot', tol=2E-3, its_an_angle=1); check(rx_amp_log, rx_amp_log_c, 'rx_amp'); check(timing_est_log, timing_est_log_c, 'timing_est'); +check(timing_valid_log, timing_valid_log_c, 'timing_valid'); +check(timing_mx_log, timing_mx_log_c, 'timing_mx'); check(coarse_foff_est_hz_log, coarse_foff_est_hz_log_c, 'coarse_foff_est_hz'); check(sample_point_log, sample_point_log_c, 'sample_point'); check(foff_hz_log, foff_hz_log_c, 'foff_est_hz'); diff --git a/codec2-dev/src/ofdm.c b/codec2-dev/src/ofdm.c index 5bf716de..bc87943a 100644 --- a/codec2-dev/src/ofdm.c +++ b/codec2-dev/src/ofdm.c @@ -149,6 +149,8 @@ static complex float vector_sum(complex float *a, int num_elements) { * * Also estimates coarse frequency offset based on sampling the pilots * phase at half symbol intervales. + * + * Unlike Octave version use states to return a few values. */ static int coarse_sync(struct OFDM *ofdm, complex float *rx, int length, float *foff_est) { @@ -156,9 +158,16 @@ static int coarse_sync(struct OFDM *ofdm, complex float *rx, int length, float * int Ncorr = length - (OFDM_SAMPLESPERFRAME + (OFDM_M + OFDM_NCP)); int Fs = OFDM_FS; int SFrame = OFDM_SAMPLESPERFRAME; - float corr[Ncorr]; + float corr[Ncorr], av_level; int i, j, k; + complex float acc = 0.0f + 0.0f * I; + for (i = 0; i < length; i++) { + acc += rx[i] * conjf(rx[i]); + } + + av_level = 2.0*sqrt(ofdm->timing_norm*crealf(acc)/length) + 1E-12; + for (i = 0; i < Ncorr; i++) { complex float temp = 0.0f + 0.0f * I; @@ -168,21 +177,27 @@ static int coarse_sync(struct OFDM *ofdm, complex float *rx, int length, float * temp = temp + (rx[i + j + SFrame] * csam); } - corr[i] = cabsf(temp); + corr[i] = cabsf(temp)/av_level; } /* find the max magnitude and its index */ - float mag = 0.0f; - int t_est = 0; + float timing_mx = 0.0f; + int timing_est = 0; for (i = 0; i < Ncorr; i++) { - if (corr[i] > mag) { - mag = corr[i]; - t_est = i; + if (corr[i] > timing_mx) { + timing_mx = corr[i]; + timing_est = i; } } + ofdm->timing_mx = timing_mx; + ofdm->timing_valid = timing_mx > OFDM_TIMING_MX_THRESH; + if (ofdm->verbose > 1) { + fprintf(stderr, " max: %f timing_est: %d timing_valid: %d\n", ofdm->timing_mx, timing_est, ofdm->timing_valid); + } + /* Coarse frequency estimation - note this isn't always used, some CPU could be saved buy putting this in another function */ @@ -197,13 +212,13 @@ static int coarse_sync(struct OFDM *ofdm, complex float *rx, int length, float * /* pilot at start of frame */ - p1 = p1 + (rx[t_est + j] * csam1); - p2 = p2 + (rx[t_est + k] * csam2); + p1 = p1 + (rx[timing_est + j] * csam1); + p2 = p2 + (rx[timing_est + k] * csam2); /* pilot at end of frame */ - p3 = p3 + (rx[t_est + j + SFrame] * csam1); - p4 = p4 + (rx[t_est + k + SFrame] * csam2); + p3 = p3 + (rx[timing_est + j + SFrame] * csam1); + p4 = p4 + (rx[timing_est + k + SFrame] * csam2); } /* Calculate sample rate of phase samples, we are sampling phase @@ -217,8 +232,8 @@ static int coarse_sync(struct OFDM *ofdm, complex float *rx, int length, float * with 0 inputs. */ *foff_est = Fs1 * cargf(conjf(p1)*p2 + conjf(p3)*p4 + 1E-12)/(2.0*M_PI); - - return t_est; + + return timing_est; } /* @@ -330,6 +345,10 @@ struct OFDM *ofdm_create(const struct OFDM_CONFIG * config) { } } + for (i = 0; i < OFDM_RXBUF; i++) { + ofdm->rxbuf[i] = 0.0f + 0.0f * I; + } + for (i = 0; i < (OFDM_NS + 3); i++) { for (j = 0; j < (OFDM_NC + 2); j++) { ofdm->rx_sym[i][j] = 0.0f + 0.0f * I; @@ -347,6 +366,8 @@ struct OFDM *ofdm_create(const struct OFDM_CONFIG * config) { ofdm->foff_est_hz = 0.0f; ofdm->sample_point = 0; ofdm->timing_est = 0; + ofdm->timing_valid = 0; + ofdm->timing_mx = 0.0f; ofdm->nin = OFDM_SAMPLESPERFRAME; /* create the OFDM waveform */ @@ -364,6 +385,8 @@ struct OFDM *ofdm_create(const struct OFDM_CONFIG * config) { for (i = 0, j = (OFDM_M - OFDM_NCP); i < OFDM_NCP; i++, j++) { ofdm->pilot_samples[i] = temp[j]; } + + // states.timing_norm = Npsam*(rate_fs_pilot_samples*rate_fs_pilot_samples'); /* Now copy the whole thing after the above */ @@ -371,6 +394,15 @@ struct OFDM *ofdm_create(const struct OFDM_CONFIG * config) { ofdm->pilot_samples[i] = temp[j]; } + /* calculate constant used to normalise timing correlation maximum */ + + complex float acc = 0.0f + 0.0f * I; + for (i = 0; i < OFDM_M+OFDM_NCP; i++) { + acc += ofdm->pilot_samples[i] * conjf(ofdm->pilot_samples[i]); + } + ofdm->timing_norm = (OFDM_M + OFDM_NCP) * crealf(acc); + //fprintf(stderr, "timing_norm: %f\n", ofdm->timing_norm); + return ofdm; /* Success */ } diff --git a/codec2-dev/src/ofdm_internal.h b/codec2-dev/src/ofdm_internal.h index 2aa3b9d8..162cf334 100644 --- a/codec2-dev/src/ofdm_internal.h +++ b/codec2-dev/src/ofdm_internal.h @@ -76,7 +76,8 @@ extern "C" { #define OFDM_MAX_SAMPLESPERFRAME (OFDM_SAMPLESPERFRAME + (OFDM_M + OFDM_NCP)/4) #define OFDM_RXBUF (3 * OFDM_SAMPLESPERFRAME + 3 * (OFDM_M + OFDM_NCP)) - +#define OFDM_TIMING_MX_THRESH 0.3 + /* Dummy struct for now, will contain constant configuration for OFDM modem */ struct OFDM_CONFIG{ int a; @@ -91,6 +92,8 @@ struct OFDM { int verbose; int sample_point; int timing_est; + int timing_valid; + float timing_mx; float coarse_foff_est_hz; int nin; @@ -99,6 +102,7 @@ struct OFDM { bool phase_est_en; complex float pilot_samples[OFDM_M + OFDM_NCP]; + float timing_norm; complex float W[OFDM_NC + 2][OFDM_M]; complex float rxbuf[OFDM_RXBUF]; complex float pilots[OFDM_NC + 2]; diff --git a/codec2-dev/unittest/tofdm.c b/codec2-dev/unittest/tofdm.c index b790292a..ce97a9aa 100644 --- a/codec2-dev/unittest/tofdm.c +++ b/codec2-dev/unittest/tofdm.c @@ -45,6 +45,8 @@ #define SAMPLE_CLOCK_OFFSET_PPM 100 #define FOFF_HZ 0.5f +#define ASCALE (2E5*1.1491/2.0) /* scale from shorts back to floats */ + /*---------------------------------------------------------------------------*\ FUNCTION....: fs_offset() @@ -63,9 +65,10 @@ static int fs_offset(COMP out[], COMP in[], int n, float sample_rate_ppm) { double tin = 0.0; int tout = 0; - while (tin < (double) n) { + while (tin < (double) (n-1)) { t1 = (int) floor(tin); t2 = (int) ceil(tin); + assert(t2 < n); f = (tin - (double) t1); @@ -75,7 +78,8 @@ static int fs_offset(COMP out[], COMP in[], int n, float sample_rate_ppm) { tout += 1; tin += 1.0 + sample_rate_ppm / 1E6; } - + //printf("n: %d tout: %d tin: %f\n", n, tout, tin); + return tout; } @@ -139,6 +143,8 @@ int main(int argc, char *argv[]) float foff_hz_log[NFRAMES]; int rx_bits_log[OFDM_BITSPERFRAME*NFRAMES]; int timing_est_log[NFRAMES]; + int timing_valid_log[NFRAMES]; + float timing_mx_log[NFRAMES]; float coarse_foff_est_hz_log[NFRAMES]; int sample_point_log[NFRAMES]; @@ -181,7 +187,7 @@ int main(int argc, char *argv[]) Demod \*---------------------------------------------------------*/ - /* Init rx with ideal timing so we can test with timing estimation disabled */ + /* Init/pre-load rx with ideal timing so we can test with timing estimation disabled */ int Nsam = samples_per_frame*NFRAMES; int prx = 0; @@ -189,20 +195,29 @@ int main(int argc, char *argv[]) int lnew; COMP rxbuf_in[max_samples_per_frame]; - + for (i=0; irxbuf[OFDM_RXBUF-nin+i] = rx_log[prx].real + I*rx_log[prx].imag; } - + int nin_tot = 0; /* disable estimators for initial testing */ - ofdm_set_verbose(ofdm, true); + ofdm_set_verbose(ofdm, false); ofdm_set_timing_enable(ofdm, true); ofdm_set_foff_est_enable(ofdm, true); ofdm_set_phase_est_enable(ofdm, true); + #ifdef TESTING_FILE + FILE *fin=fopen("/home/david/codec2-dev/build_linux/src/ofdm_c_test.raw", "rb"); + assert(fin != NULL); + int Nbitsperframe = ofdm_get_bits_per_frame(ofdm); + int Nmaxsamperframe = ofdm_get_max_samples_per_frame(); + short rx_scaled[Nmaxsamperframe]; + COMP rx[Nmaxsamperframe]; + #endif + for(f=0; frxbuf[i]); rxbuf_log[OFDM_RXBUF*f+i].imag = cimagf(ofdm->rxbuf[i]); - } + } for (i = 0; i < (OFDM_NS + 3); i++) { for (j = 0; j < (OFDM_NC + 2); j++) { @@ -271,6 +303,8 @@ int main(int argc, char *argv[]) foff_hz_log[f] = ofdm->foff_est_hz; timing_est_log[f] = ofdm->timing_est + 1; /* offset by 1 to match Octave */ + timing_valid_log[f] = ofdm->timing_valid; /* offset by 1 to match Octave */ + timing_mx_log[f] = ofdm->timing_mx; /* offset by 1 to match Octave */ coarse_foff_est_hz_log[f] = ofdm->coarse_foff_est_hz; sample_point_log[f] = ofdm->sample_point + 1; /* offset by 1 to match Octave */ @@ -296,6 +330,8 @@ int main(int argc, char *argv[]) octave_save_float(fout, "rx_amp_log_c", (float*)rx_amp_log, 1, OFDM_ROWSPERFRAME*OFDM_NC*NFRAMES, OFDM_ROWSPERFRAME*OFDM_NC*NFRAMES); octave_save_float(fout, "foff_hz_log_c", foff_hz_log, NFRAMES, 1, 1); octave_save_int(fout, "timing_est_log_c", timing_est_log, NFRAMES, 1); + octave_save_int(fout, "timing_valid_log_c", timing_valid_log, NFRAMES, 1); + octave_save_float(fout, "timing_mx_log_c", timing_mx_log, NFRAMES, 1, 1); octave_save_float(fout, "coarse_foff_est_hz_log_c", coarse_foff_est_hz_log, NFRAMES, 1, 1); octave_save_int(fout, "sample_point_log_c", sample_point_log, NFRAMES, 1); octave_save_complex(fout, "rx_np_log_c", (COMP*)rx_np_log, 1, OFDM_ROWSPERFRAME*OFDM_NC*NFRAMES, OFDM_ROWSPERFRAME*OFDM_NC*NFRAMES);