%run_curves
%run_curves_estimators
%acquisition_histograms(0, 0)
-acquisition_test(10, -3, 5)
+acquisition_test(10, 4, 5)
+ 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
%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)
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
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
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----------------------------------
% 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
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;
Nbitsperframe = states.Nbitsperframe;
Nrowsperframe = states.Nrowsperframe;
Nsamperframe = states.Nsamperframe;
+timing_mx_thresh = 0.3;
W = states.W;
w = states.w;
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);
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 = [];
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];
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]);
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])
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');
*
* 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) {
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;
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 */
/* 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
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;
}
/*
}
}
+ 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;
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 */
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 */
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 */
}
#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;
int verbose;
int sample_point;
int timing_est;
+ int timing_valid;
+ float timing_mx;
float coarse_foff_est_hz;
int nin;
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];
#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()
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);
tout += 1;
tin += 1.0 + sample_rate_ppm / 1E6;
}
-
+ //printf("n: %d tout: %d tin: %f\n", n, tout, tin);
+
return tout;
}
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];
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;
int lnew;
COMP rxbuf_in[max_samples_per_frame];
-
+
for (i=0; i<nin; i++,prx++) {
ofdm->rxbuf[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; f<NFRAMES; f++) {
/* For initial testng, timing est is off, so nin is always
fixed. TODO: we need a constant for rxbuf_in[] size that
}
assert(prx <= max_samples_per_frame*NFRAMES);
+ #ifdef TESTING_FILE
+ fread(rx_scaled, sizeof(short), nin, fin);
+ for(i=0; i<nin; i++) {
+ rx[i].real = (float)rx_scaled[i]/ASCALE;
+ rx[i].imag = 0.0;
+ }
+ ofdm_demod(ofdm, rx_bits, rx);
+ #else
ofdm_demod(ofdm, rx_bits, rxbuf_in);
-
+ #endif
+
+ int Nerrs = 0;
+ for(i=0; i<OFDM_BITSPERFRAME; i++) {
+ if (test_bits_ofdm[i] != rx_bits[i]) {
+ Nerrs++;
+ }
+ }
+ //printf("f: %d Nerr: %d\n", f, Nerrs);
+
/* rx vector logging -----------------------------------*/
assert(nin_tot < samples_per_frame*NFRAMES);
for(i=0; i<OFDM_RXBUF; i++) {
rxbuf_log[OFDM_RXBUF*f+i].real = crealf(ofdm->rxbuf[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++) {
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 */
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);
projects / freetel-svn-tracking.git / commitdiff