#ifdef PC_BUILD #include "../gps-test-tool/main.h" #else #include "main.h" #endif #define KALMAN_Q 8.5e-6 #define KALMAN_R 4e-5 #define KALMAN_ERR_MAX 6e-4 __flash const char xml_header[] = "\n" "\n"; __flash const char xml_trk_start[] = "\t\n"; __flash const char xml_trkseg_end[] = "\t\t\n"; __flash const char xml_trkseg_start[] = "\t\t\n"; FIL gpx_file; static char buf[sizeof(xml_header)+sizeof(xml_trk_start)+2]; struct kalman_s { unsigned char initialized; float x_est_last; float P_last; float Q; float R; float K; }; #define PREV_POINTS_LENGTH 4 #define AVG_COUNT 3 #define MIN_DIST_DELTA 2.0 /* Elevation dead-band + smoothing, applied to the filtered ascent/descent totals */ #define ELEV_SMOOTH_TAU 30.0 /* seconds */ #define ELEV_DEADBAND 5.0 /* meters */ /* Large-jump detection thresholds. POS_JUMP_THRESHOLD is informational * logging only - the distance-diff spike check above already rejects * inconsistent position jumps regardless of size. ALT_JUMP_THRESHOLD both * logs and (combined with the same spike shape) rejects an inconsistent * altitude jump - unlike position there was previously no altitude * consistency check at all, letting a poor-geometry fix's bogus altitude * (right sats/HDOP, wrong solution - e.g. just after reacquiring signal) * straight into the accepted track. */ #define POS_JUMP_THRESHOLD 100.0 /* meters between consecutive raw fixes */ #define ALT_JUMP_THRESHOLD 50.0 /* meters between consecutive raw fixes */ /* Post-gap settling: a freshly reacquired fix's position/altitude solution * can keep drifting for several seconds after reacquisition even while * reported quality (satellite count, HDOP) already looks good throughout - * seen on a real recorded track exiting an underground space: HDOP was * already 0.81/17 sats just 7s after reacquiring, yet elevation kept * changing by several meters per sample for another ~10s after that before * settling to normal (sub-meter) noise. Sat count/HDOP were tried as a * "has it converged" signal and rejected: both fluctuate (sat count can * drop as easily as rise while the receiver's tracking set changes) in a * way that gave false "stable" readings well before the position had * actually settled. A fixed wall-clock window matches what the real data * showed converging by, without depending on a quality signal that proved * unreliable for this specific purpose in practice. * Withhold distance/elevation/GPX credit for points within this window * after a gap, while still updating the Kalman/last-position state normally * so a genuine tunnel-with-real-displacement case isn't penalized beyond * losing GAP_SETTLE_SECONDS worth of distance right at the exit. */ #define GAP_SETTLE_SECONDS 18 /* real data showed convergence taking ~15-17s after reacquisition */ struct prev_points_s { struct location_s data[PREV_POINTS_LENGTH]; unsigned char start; unsigned char count; }; struct avg_store_s { float lat; float lon; time_t time; }; struct elevation_s { unsigned char initialized; float smoothed_alt; float baseline; time_t last_time; }; static struct gpx_s { struct prev_points_s prev_points; unsigned char avg_count; unsigned char paused; unsigned char point_count; struct avg_store_s avg_store; struct location_s last_saved; struct location_s last_distance_point; /* Last accepted point for distance calculation */ struct kalman_s kalman[2]; struct elevation_s elevation; time_t gap_settle_until; /* points before this time are held back post-gap; 0 = not settling */ } gpx; float kalman_predict(struct kalman_s *k, float data); void kalman_init(struct kalman_s *k); float elevation_process(float alt, time_t time); void add_elevation_filtered(float amount, unsigned char is_gain); void prev_points_append(struct location_s *new){ gpx.prev_points.data[(gpx.prev_points.start + gpx.prev_points.count)%PREV_POINTS_LENGTH] = *new; if(++gpx.prev_points.count > PREV_POINTS_LENGTH){ gpx.prev_points.count--; gpx.prev_points.start++; gpx.prev_points.start %= PREV_POINTS_LENGTH; } } struct location_s *prev_points_get(unsigned char index){ unsigned char i, addr = gpx.prev_points.start; for(i=0; i tag */ strcpy_P(buf, xml_header); strcat_P(buf, xml_trk_start); return f_write(file, buf, strlen(buf), &bw); } void gpx_save_single_point(struct location_s *loc) { FIL gpx; UINT bw; unsigned char err = 0; char *time = get_iso_time(loc->time, 1); iso_time_to_filename(time); xsprintf(buf, PSTR("%s-POINT.GPX"), time); xprintf(PSTR("Writing single point in %s\r\n"), buf); if ((err = f_open(&gpx, buf, FA_WRITE | FA_OPEN_ALWAYS))) { f_close(&gpx); // System.status = STATUS_FILE_OPEN_ERROR; xputs_P(PSTR("File open error\r\n")); return; /* Failed to open file */ } strcpy_P(buf, xml_header); err |= f_write(&gpx, buf, strlen(buf), &bw); xsprintf(buf, PSTR("\t\n\n"), loc->lat, loc->lon); err |= f_write(&gpx, buf, strlen(buf), &bw); err |= f_close(&gpx); if (err) { /* TODO */ } } unsigned char is_paused(void) { return tracking_is_paused(); } unsigned char gpx_write(struct location_s *loc, FIL *file) { unsigned int bw; const char *time; unsigned char paused = tracking_is_paused(); if (paused) { if (!gpx.paused) { strcpy_P(buf, xml_trkseg_end); gpx.paused = 1; gpx.point_count = 0; } else { return 0; /* nothing to store */ } } else { if (gpx.paused) { strcpy_P(buf, xml_trkseg_start); f_write(file, buf, strlen(buf), &bw); gpx.paused = 0; } time = get_iso_time(loc->time, 0); xsprintf(buf, PSTR("\t\t\t\n\t\t\t\t%.2f\n\t\t\t\t\n"), loc->lat, loc->lon, loc->alt, time); strcat_P(buf, PSTR("\t\t\t\n")); } { unsigned char ret; unsigned long int t0 = get_uptime_ms(); ret = f_write(file, buf, strlen(buf), &bw); io_mark(t0); return ret; } } unsigned char gpx_close(FIL *file) { unsigned int bw; buf[0] = '\0'; if (!gpx.paused) strcpy_P(buf, xml_trkseg_end); strcat_P(buf, PSTR("\t\n\n")); f_write(file, buf, strlen(buf), &bw); return f_close(file); } /* Call once, whenever fix_gap_detected() reports a gap on the point about to * be passed to gpx_process_point() - see fix_gap_detected()'s comment for why * a gap (regardless of cause: fix loss, a 2D fix, poor HDOP, or the receiver * simply not producing lines) needs handling at all. Reset the Kalman/window * state that would otherwise splice the old and new fixes together, and * re-run the existing skip_points warm-up: a freshly reacquired fix can carry * a wildly wrong position/altitude for a few seconds while it converges (seen * after a real signal loss - a nominally-3D, in-range-HDOP fix still a couple * hundred meters off in altitude until more satellites came in), the same way * a fix has never been trusted right after power-on. */ void gpx_reset_gap(void) { if (get_flag(CONFFLAG_VERBOSE_LOG)) xputs_P(PSTR("FIX GAP: resetting warm-up/filter state\r\n")); kalman_init(&gpx.kalman[0]); kalman_init(&gpx.kalman[1]); gpx.prev_points.count = 0; gpx.avg_count = 0; gpx.avg_store.lat = 0; gpx.avg_store.lon = 0; gpx.avg_store.time = 0; gpx.point_count = 0; /* Without this, the first point(s) accepted after the gap (still inside * the prev_points warm-up above, so exempt from the distdiff/altdiff * spike checks until that window refills) would add one uncontrolled * distance/elevation segment straight from the pre-gap position - which * may be stale by an arbitrary, unknown amount of real movement during * the gap - with none of the usual spike protection. Clearing it forces * distance/elevation accumulation to (re)start fresh from the first * post-gap point actually reached once warm-up completes, the same way * it starts fresh after power-on. */ gpx.last_distance_point.lat = NAN; /* -1 (all-ones): "settling, but the deadline isn't computed yet" - the * real deadline needs the first post-gap point's own timestamp, which * isn't available here (this runs before gpx_process_point() sees that * point), so gpx_process_point() computes gap_settle_until = that * point's time + GAP_SETTLE_SECONDS on its first point after this. */ gpx.gap_settle_until = (time_t)-1; } void gpx_process_point(struct location_s *loc, FIL *file){ float lon_est, lon_err, lat_est, lat_err, dist; struct location_s *ptr; static struct location_s nloc; static struct location_s filtered_loc; if (tracking_is_paused()) System.points_paused++; if (gpx.point_count < System.conf.skip_points) { /* Skipping initial points */ gpx.point_count++; System.points_skipped++; return; } if (gpx.gap_settle_until) { if (gpx.gap_settle_until == (time_t)-1) gpx.gap_settle_until = loc->time + GAP_SETTLE_SECONDS; if (loc->time < gpx.gap_settle_until) { if (get_flag(CONFFLAG_VERBOSE_LOG)) xputs_P(PSTR("GAP SETTLE: waiting for fix to stabilize\r\n")); System.points_skipped++; return; } gpx.gap_settle_until = 0; /* window elapsed - resume normal processing below */ } /* Always apply Kalman filtering for distance calculation */ lat_est = kalman_predict(&gpx.kalman[0], loc->lat); lon_est = kalman_predict(&gpx.kalman[1], loc->lon); filtered_loc.lat = lat_est; filtered_loc.lon = lon_est; filtered_loc.time = loc->time; filtered_loc.alt = loc->alt; if (get_flag(CONFFLAG_DISABLE_FILTERS)) { /* Write unfiltered data to GPX, but always calculate distance/elevation from filtered data */ if (get_flag(CONFFLAG_VERBOSE_LOG)) xputs_P(PSTR("Write with filters disabled\r\n")); gpx_write(loc, file); System.points_written++; System.points_accepted++; /* Calculate distance and elevation from filtered points */ if (!isnan(gpx.last_distance_point.lat)) { float ele_change; dist = distance(&gpx.last_distance_point, &filtered_loc); if (dist > POS_JUMP_THRESHOLD) log_jump(0, dist); add_distance(dist); ele_change = filtered_loc.alt - gpx.last_distance_point.alt; if (fabs(ele_change) > ALT_JUMP_THRESHOLD) log_jump(1, ele_change); add_elevation(ele_change); } elevation_process(filtered_loc.alt, filtered_loc.time); gpx.last_distance_point = filtered_loc; } else { /* Apply Kalman error check. KALMAN_ERR_MAX is a threshold in degrees * of latitude - a degree of longitude covers less ground the further * from the equator, by a factor of cos(latitude), so lon_err is * scaled up here to the equivalent latitude-degree error before the * same threshold is applied to both. Without this, the effective * east-west tolerance in meters would shrink at higher latitudes * (e.g. roughly half as permissive at 60 degrees as at the equator), * causing spurious rejects that have nothing to do with fix quality. */ lat_err = fabs(loc->lat - lat_est); lon_err = fabs(loc->lon - lon_est) * cosf(loc->lat * M_PI / 180.0); // xprintf(PSTR("lat_err: %e, lon_err: %e, limit: %e\r\n"), lat_err, lon_err, (float)KALMAN_ERR_MAX); if(lat_err > KALMAN_ERR_MAX || lon_err > KALMAN_ERR_MAX){ if (get_flag(CONFFLAG_VERBOSE_LOG)) xputs_P(PSTR("KALMAN REJECT\r\n")); log_reject(REJECT_REASON_KALMAN); return; } prev_points_append(&filtered_loc); if(gpx.prev_points.count == PREV_POINTS_LENGTH){ float dist12 = distance(prev_points_get(0), prev_points_get(1)); float dist34 = distance(prev_points_get(2), prev_points_get(3)); float dist32 = distance(prev_points_get(2), prev_points_get(1)); /* Same spike test as the distance check below, applied to altitude: * a poor-geometry fix (e.g. right after reacquiring signal) can pass * the satellite-count/HDOP gate in nmea.c yet still carry a wildly * wrong altitude, so catch it here the same way a position spike is * caught - a jump surrounded by two much smaller ones on both sides. * TODO: GSA's VDOP would be a more direct altitude-quality gate than * this shape heuristic, but VDOP is only in GSA, which is parsed * after GGA within the epoch (see nmea.c's gp_gga_parse HDOP * comment) - using it here would need buffering a point across the * epoch boundary until VDOP arrives. Left as a future refinement. */ float alt12 = fabs(prev_points_get(0)->alt - prev_points_get(1)->alt); float alt34 = fabs(prev_points_get(2)->alt - prev_points_get(3)->alt); float alt32 = fabs(prev_points_get(2)->alt - prev_points_get(1)->alt); if (get_flag(CONFFLAG_VERBOSE_LOG)) xprintf(PSTR("New distance: %fm\r\n"), dist32); /* A real spike must stand out on BOTH sides - bigger than the * older segment AND bigger than the newer one. The previous * "dist34 > dist12 && dist32 > dist12" shape only checked dist32 * against dist12 twice from different angles, so it stayed true * for every step of a sustained, monotonic acceleration (dist12 < * dist32 < dist34, completely normal e.g. when picking up pace * while hiking) - rejecting real point after real point until the * trend broke, then crediting all of that accumulated real * distance to whichever single point finally got through, as one * artificially fast-looking GPX segment. Comparing dist32 against * both neighbors instead lets a genuine acceleration/deceleration * ramp through cleanly while still catching an isolated outlier * surrounded by two much smaller segments. */ if(dist32 > dist12 && dist32 > dist34){ if (get_flag(CONFFLAG_VERBOSE_LOG)) xputs_P(PSTR("DISTANCE DIFF REJECT\r\n")); log_reject(REJECT_REASON_DISTDIFF); return; } if(alt32 > ALT_JUMP_THRESHOLD && alt32 > alt12 && alt32 > alt34){ if (get_flag(CONFFLAG_VERBOSE_LOG)) xputs_P(PSTR("ALTITUDE DIFF REJECT\r\n")); log_reject(REJECT_REASON_ALTDIFF); return; } if (dist32 > POS_JUMP_THRESHOLD) log_jump(0, dist32); ptr = prev_points_get(PREV_POINTS_LENGTH - 2); } else { if(gpx.prev_points.count >= PREV_POINTS_LENGTH-2){ ptr = prev_points_get(gpx.prev_points.count - 2); if (get_flag(CONFFLAG_VERBOSE_LOG)) xputs_P(PSTR("NEW\r\n")); } else { return; } } if(distance(&gpx.last_saved, ptr) < MIN_DIST_DELTA){ if (get_flag(CONFFLAG_VERBOSE_LOG)) xputs_P(PSTR("Too small position change REJECT\r\n")); log_reject(REJECT_REASON_MINDIST); return; } if (get_flag(CONFFLAG_VERBOSE_LOG)) xputs_P(PSTR("ACCEPT\r\n")); log_reject_flush(); /* close out any pending reject burst now that good data has resumed */ System.points_accepted++; /* Calculate distance and elevation for accepted point */ if (!isnan(gpx.last_distance_point.lat)) { float ele_change; dist = distance(&gpx.last_distance_point, ptr); add_distance(dist); ele_change = ptr->alt - gpx.last_distance_point.alt; if (fabs(ele_change) > ALT_JUMP_THRESHOLD) log_jump(1, ele_change); add_elevation(ele_change); } elevation_process(ptr->alt, ptr->time); gpx.last_distance_point = *ptr; gpx.avg_store.lat += ptr->lat; gpx.avg_store.lon += ptr->lon; /* Use the newest accepted point's time, not e.g. the midpoint's: the * averaged lat/lon reflects real movement up through this instant, so * timestamping it any earlier claims that movement happened faster * than it did. This matters most exactly when the accepted points * feeding this average were NOT evenly spaced in real time (e.g. a * reject streak between two of them) - an earlier timestamp there * would understate the elapsed time for the position delta and * inflate the speed an external tool derives from consecutive GPX * points. */ gpx.avg_store.time = ptr->time; if(++gpx.avg_count == AVG_COUNT){ nloc.lat = gpx.avg_store.lat / AVG_COUNT; nloc.lon = gpx.avg_store.lon / AVG_COUNT; nloc.time = gpx.avg_store.time; nloc.alt = gpx.elevation.smoothed_alt; /* filtered (smoothed) altitude, for the filtered GPX write */ gpx.avg_count = 0; gpx.avg_store.lat = 0; gpx.avg_store.lon = 0; gpx.avg_store.time = 0; gpx.last_saved = nloc; gpx_write(&nloc, file); System.points_written++; } } if (System.time_start == 0) System.time_start = utc; } void kalman_init(struct kalman_s *k){ k->initialized = 0; k->P_last = 0; //the noise in the system k->Q = KALMAN_Q; // process variance k->R = KALMAN_R; // measurement variance k->K = 0; } float kalman_predict(struct kalman_s *k, float data){ if(!k->initialized){ //initial values for the kalman filter k->x_est_last = data; k->initialized = 1; return data; } //do a prediction float x_temp_est = k->x_est_last; float P_temp = k->P_last + k->Q; //calculate the Kalman gain k->K = P_temp * (1.0/(P_temp + k->R)); //correct float x_est = x_temp_est + k->K * (data - x_temp_est); k->P_last = (1 - k->K) * P_temp; k->x_est_last = x_est; return x_est; } #define R_EARTH 6371e3 // m float distance(struct location_s *pos1, struct location_s *pos2){ float lat1 = pos1->lat * M_PI / 180.0; float lat2 = pos2->lat * M_PI / 180.0; float dlat = (pos2->lat - pos1->lat) * M_PI / 180.0; float dlon = (pos2->lon - pos1->lon) * M_PI / 180.0; float a = sinf(dlat/2.0) * sinf(dlat/2.0) + cosf(lat1) * cosf(lat2) * sinf(dlon/2.0) * sinf(dlon/2.0); float c = 2 * atan2f(sqrtf(a), sqrtf(1-a)); float ret = R_EARTH * c; return ret; } /* Initial bearing (degrees, 0-360, 0=north) from pos1 to pos2 - used by the * auto-pause resume-by-distance check (autopause.c) to tell a real, directed * departure from a smooth-but-directionless position drift (e.g. a satellite * geometry change reshaping the fix while the receiver doesn't move at all). */ float bearing(struct location_s *pos1, struct location_s *pos2){ float lat1 = pos1->lat * M_PI / 180.0; float lat2 = pos2->lat * M_PI / 180.0; float dlon = (pos2->lon - pos1->lon) * M_PI / 180.0; float y = sinf(dlon) * cosf(lat2); float x = cosf(lat1) * sinf(lat2) - sinf(lat1) * cosf(lat2) * cosf(dlon); float deg = atan2f(y, x) * 180.0 / M_PI; if (deg < 0) deg += 360.0; return deg; } void add_distance(float dist) { unsigned char paused = tracking_is_paused(); if (!paused) System.distance += (dist+0.005)*100.0; if (get_flag(CONFFLAG_VERBOSE_LOG)) xprintf(PSTR("Distance: %.2f m; sum: %.2f m\r\n"), (double)dist, (double)System.distance/100.0); } /* Unfiltered (raw), per-point gain/loss - kept only for comparison against the * filtered (dead-band+smoothed) totals in the periodic status line/session summary. */ void add_elevation(float ele_change) { unsigned char paused = tracking_is_paused(); if (!paused) { if (ele_change > 0) { System.elevation_gain_raw += (ele_change+0.05)*10.0; } else if (ele_change < 0) { System.elevation_loss_raw += (-ele_change+0.05)*10.0; } } if (get_flag(CONFFLAG_VERBOSE_LOG)) xprintf(PSTR("Elevation change: %.1f m; raw gain: %.1f m, raw loss: %.1f m\r\n"), (double)ele_change, (double)System.elevation_gain_raw/10.0, (double)System.elevation_loss_raw/10.0); } /* Filtered (dead-band + smoothed) gain/loss, used for display, GPX ascent stats * and the periodic status line/session summary. */ void add_elevation_filtered(float amount, unsigned char is_gain) { unsigned long int dm = (unsigned long int)(amount*10.0 + 0.5); if (is_gain) System.elevation_gain += dm; else System.elevation_loss += dm; } /* Exponential smoothing (~30s time constant) followed by a 5m dead-band on the * result, so that a step in the smoothed altitude only counts once it clears * the dead-band, and only the amount past the dead-band edge is credited. */ float elevation_process(float alt, time_t time) { unsigned char paused = tracking_is_paused(); unsigned long int alt_dm; float delta; if (!gpx.elevation.initialized) { gpx.elevation.smoothed_alt = alt; gpx.elevation.baseline = alt; gpx.elevation.last_time = time; gpx.elevation.initialized = 1; } else { float dt = (float)(time - gpx.elevation.last_time); float alpha; if (dt <= 0) dt = 1.0; gpx.elevation.last_time = time; alpha = dt / (ELEV_SMOOTH_TAU + dt); gpx.elevation.smoothed_alt += alpha * (alt - gpx.elevation.smoothed_alt); } if (gpx.elevation.smoothed_alt > 0) { alt_dm = (unsigned long int)(gpx.elevation.smoothed_alt*10.0 + 0.5); if (alt_dm > System.alt_max) System.alt_max = alt_dm; } if (paused) { /* Discard the dead-band reference drift accumulated while paused, * the same way distance/raw elevation data is discarded when paused. */ gpx.elevation.baseline = gpx.elevation.smoothed_alt; } else { delta = gpx.elevation.smoothed_alt - gpx.elevation.baseline; if (delta > ELEV_DEADBAND) { add_elevation_filtered(delta - ELEV_DEADBAND, 1); gpx.elevation.baseline = gpx.elevation.smoothed_alt - ELEV_DEADBAND; } else if (delta < -ELEV_DEADBAND) { add_elevation_filtered(-delta - ELEV_DEADBAND, 0); gpx.elevation.baseline = gpx.elevation.smoothed_alt + ELEV_DEADBAND; } } return gpx.elevation.smoothed_alt; }