Changeset 7065


Ignore:
Timestamp:
Nov 28, 2014 11:37:32 AM (5 years ago)
Author:
toivo
Message:

Add a function for interpolation using Neville's Algorithm to jmi_delay.c and use it for linear interpolation, for #4097.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/RuntimeLibrary/src/jmi/jmi_delay.c

    r7054 r7065  
    505505}
    506506
     507/* Interpolate the minimum degree polynomial through the given points using Neville's Algorithm */
     508#define MAX_NEVILLE_PTS 16
     509static jmi_real_t neville_evaluate(jmi_delaybuffer_t *buffer, jmi_real_t t, int first_index, int n_points) {
     510    jmi_real_t work[MAX_NEVILLE_PTS];
     511    jmi_real_t ts[MAX_NEVILLE_PTS];
     512    jmi_delay_point_t *buf = buffer->buf;
     513    int i, n;
     514    if (n_points > MAX_NEVILLE_PTS) return -1; /* todo: error */
     515
     516    /* Copy the initial points */
     517    for (i=0; i < n_points; i++) {
     518        int pos = index2pos(buffer, i + first_index);
     519        work[i] = buf[pos].y;
     520        ts[i] = buf[pos].t;
     521    }
     522
     523    /* Evaluate the intermediate results */
     524    /* Loop over the polynomial order n */
     525    for (n=1; n < n_points; n++) {
     526        /* Loop over the interpolating polynomials of order n */
     527        /* work[i] contains p[i:i+n-1](t) before this, and p[i:i+n](t) after */
     528        for (i=0; i < n_points-n; i++) {
     529            int j = i+n;
     530            work[i] = ((ts[j] - t)*work[i] + (t-ts[i])*work[i+1])/(ts[j]-ts[i]);
     531        }
     532    }
     533    return work[0];
     534}
     535
    507536static jmi_real_t evaluate(jmi_delaybuffer_t *buffer, jmi_boolean at_event,
    508537                           jmi_real_t tr, jmi_delay_position_t *position) {
     
    514543    if (update_position(buffer, at_event, tr, position) < 0) return -1; /* todo: error */
    515544
    516     /* Linear interpolation */
     545    /* Clamp to initial value if before initial time. Todo: better way? */
     546   
     547    if ((position->curr_interval == buffer->head_index) && (tr < buf[buffer->head].t)) {
     548        return buf[buffer->head].y;
     549    }
     550   
     551
     552    return neville_evaluate(buffer, tr, position->curr_interval, 2);
     553
     554    /* Linear interpolation. Todo: remove / use as special case? */
     555    /*
    517556    {
    518557        int lpos = index2pos(buffer, position->curr_interval);
     
    526565        return y0 + (y1-y0)*(tr-t0)/(t1-t0);
    527566    }
     567    */
    528568}
    529569
Note: See TracChangeset for help on using the changeset viewer.