This specification defines an interface that provides signal processing and mathematical functions for 32-bit floating point [[!TYPED-ARRAYS]] (Float32Array).

Introduction

This specification defines methods that operate on entire arrays rather than single elements at a time. The methods are useful for implementing performance critical array processing, such as digital signal processing.

A fundamental principle of all methods defined in this specification is that they operate on existing arrays rather than creating new ones. Thus, it is possible to use the methods in applications that are sensitive to heap allocation and garbage collection.

See for some examples of using the API.

Definitions

Terms

The generic term TypedArray is used to indicate any valid typed array view type.

The term NaN is short for Not a Number, and is a numeric value representing an undefined or unrepresentable value.

The term RMS is short for Root Mean Square, which is a measure of signal energy. For a given signal it is calculated as the square root of the arithmetic mean of the squares of the absolute values of the signal.

The term single precision means binary 32-bit IEEE 754 floating point.

The term double precision means binary 64-bit IEEE 754 floating point.

Pseudo code functions

Much of this specification is written in terms of mathematical expressions. In such expressions, the following pseudo function definitions MUST apply:

clamp(x, xmin, xmax)
  • Returns x if xmin ≤ x ≤ xmax.
  • Returns xmin if x < xmin.
  • Returns xmax if x > xmax.
  • Returns NaN if any of x, xmin or xmax is NaN.
modulo(x, y)

Returns x - floor(x/y) * y.

sign(x)
  • Returns -1 if x < 0 or x = -0.
  • Returns 1 if x > 0 or x = +0.
  • Returns NaN if x is NaN.

Furthermore, each of the following pseudo functions MUST use the same definition as the Math object method with the same name, as specified in [[!ECMA-262]], with the exception that single precision arithmetic SHOULD be used instead of double precision arithmetic: abs, acos, asin, atan, atan2, ceil, cos, exp, floor, log, min, max, pow, random, round, sin, sqrt, tan.

Complex signals

Complex valued signals are represented as two separate arrays: one containing the real part of the signal (suffixed Real in the IDL), and one containing the imaginary part of the signal (suffixed Imag in the IDL).

Numerical accuracy

Unless otherwise noted, numerical computations that are defined in this specification MUST meet the following requirements:

The rationale is that computational performance is important. Conforming clients should be able to utilize as many optimization techniques as possible when implementing the API, including (but not limited to) using hardware specific instructions and data types that may or may not fully support IEEE 754.

As an example of reordering of operations, a conforming client may choose to implement sqrt(x) as x * (1 / sqrt(x)).

Overlapping operation

Here we should specify approximately:

Overlapping operations are not allowed except for the special case where the destination array is the same as one of the source arrays (may differ for different methods).

The rationale is that while an implementation could allow overlapping operations for all methods, it would come at a performance and/or complexity cost.

ArrayMath interface

void add (Float32Array dst, Float32Array x, Float32Array y)
For each k in the range 0 ≤ k < min(dst.length, x.length, y.length), the method MUST compute dst[k] = x[k] + y[k].
void add (Float32Array dst, double x, Float32Array y)
For each k in the range 0 ≤ k < min(dst.length, y.length), the method MUST compute dst[k] = x + y[k].
void sub (Float32Array dst, Float32Array x, Float32Array y)
For each k in the range 0 ≤ k < min(dst.length, x.length, y.length), the method MUST compute dst[k] = x[k] - y[k].
void sub (Float32Array dst, double x, Float32Array y)
For each k in the range 0 ≤ k < min(dst.length, y.length), the method MUST compute dst[k] = x - y[k].
void mul (Float32Array dst, Float32Array x, Float32Array y)
For each k in the range 0 ≤ k < min(dst.length, x.length, y.length), the method MUST compute dst[k] = x[k] * y[k].
void mul (Float32Array dst, double x, Float32Array y)
For each k in the range 0 ≤ k < min(dst.length, y.length), the method MUST compute dst[k] = x * y[k].
void mulCplx (Float32Array dstReal, Float32Array dstImag, Float32Array xReal, Float32Array xImag, Float32Array yReal, Float32Array yImag)
This method performs multiplication of two complex arrays.

