barycentric_interpolator.ocl

#define T double
kernel void barycentric_interpolator(
	global const T* xp,
	global T* yp,
	int count,
	global const T* xs,
	global const T* ws,
	global const T* ys
) {
	int gid = get_global_id(0);
	const T x = xp[gid];
	T u = 0;
	T v = 0;

	for(int i = 0; i < count; i++) {
		const T xi = xs[i];
		const T wi = ws[i];
		const T yi = ys[i];

		const T ai = wi / (xi != x ? xi - x : 1);
		u += yi * ai;
		v += ai;
	}

	yp[gid] = u / v;
}

stern_brocot_tone_barycentric.py

import sys
import argparse
import numpy as np
import copy
import pyopencl as cl
import pyopencl.array
import pyopencl.scan
import scipy.io.wavfile
from scipy.interpolate import BarycentricInterpolator

def main():
    parser = argparse.ArgumentParser(description='Generate Stern Brocot tone.')
    parser.add_argument('outname', help='output file')
    parser.add_argument('-t', '--level-time', default=2, help='time per level in seconds', type=float)
    parser.add_argument('-a', '--amplitude', default=0.8, help='maximum amplitude', type=float)
    parser.add_argument('-r', '--rate', default=44100, help='sample rate', type=int)
    parser.add_argument('-d', '--depth', default=6, help='number of levels', type=int)
    parser.add_argument('-p0', default='1/1', help='bottom ratio',type=str)
    parser.add_argument('-p1', default='1/0', help='top ratio',type=str)
    parser.add_argument('--octave', default=4, help='base octave', type=int)
    parser.add_argument('--note', default=0, help='base note', type=int)
    parser.add_argument('--logspace', action='store_true', help='use log space branches')
    parser.add_argument('--verbose', '-v', action='count', default=0)
    args = parser.parse_args()

    depth = args.depth
    rate = args.rate
    total_time = args.level_time * (args.depth+1)

    pq0 = tuple(map(int,args.p0.split('/')))
    pq1 = tuple(map(int,args.p1.split('/')))
    dtype_name = 'double'
    float_dtype = np.float64
    note = args.octave + args.note / 12
    fund_freq = 440 * pow(2, -1 + 3 / 12 - 4 + note)
    total_samples = round(total_time * rate)
    wav = np.zeros(total_samples, dtype=float_dtype)

    def rec(pq0, pq1, depth, interpolator=None, path=[]):
        level = len(path)
        p0, q0 = pq0
        p1, q1 = pq1
        pq2 = p0 + p1, q0 + q1
        p2, q2 = pq2

        x2 = float(p2) / float(q2)
        if interpolator is None:
            interpolator2 = BarycentricInterpolator([level], [x2])
        else:
            interpolator2 = copy.deepcopy(interpolator)
            interpolator2.add_xi([level], [x2])
        if level - 1 >= depth:
            interpolator2.add_xi([level + 1], [x2])
            yield interpolator2
            return

        yield from rec(pq0, pq2, depth, interpolator2, path + [0])
        yield from rec(pq2, pq1, depth, interpolator2, path + [1])

    mf = cl.mem_flags
    ctx = cl.create_some_context()

    with open('barycentric_interpolator.ocl', 'rt') as infile:
        program_source = infile.read()
    program = cl.Program(ctx, program_source).build()

    scanner = pyopencl.scan.InclusiveScanKernel(
        ctx=ctx,
        dtype=float_dtype,
        scan_expr='a+b',
        neutral='0',
    )
    if args.logspace:
        exp_elementwise = pyopencl.elementwise.ElementwiseKernel(
            context=ctx,
            arguments=f'__global {dtype_name} * ys',
            operation='ys[i] = exp(ys[i]);',
        )

    elementwise = pyopencl.elementwise.ElementwiseKernel(
        context=ctx,
        arguments=f'{dtype_name} c0, {dtype_name} c1, __global {dtype_name} * ys, __global {dtype_name} * wav',
        operation='wav[i] += sin(c0 * ys[i] + c1);',
    )

    with cl.CommandQueue(ctx) as queue:
        wav_ocl_arr = cl.array.Array(queue, (total_samples,), dtype=float_dtype)
        cl.enqueue_fill_buffer(queue, wav_ocl_arr.data, offset=0, pattern=float_dtype(0), size=wav.nbytes)

        xs = np.linspace(0, depth + 1, len(wav), dtype=float_dtype)
        xs_ocl = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=xs)
        ys_ocl_arr = cl.array.Array(queue, (total_samples,), dtype=float_dtype)
        for i, interpolator in enumerate(rec(pq0, pq1, depth, BarycentricInterpolator([-1], [1]))):
            phase = np.random.random() * 2 * np.pi
            xi = interpolator.xi.flatten().astype(float_dtype)
            wi = interpolator.wi.flatten().astype(float_dtype)
            yi = interpolator.yi.flatten().astype(float_dtype)
            if args.logspace:
                yi = np.log(yi)

            xi_ocl = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=xi)
            wi_ocl = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=wi)
            yi_ocl = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=yi)
            program.barycentric_interpolator(
                queue,
                (len(xs),),
                None,
                xs_ocl,
                ys_ocl_arr.data,
                np.int32(len(interpolator.xi)),
                xi_ocl,
                wi_ocl,
                yi_ocl,
            ).wait()
            if args.logspace:
                exp_elementwise(ys_ocl_arr, queue=queue)

            scanner(ys_ocl_arr, queue=queue)
            knl = elementwise(float_dtype(fund_freq * 2 * np.pi / rate), float_dtype(phase), ys_ocl_arr, wav_ocl_arr, queue=queue)
            knl.wait()
            if args.verbose >= 1:
                print(f'{i/(2**(depth+1)-1):2.2%} complete')
        cl.enqueue_copy(queue, wav, wav_ocl_arr.data)
    wav *= args.amplitude / np.max(np.abs(wav))
    scipy.io.wavfile.write(args.outname, rate, wav)


if __name__ == '__main__':
    main()