nicholas@2642: // Taken for JS-Xtract
nicholas@2642: // https://github.com/nickjillings/js-xtract
nicholas@2642: 
nicholas@2642: function xtract_resample(data, p, q, n) {
nicholas@2642:     // Same function call as matlab:
nicholas@2642:     // data is the array
nicholas@2642:     // p is the target sample rate
nicholas@2642:     // q is the source sample rate
nicholas@2642:     // n is the desired filter order, n==1024 if undefined
nicholas@2642: 
nicholas@2642:     function filter(N, c) {
nicholas@2642:         var c_b, Re, Im, b;
nicholas@2642:         c_b = Math.floor(c * N);
nicholas@2642:         Re = new Float64Array(N);
nicholas@2642:         Im = new Float64Array(N);
nicholas@2642:         var i, j;
nicholas@2642:         for (i = 0; i < c_b; i++) {
nicholas@2642:             Re[i] = 1;
nicholas@2642:         }
nicholas@2642:         for (i = N - c_b + 1; i < N; i++) {
nicholas@2642:             Re[i] = 1;
nicholas@2642:         }
nicholas@2642:         inverseTransform(Re, Im);
nicholas@2642:         // Scale and shift into Im
nicholas@2642:         for (i = 0; i < N; i++) {
nicholas@2642:             j = (i + (N >> 1)) % N;
nicholas@2642:             Im[i] = Re[j] / N;
nicholas@2642:             // Apply compute blackman-harris to Im
nicholas@2642:             var rad = (Math.PI * i) / (N);
nicholas@2642:             Im[i] *= 0.35875 - 0.48829 * Math.cos(2 * rad) + 0.14128 * Math.cos(4 * rad) - 0.01168 * Math.cos(6 * rad);
nicholas@2642:         }
nicholas@2642:         return Im;
nicholas@2642:     }
nicholas@2642: 
nicholas@2642:     function polyn(data, K) {
nicholas@2642:         var N = data.length;
nicholas@2642:         var x = [0, data[0], data[1]];
nicholas@2642:         var dst = new Float64Array(K);
nicholas@2642:         var ratio = K / N;
nicholas@2642:         var tinc = 1 / ratio;
nicholas@2642:         var n = 0,
nicholas@2642:             t = 0,
nicholas@2642:             k;
nicholas@2642:         for (k = 0; k < K; k++) {
nicholas@2642:             if (t == n) {
nicholas@2642:                 // Points lie on same time
nicholas@2642:                 dst[k] = data[n];
nicholas@2642:             } else {
nicholas@2642:                 var y = (t - n - 1) * (t - n) * x[0] - 2 * (t - n - 1) * (t - n + 1) * x[1] + (t - n) * (t - n + 1) * x[2];
nicholas@2642:                 dst[k] = y / 2;
nicholas@2642:             }
nicholas@2642:             t += tinc;
nicholas@2642:             if (t >= n + 1) {
nicholas@2642:                 n = Math.floor(t);
nicholas@2642:                 x[0] = data[n - 1];
nicholas@2642:                 x[1] = data[n];
nicholas@2642:                 if (n + 1 < N) {
nicholas@2642:                     x[2] = data[n + 1];
nicholas@2642:                 } else {
nicholas@2642:                     x[2] = 0.0;
nicholas@2642:                 }
nicholas@2642:             }
nicholas@2642:         }
nicholas@2642:         return dst;
nicholas@2642:     }
nicholas@2642: 
nicholas@2642:     function zp(a) {
nicholas@2642:         var b = new Float64Array(a.length * 2);
nicholas@2642:         for (var n = 0; n < a.length; n++) {
nicholas@2642:             b[n] = a[n];
nicholas@2642:         }
nicholas@2642:         return b;
nicholas@2642:     }
nicholas@2642: 
nicholas@2642:     function overlap(X, b) {
nicholas@2642:         var i, f;
nicholas@2642:         var Y = new Float64Array(X.length);
nicholas@2642:         var N = b.length;
nicholas@2642:         var N2 = 2 * N;
nicholas@2642:         var B = {
nicholas@2642:             real: zp(b),
nicholas@2642:             imag: new Float64Array(N * 2)
nicholas@2642:         }
nicholas@2642:         transform(B.real, B.imag);
nicholas@2642:         var Xi = X.xtract_get_data_frames(N, N, false);
nicholas@2642:         var Yi = Y.xtract_get_data_frames(N, N, false);
nicholas@2642:         var x_last = new Float64Array(N);
nicholas@2642:         var y_last = new Float64Array(N);
nicholas@2642:         var w = new Float64Array(N2);
nicholas@2642:         for (i = 0; i < N2; i++) {
nicholas@2642:             var rad = (Math.PI * i) / (N2);
nicholas@2642:             w[i] = 0.35875 - 0.48829 * Math.cos(2 * rad) + 0.14128 * Math.cos(4 * rad) - 0.01168 * Math.cos(6 * rad);
nicholas@2642:         }
nicholas@2642:         var xF = {
nicholas@2642:             real: new Float64Array(N2),
nicholas@2642:             imag: new Float64Array(N2)
nicholas@2642:         }
nicholas@2642:         var yF = {
nicholas@2642:             real: new Float64Array(N2),
nicholas@2642:             imag: new Float64Array(N2)
nicholas@2642:         }
nicholas@2642:         for (f = 0; f < Xi.length; f++) {
nicholas@2642:             for (i = 0; i < N; i++) {
nicholas@2642:                 xF.real[i] = x_last[i] * w[i];
nicholas@2642:                 xF.real[i + N] = Xi[f][i] * w[i + N];
nicholas@2642:                 x_last[i] = Xi[f][i];
nicholas@2642:                 xF.imag[i] = 0;
nicholas@2642:                 xF.imag[i + N] = 0;
nicholas@2642:             }
nicholas@2642:             transform(xF.real, xF.imag);
nicholas@2642:             for (i = 0; i < N2; i++) {
nicholas@2642:                 yF.real[i] = xF.real[i] * B.real[i] - xF.imag[i] * B.imag[i];
nicholas@2642:                 yF.imag[i] = xF.imag[i] * B.real[i] + xF.real[i] * B.imag[i];
nicholas@2642:             }
nicholas@2642:             transform(yF.imag, yF.real);
nicholas@2642:             // Perform fft_shift and scale
nicholas@2642:             for (i = 0; i < N; i++) {
nicholas@2642:                 var h = yF.real[i + N] / N;
nicholas@2642:                 yF.real[i + N] = yF.real[i] / N;
nicholas@2642:                 yF.real[i] = h;
nicholas@2642:             }
nicholas@2642:             for (i = 0; i < N; i++) {
nicholas@2642:                 Yi[f][i] = (yF.real[i] + y_last[i]);
nicholas@2642:                 y_last[i] = yF.real[i + N];
nicholas@2642:             }
nicholas@2642:         }
nicholas@2642:         return Y;
nicholas@2642:     }
nicholas@2642: 
nicholas@2642:     // Determine which way to go
nicholas@2642:     var b, N = data.length;
nicholas@2642:     if (typeof n != "number" || n <= 0) {
nicholas@2642:         n = 512;
nicholas@2642:     }
nicholas@2642:     if (p == q) {
nicholas@2642:         return data;
nicholas@2642:     }
nicholas@2642:     var ratio = (p / q);
nicholas@2642:     var K = Math.floor(N * ratio);
nicholas@2642:     var dst;
nicholas@2642:     if (p > q) {
nicholas@2642:         // Upsampling
nicholas@2642:         // 1. Expand Data range
nicholas@2642:         dst = polyn(data, K);
nicholas@2642:         // 2. Filter out spurious energy above q
nicholas@2642:         var b = filter(n, 1 / ratio);
nicholas@2642:         overlap(data, b);
nicholas@2642:     } else {
nicholas@2642:         // Downsampling
nicholas@2642:         // 1. Filter out energy above p
nicholas@2642:         var b = filter(n, ratio / 2);
nicholas@2642:         var ds1 = overlap(data, b);
nicholas@2642:         // 2. Decrease data range
nicholas@2642:         dst = polyn(ds1, K);
nicholas@2642:     }
nicholas@2642:     return dst;
nicholas@2642: }