For each k in the range 0 ≤ k < min(dstReal.length, dstImag.length, xReal.length, xImag.length, yReal.length, yImag.length), the method MUST compute: 

dstReal[k] = xReal[k] * yReal[k] - xImag[k] * yImag[k]
dstImag[k] = xReal[k] * yImag[k] + xImag[k] * yReal[k]
void mulCplx (Float32Array dstReal, Float32Array dstImag, double xReal, double xImag, Float32Array yReal, Float32Array yImag)
This method performs multiplication of a complex scalar and a complex array.

For each k in the range 0 ≤ k < min(dstReal.length, dstImag.length, yReal.length, yImag.length), the method MUST compute: 

dstReal[k] = xReal * yReal[k] - xImag * yImag[k]
dstImag[k] = xReal * yImag[k] + xImag * yReal[k]
void div (Float32Array dst, Float32Array x, Float32Array y)
For each k in the range 0 ≤ k < min(dst.length, x.length, y.length), the method MUST compute dst[k] = x[k] / y[k].
void div (Float32Array dst, double x, Float32Array y)
For each k in the range 0 ≤ k < min(dst.length, y.length), the method MUST compute dst[k] = x / y[k].
void divCplx (Float32Array dstReal, Float32Array dstImag, Float32Array xReal, Float32Array xImag, Float32Array yReal, Float32Array yImag)
This method performs division of two complex arrays.

For each k in the range 0 ≤ k < min(dstReal.length, dstImag.length, xReal.length, xImag.length, yReal.length, yImag.length), the method MUST compute: 

denom = yReal[k] * yReal[k] + yImag[k] * yImag[k]
dstReal[k] = (xReal[k] * yReal[k] + xImag[k] * yImag[k]) / denom
dstImag[k] = (xImag[k] * yReal[k] - xReal[k] * yImag[k]) / denom
void divCplx (Float32Array dstReal, Float32Array dstImag, double xReal, double xImag, Float32Array yReal, Float32Array yImag)
This method performs division of a complex scalar and a complex array.

For each k in the range 0 ≤ k < min(dstReal.length, dstImag.length, yReal.length, yImag.length), the method MUST compute: 

denom = yReal[k] * yReal[k] + yImag[k] * yImag[k]
dstReal[k] = (xReal * yReal[k] + xImag * yImag[k]) / denom
dstImag[k] = (xImag * yReal[k] - xReal * yImag[k]) / denom
void madd (Float32Array dst, Float32Array x, Float32Array y, Float32Array z)
For each k in the range 0 ≤ k < min(dst.length, x.length, y.length, z.length), the method MUST compute dst[k] = x[k] * y[k] + z[k].
void madd (Float32Array dst, double x, Float32Array y, Float32Array z)
For each k in the range 0 ≤ k < min(dst.length, y.length, z.length), the method MUST compute dst[k] = x * y[k] + z[k].
void abs (Float32Array dst, Float32Array x)
For each k in the range 0 ≤ k < min(dst.length, x.length), the method MUST compute dst[k] = abs(x[k]).
void absCplx (Float32Array dst, Float32Array xReal, Float32Array xImag)
This method calculates the absolute value (i.e. magnitude) of a complex array.

For each k in the range 0 ≤ k < min(dst.length, xReal.length, xImag.length), the method MUST compute dst[k] = sqrt(xReal[k]2 + xImag[k]2).

