Trigonometric Interpolation returns NaN

Question

I'm a musician, who's new to programming. I use JavaScript inside Max Msp (hence the bang() and post() functions) to create a trigonometric interpolation, interpolating between given equidistant points (for testing, only values of sine from [0, 2π) and returning values from the same points). When I run the code, it returns NaN, except for x = 0, as my tau() function returns only 1 in this special case. Could it be, that it has something to do with summing Math.sin results?

var f = new Array(9);
var TWO_PI = 2*Math.PI;

bang();

function bang() {
    for(var i = 0; i < f.length; i++) {
    f[i] = Math.sin(i/f.length*TWO_PI);
    //post("f[" + i + "]: " + Math.round(f[i]*1000)/1000 + "\n");
    }
    
    var points = new Array(f.length);
    for(var i = 0; i < points.length; i++) {
        var idx = i/points.length*TWO_PI;
        points[i] = [i, p(idx)];
    //post("p[" + points[i][0] + "]: " + Math.round(points[i][1]*1000)/1000 + "\n");
    }
    
    console.log("p(2): " + p(2/points.length*TWO_PI) + "\n");
}

function p(x) {
  var result = 0;
  for(var k = 0; k < f.length; k++) {
    result += f[k]*tau(k, x);
  }
  return result;
}

function tau(k, x) {
    var dividend = sinc(1/2*f.length*(x-k/f.length*TWO_PI));
    var divisor = sinc(1/2*(x-k/f.length*TWO_PI));
    var result = dividend/divisor;
    if(f.length%2 == 0) result *= Math.cos(1/2*(x-k/f.length*TWO_PI));
    
    if(x == 0) return 1;
    return result;
}

function sinc(x) {
    return Math.sin(x)/x;
}

Answer 1

In your tau function, if x equals k / f.length * TWO_PI (which it will since x is multiples of 1 / points.length * TWO_PI) your sinc function divides by 0, making divisor equal to NaN, which then propagates.

Answer 2

You have to be a bit careful in implementing sinc to avoid dividing by 0. One way is to say that if x is small enough we can replace sin(x) by the first few terms of its taylor series, and all the terms are divisible by x.

I don't know javascript but here is the function in C in case it is of use

#define SINC_EPS (1e-6)
// for small x, 
// missing sinc terms start with pow(x,4)/120, and value close to 1
// so the error too small to be seen in a double
double  sinc( double x)
{       if ( fabs(x) < SINC_EPS)
        { return 1.0 - x*x/6.0;
        }
        else 
        { return sin(x)/x;
        }
}