void acos (Float32Array dst, Float32Array x)
For each k in the range 0 ≤ k < min(dst.length, x.length), the method MUST compute dst[k] = acos(x[k]).
void asin (Float32Array dst, Float32Array x)
For each k in the range 0 ≤ k < min(dst.length, x.length), the method MUST compute dst[k] = asin(x[k]).
void atan (Float32Array dst, Float32Array x)
For each k in the range 0 ≤ k < min(dst.length, x.length), the method MUST compute dst[k] = atan(x[k]).
void atan2 (Float32Array dst, Float32Array y, Float32Array x)
For each k in the range 0 ≤ k < min(dst.length, x.length, y.length), the method MUST compute dst[k] = atan2(y[k], x[k]).
void ceil (Float32Array dst, Float32Array x)
For each k in the range 0 ≤ k < min(dst.length, x.length), the method MUST compute dst[k] = ceil(x[k]).
void cos (Float32Array dst, Float32Array x)
For each k in the range 0 ≤ k < min(dst.length, x.length), the method MUST compute dst[k] = cos(x[k]).
void exp (Float32Array dst, Float32Array x)
For each k in the range 0 ≤ k < min(dst.length, x.length), the method MUST compute dst[k] = exp(x[k]).
void floor (Float32Array dst, Float32Array x)
For each k in the range 0 ≤ k < min(dst.length, x.length), the method MUST compute dst[k] = floor(x[k]).
void log (Float32Array dst, Float32Array x)
For each k in the range 0 ≤ k < min(dst.length, x.length), the method MUST compute dst[k] = log(x[k]).
double max (Float32Array x)
The method MUST return the largest value among all the elements of the array x. If the array is empty, it MUST return −∞.
double min (Float32Array x)
The method MUST return the smallest value among all the elements of the array x. If the array is empty, it MUST return ∞.
void pow (Float32Array dst, Float32Array x, Float32Array y)
For each k in the range 0 ≤ k < min(dst.length, x.length, y.length), the method MUST compute dst[k] = pow(x[k], y[k]).
void pow (Float32Array dst, Float32Array x, double y)
For each k in the range 0 ≤ k < min(dst.length, x.length), the method MUST compute dst[k] = pow(x[k], y).
void random (Float32Array dst, optional double low=0, optional double high=1)

For each k in the range 0 ≤ k < dst.length, the method MUST compute dst[k] = random() * (high - low) + low.

void round (Float32Array dst, Float32Array x)
For each k in the range 0 ≤ k < min(dst.length, x.length), the method MUST compute dst[k] = round(x[k]).
void sin (Float32Array dst, Float32Array x)
For each k in the range 0 ≤ k < min(dst.length, x.length), the method MUST compute dst[k] = sin(x[k]).
void sqrt (Float32Array dst, Float32Array x)
For each k in the range 0 ≤ k < min(dst.length, x.length), the method MUST compute dst[k] = sqrt(x[k]).
void tan (Float32Array dst, Float32Array x)
For each k in the range 0 ≤ k < min(dst.length, x.length), the method MUST compute dst[k] = tan(x[k]).
void clamp (Float32Array dst, Float32Array x, double xMin, double xMax)
For each k in the range 0 ≤ k < min(dst.length, x.length), the method MUST compute dst[k] = clamp(x[k], xMin, xMax).
void fract (Float32Array dst, Float32Array x)
For each k in the range 0 ≤ k < min(dst.length, x.length), the method MUST compute dst[k] = modulo(x[k], 1).
void fill (Float32Array dst, double value)
For each k in the range 0 ≤ k < dst.length, the method MUST perform the assignment dst[k] = value.
void ramp (Float32Array dst, double first, double last)

If dst contains at least one element, the method MUST set the first element to first.

For each k in the range 1 ≤ k < dst.length, the method MUST compute dst[k] = first + k * ((last - first) / (dst.length - 1)).

void sign (Float32Array dst, Float32Array x)
For each k in the range 0 ≤ k < min(dst.length, x.length), the method MUST compute dst[k] = sign(x[k]).
double sum (Float32Array x)

The method MUST return the sum of all the elements of the array x. If the array is empty, it MUST return 0 (zero).

It is RECOMMENDED that the method is implemented using a summation algorithm with an accumulated round-off error growth equal to or better than that of pairwise summation.

void sampleLinear (Float32Array dst, Float32Array x, Float32Array t)

This method implements linear interpolation of the source array x. The array is sampled at the points given in the array t, and stored in the destination array dst.

Values in t are clamped to the valid index range of the array x, i.e. [0, x.length).

The method MUST perform one of the following two operations:

  1. If the array x is empty (has zero elements), the method MUST set each element dst[k] to zero, where k is in the range 0 ≤ k < min(dst.length, t.length).
  2. If the array x has at least one element, the method MUST perform the following operation for each k in the range 0 ≤ k < min(dst.length, t.length): 
    t' = clamp(t[k], 0, x.length - 1)
idx = floor(t')
w = t' - idx
p1 = x[idx]
p2 = x[min(idx + 1, x.length - 1)]
dst[k] = (1 - w) * p1 + w * p2
void sampleLinearRepeat (Float32Array dst, Float32Array x, Float32Array t)

This method implements linear interpolation of the source array x. The array is sampled at the points given in the array t, and stored in the destination array dst.

Values in t are taken modulo x.length.

The method MUST perform one of the following two operations:

  1. If the array x is empty (has zero elements), the method MUST set each element dst[k] to zero, where k is in the range 0 ≤ k < min(dst.length, t.length).
  2. If the array x has at least one element, the method MUST perform the following operation for each k in the range 0 ≤ k < min(dst.length, t.length): 
    t' = modulo(t[k], x.length)
idx = floor(t')
w = t' - idx
p1 = x[idx]
p2 = x[modulo(idx + 1, x.length)]
dst[k] = (1 - w) * p1 + w * p2
void sampleCubic (Float32Array dst, Float32Array x, Float32Array t)

This method implements cubic Catmull-Rom spline interpolation of the source array x. The array is sampled at the points given in the array t, and stored in the destination array dst.

Values in t are clamped to the valid index range of the array x, i.e. [0, x.length).

The method MUST perform one of the following two operations:

  1. If the array x is empty (has zero elements), the method MUST set each element dst[k] to zero, where k is in the range 0 ≤ k < min(dst.length, t.length).
  2. If the array x has at least one element, the method MUST perform the following operation for each k in the range 0 ≤ k < min(dst.length, t.length): 
    t' = clamp(t[k], 0, x.length - 1)
idx = floor(t')
w = t' - idx
w2 = w*w
w3 = w*w*w
h1 =  2*w3 - 3*w2 + 1
h2 = -2*w3 + 3*w2
h3 = 0.5 * (w3 - 2*w2 + w)
h4 = 0.5 * (w3 -   w2)
p1 = x[max(idx - 1, 0)]
p2 = x[idx]
p3 = x[min(idx + 1, x.length - 1)]
p4 = x[min(idx + 2, x.length - 1)]
dst[k] = h1 * p2 + h2 * p3 + h3 * (p3 - p1) + h4 * (p4 - p2)
void sampleCubicRepeat (Float32Array dst, Float32Array x, Float32Array t)

This method implements cubic Catmull-Rom spline interpolation of the source array x. The array is sampled at the points given in the array t, and stored in the destination array dst.

Values in t are taken modulo x.length.

The method MUST perform one of the following two operations:

  1. If the array x is empty (has zero elements), the method MUST set each element dst[k] to zero, where k is in the range 0 ≤ k < min(dst.length, t.length).
  2. If the array x has at least one element, the method MUST perform the following operation for each k in the range 0 ≤ k < min(dst.length, t.length): 
    t' = modulo(t[k], x.length)
idx = floor(t')
w = t' - idx
w2 = w*w
w3 = w*w*w
h1 =  2*w3 - 3*w2 + 1
h2 = -2*w3 + 3*w2
h3 = 0.5 * (w3 - 2*w2 + w)
h4 = 0.5 * (w3 -   w2)
p1 = x[modulo(idx - 1, x.length)]
p2 = x[idx]
p3 = x[modulo(idx + 1, x.length)]
p4 = x[modulo(idx + 2, x.length)]
dst[k] = h1 * p2 + h2 * p3 + h3 * (p3 - p1) + h4 * (p4 - p2)
void pack (TypedArray dst, unsigned long offset, unsigned long stride, Float32Array src1, optional Float32Array src2, optional Float32Array src3, optional Float32Array src4)
This method packs (interleaves) the unpacked data in src1, src2, src3 and src4, and stores the data interleaved in the dst array, which may be of any Typed Array type.

This method MUST perform the following steps:

  1. If any of the arguments src2, src3 or src4, provided that it was passed as an argument, has a length that differs from src1.length, the method MUST throw a "NotSupportedError" exception.
  2. If stride is less than 1, the method MUST throw a "NotSupportedError" exception.
  3. Let dstCount be floor((dst.length - offset) / stride).
  4. For each k in the range 0 ≤ k < min(dstCount, src1.length), the method MUST perform the following operation, where type conversions follow the rules of [[!WEBIDL]]: 
    dst[offset + stride*k] = src1[k];
if (src2 argument is defined)
  dst[offset + stride*k + 1] = src2[k];
if (src3 argument is defined)
  dst[offset + stride*k + 2] = src3[k];
if (src4 argument is defined)
  dst[offset + stride*k + 3] = src4[k];
void unpack (TypedArray src, unsigned long offset, unsigned long stride, Float32Array dst1, optional Float32Array dst2, optional Float32Array dst3, optional Float32Array dst4)
This method unpacks (de-interleaves) the packed data in src, which may be of any Typed Array type, and stores the de-interleaved data in the arrays dst1, dst2, dst3 and dst4.

This method MUST perform the following steps:

  1. If any of the arguments dst2, dst3 or dst4, provided that it was passed as an argument, has a length that differs from dst1.length, the method MUST throw a "NotSupportedError" exception.
  2. If stride is less than 1, the method MUST throw a "NotSupportedError" exception.
  3. Let srcCount be floor((src.length - offset) / stride).
  4. For each k in the range 0 ≤ k < min(srcCount, dst1.length), the method MUST perform the following operation, where type conversions follow the rules of [[!WEBIDL]]: 
    dst1[k] = src[offset + stride*k];
if (dst2 argument is defined)
  dst2[k] = src[offset + stride*k + 1];
if (dst3 argument is defined)
  dst3[k] = src[offset + stride*k + 2];
if (dst4 argument is defined)
  dst4[k] = src[offset + stride*k + 3];

Filter interface

A Filter object implements FIR (finite impulse response) and IIR (infinite impulse response) filters of any order. The type and order of the filter is defined by setting the a and b filter coefficients.

If a is an empty array (if it has zero elements), the object will effectively implement a convolution between the signal x and the impulse response b.

A Filter object maintains intra-call state in order to facilitate continuous filter operation across several buffers.

We should add versions of setA() and setB() that support regular ECMAScript arrays. Same thing with the constructor (to be able to pass arrays directly instead of using setA() / setB()).

A Filter object MUST internally store filter coefficients — hereafter referred to as the arrays a (for the recursive coefficients) and b (for the convolution coefficients).

A Filter object MUST maintain a history state — hereafter referred to as the history — of the input and output signals. This state must keep enough data from previously processed input and output signals to make the filtering function well defined.

At a minimum, the history must contain the b.length - 1 last processed input samples and the a.length last processed output samples from previous calls to the filter() method.

The Filter() constructor, when invoked, MUST return a newly created Filter object. Furthermore, the constructor MUST perform the following steps before returning the new Filter object:

  1. If bSize is less than 1, a "NotSupportedError" exception MUST be thrown.
  2. Create the internal array b of length bSize, where the first element is assigned the value 1, and the rest of the elements are set to 0 (zero).
  3. Create the internal array a of length aSize, and set all the elements to 0 (zero).
  4. Clear the history, as if clearHistory() was called.

If there was not enough memory to create the Filter object, a "NotSupportedError" exception MUST be thrown.

void filter (Float32Array dst, Float32Array x)
Applies the configured filter to the array x and stores the result in the array dst.

The method MUST perform the following steps:

  1. If the length of x differs from the length of dst, a "NotSupportedError" exception MUST be thrown.
  2. Perform the filtering operation according to the following pseudo code: 
    for k = 0 to x.length - 1 do
  res = 0;
  for m = 0 to b.length - 1 do
    res = res + b[m] * x[k - m];
  end
  for m = 0 to a.length - 1 do
    res = res - a[m] * dst[k - 1 - m];
  end
  dst[k] = res;
end
    For negative indexes into the arrays x and dst, the internal history state from previous calls to this method MUST be used.
  3. Update the internal history state so that it includes the last min(b.length-1, x.length) samples from the array x and the last min(a.length, dst.length) samples from the array dst.
void clearHistory ()
This method MUST clear the internal history state, so that the next filter() call acts as if at least b.length - 1 input samples with the value 0 (zero) have previously been processed, and at least a.length output samples with the value 0 (zero) have previously been generated.
void setB (Float32Array values)
This method MUST set the first min(b.length, values.length) coefficients of the internal b array of the filter object.
void setA (Float32Array values)
This method MUST set the first min(a.length, values.length) coefficients of the internal a array of the filter object.

FFT interface

An FFT object implements forward and inverse discrete Fourier transforms.
The element layout of the frequency domain (forward transform) has to be defined.

The FFT() constructor, when invoked, MUST return a newly created FFT object.

If the constructor was passed an argument, then the FFT object's size property MUST be set to the value of that argument. Otherwise, the object's size property MUST be set to 256.

If a value of less than 1 is passed as the argument to the FFT() constructor, a "NotSupportedError" exception MUST be thrown.

If there was not enough memory to create the FFT object, a "NotSupportedError" exception MUST be thrown.

The size property tells the number of frequency bins to use in subsequent calls to the transformation methods of the FFT object: forward, forwardCplx, inverse and inverseCplx.

void forward (Float32Array dstReal, Float32Array dstImag, Float32Array x)

This method MUST perform the same operation as if forwardCplx was called using the following syntax: forwardCplx (dstReal, dstImag, x, dummy), where dummy is a zero-filled Float32Array of length size.

void forwardCplx (Float32Array dstReal, Float32Array dstImag, Float32Array xReal, Float32Array xImag)

This method MUST compute the forward Fourier transform of the first size elements of the complex input signal (xReal and xImag), and store the result in the first size elements of the complex output signal (dstReal and dstImag).

The generated output signal MUST be scaled such that the RMS of it equals the RMS of the input signal.

If any of the arrays xReal, xImag, dstReal and dstImag has less than size elements, the method MUST throw a "NotSupportedError" exception.

void inverse (Float32Array dst, Float32Array xReal, Float32Array xImag)

This method MUST perform the same operation as if inverseCplx was called using the following syntax: inverseCplx (dst, dummy, xReal, xImag), where dummy is a Float32Array of length size.

void inverseCplx (Float32Array dstReal, Float32Array dstImag, Float32Array xReal, Float32Array xImag)

This method MUST compute the inverse Fourier transform of the first size elements of the complex input signal (xReal and xImag), and store the result in the first size elements of the complex output signal (dstReal and dstImag).

The generated output signal MUST be scaled such that the RMS of it equals the RMS of the input signal.

If any of the arrays xReal, xImag, dstReal and dstImag has less than size elements, the method MUST throw a "NotSupportedError" exception.

readonly attribute unsigned long size

Reading this attribute MUST return the size of the Fourier transform, i.e. the number of frequency bins that are used during a transform.

Examples

Below are some examples of how to use different interfaces and methods provided by this specification. Unless otherwise noted, they are written in ECMAScript (with [[!TYPED-ARRAYS]] support).

Some of the examples create temporary objects and arrays. In real world applications, it is expected that such objects are kept as members of long-lived objects in order to minimize heap allocation activity and garbage collection.

Generate a sine signal with 20 cycles: 

var y = new Float32Array(1000);
ArrayMath.ramp(y, 0, 2*Math.PI*20);
ArrayMath.sin(y, y);

Generate a low-pass filtered noise signal, using a first order IIR filter: 

var x = new Float32Array(1000);
ArrayMath.random(x, -1, 1);

var y = new Float32Array(x.length);
var filt = new Filter(1, 1);
filt.setB(new Float32Array([0.01]));
filt.setA(new Float32Array([-0.99]));
filt.filter(y, x);  // y[k] = 0.01 * x[k] + 0.99 * y[k-1]

Calculate the RMS (root mean square) of a signal: 

function rms (x) {
  var y = new Float32Array(x.length);
  ArrayMath.mul(y, x, x);
  return Math.sqrt(ArrayMath.sum(y) / y.length);
}

Normalize a signal to a given peak amplitude: 

function normalize (x, magnitude) {
  var peak = Math.max(Math.abs(ArrayMath.max(x)), Math.abs(ArrayMath.min(x)));
  if (peak > 0)
    ArrayMath.mul(x, magnitude / peak, x);
}

Create a Hamming window: 

function hamming (size) {
  var w = new Float32Array(size);
  ArrayMath.ramp(w, 0, 2*Math.PI);
  ArrayMath.cos(w, w);
  ArrayMath.mul(w, -0.46, w);
  ArrayMath.add(w, 0.54, w);
  return w;
}