importing some org files

Paul Batchelor committed 6 years ago

Commit 5f348381 parent 8888fbf12 files · +3764 −1Browse files

Unified · Split

12 files

README.md +2 −1
bigverb.org +1035 −0
dcblocker.org +118 −0
fmpair.org +576 −0
index.org +21 −0
modalres.org +306 −0
osc.org +409 −0
peakeq.org +358 −0
phasewarp.org +104 −0
rline.org +403 −0
scrubber.org +56 −0
vardelay.org +376 −0

@@ -1 +1,2 @@
1	# sndkit	1	# sndkit
		2	sndkit is a sonic toolkit for everywhere.

bigverb.org+1035−0

@@ -0,0 +1,1035 @@
	1	#+TITLE: Bigverb
	2	* Introduction
	3	Bigverb is an implementation of a classic waveguide
	4	reverberator. It is composed up
	5	of 8 feedback delay lines running in parallel, in what is
	6	known as a [[https://ccrma.stanford.edu/~jos/cfdn/Feedback_Delay_Networks.html][Feedback Delay Network]]. It is largely inspired by
	7	the works Julius O. Smith, and the work he did in using
	8	waveguide networks. In particular, this algorithm is based
	9	on ideas in a 1985 paper called "A New Approach to Digital
	10	Reverberation using Waveguide Networks".
	11
	12	Bigverb, as the name implies, has a very nice warm sound
	13	when set to be a very large space, and is ideal for use in
	14	ambient style music. Bigverb is not designed to sound
	15	realistic, and is not an ideal choice for anything that
	16	needs room modelling.
	17	* Distinct Characterstics
	18	The "sound" of this algorithm comes from a few major
	19	things.
	20	** Use of Jitter Modulation
	21	Every delay line has a bit of jitter applied to
	22	the delay time, which causes a little of pitch modulation.
	23	Several schroeder and FDN-style reverb designs use this
	24	kind of modulation, such as the Mutable instruments
	25	clouds reverb, and the tank reverb described in the Dattorro
	26	paper "Effect design, part 1: reverberator and other filters".
	27	Such modulation creates the illusion of the reverb being
	28	a higher-order system than it actual (aka: a more complex
	29	reverb algorithm with more delay lines).
	30	While these designs tend to use a sinusoidal LFO on only a
	31	few delay lines, the Costello topology applies non-correlated
	32	modulation to every single delay line. As a result, you get
	33	a very dense-sounding reverb output.
	34	** Parameter Tunings
	35	The parameter tunings of this reverb add to the
	36	particular character of this sound. All units are set in
	37	such a way so that they can be stored as integer values.
	38
	39	The tuning consists of 8 parameter sets, corresponding to
	40	each delay line "unit". A set has the following parameters: delay time
	41	(=delay=) (in units of samples with a samplerate of
	42	44.1kHz), the amount of variation in delay time to apply
	43	(=drift=) (in units of tenth milliseconds, in order to use
	44	integer values), the random variation frequency (=randfreq=)
	45	(in Hz multiplied by 1000 to avoid decimals),
	46	as well as a starting seed for the internal RNG.
	47
	48	The credit for these parameter tunings (and this topology)
	49	go to Sean Costello, who also originally designed this
	50	reverb algorithm in 1999. This work would later pave the way
	51	for his famous commercial reverb plugins and algorithms
	52	published by his company ValhallaDSP.
	53
	54	#+NAME: parameter_set
	55	#+BEGIN_SRC c
	56	struct bigverb_paramset {
	57	int delay; /* in samples, 44.1 kHz */
	58	int drift; /* 1/10 milliseconds */
	59	int randfreq; /* Hertz * 1000 */
	60	int seed;
	61	};
	62
	63	static const struct bigverb_paramset params[8] = {
	64	{0x09a9, 0x0a, 0xc1c, 0x07ae},
	65	{0x0acf, 0x0b, 0xdac, 0x7333},
	66	{0x0c91, 0x11, 0x456, 0x5999},
	67	{0x0de5, 0x06, 0xf85, 0x2666},
	68	{0x0f43, 0x0a, 0x925, 0x50a3},
	69	{0x101f, 0x0b, 0x769, 0x5999},
	70	{0x085f, 0x11, 0x37b, 0x7333},
	71	{0x078d, 0x06, 0xc95, 0x3851}
	72	};
	73	#+END_SRC
	74	** Delay Line Configuration
	75	Normally, this family of reverbs is a combination of allpass
	76	filters and comb filters in various series/parallel
	77	combinations. This implementation has a far simpler
	78	topology: it is essentially (but not quite) 8 feedback delay
	79	lines in parallel with a 1-pole lowpass IIR filter for tone
	80	control. If it weren't for the massive amount delay
	81	modulation, this would be a pretty metallic sounding reverb!
	82	* Top Level Files and Structs
	83	** Files to Tangle To
	84	This Bigverb document tanlges to two files:
	85	=bigverb.c= and =bigverb.h=.
	86
	87	=bigverb.h= contains typedefs, structs, and forward
	88	function declarations. To enable structs, define the
	89	macro =SK_BIGVERB_PRIV=.
	90
	91	#+NAME: bigverb.c
	92	#+BEGIN_SRC c :tangle bigverb.c
	93	#include <stdlib.h>
	94	#include <stdio.h>
	95	#include <math.h>
	96	#define SK_BIGVERB_PRIV
	97	#include "bigverb.h"
	98	<<parameter_set>>
	99	<<fdelay_constants>>
	100	<<static_funcdefs>>
	101	<<funcs>>
	102	#+END_SRC
	103
	104	#+NAME: bigverb.h
	105	#+BEGIN_SRC c :tangle bigverb.h
	106	#ifndef SK_BIGVERB_H
	107	#define SK_BIGVERB_H
	108
	109	#ifndef SKFLT
	110	#define SKFLT float
	111	#endif
	112	<<typedefs>>
	113	<<funcdefs>>
	114
	115	#ifdef SK_BIGVERB_PRIV
	116	<<main_struct>>
	117	#endif
	118	#endif
	119	#+END_SRC
	120	** The Bigverb Struct
	121	An instance of Bigverb is contained inside of a struct
	122	called =sk_bigverb=.
	123	#+NAME: typedefs
	124	#+BEGIN_SRC c
	125	typedef struct sk_bigverb sk_bigverb;
	126	#+END_SRC
	127
	128	#+NAME: main_struct
	129	#+BEGIN_SRC c
	130	<<delay_struct>>
	131	struct sk_bigverb {
	132	int sr;
	133	<<sk_bigverb>>
	134	};
	135	#+END_SRC
	136	* Setup and Cleanup
	137	A new instance of bigverb is created with
	138	=sk_bigverb_new=. The only argument required is the sampling
	139	rate. If something goes wrong, this will return =NULL=.
	140
	141	#+NAME: funcdefs
	142	#+BEGIN_SRC c
	143	sk_bigverb * sk_bigverb_new(int sr);
	144	#+END_SRC
	145
	146	#+NAME: funcs
	147	#+BEGIN_SRC c
	148	sk_bigverb * sk_bigverb_new(int sr)
	149	{
	150	sk_bigverb *bv;
	151
	152	bv = calloc(1, sizeof(sk_bigverb));
	153
	154	bv->sr = sr;
	155	<<init_variables>>
	156	<<setup_delay_lines>>
	157
	158	return bv;
	159	}
	160	#+END_SRC
	161
	162	When it is done being used, bigverb must be cleanly freed
	163	with =sk_bigverb_del=.
	164
	165	#+NAME: funcdefs
	166	#+BEGIN_SRC c
	167	void sk_bigverb_del(sk_bigverb *bv);
	168	#+END_SRC
	169
	170	#+NAME: funcs
	171	#+BEGIN_SRC c
	172	void sk_bigverb_del(sk_bigverb *bv)
	173	{
	174	<<cleanup>>
	175	free(bv);
	176	bv = NULL;
	177	}
	178	#+END_SRC
	179	* High level parameters
	180	High level parametric control of bigverb includes
	181	"size" and "cutoff". Set parameters before computing
	182	audio. These are are just floating point values contained
	183	in the struct that can be indirectly set with setters in
	184	situations where the struct is opaque.
	185	** Size
	186	Set the reverb size with =sk_bigverb_size=
	187
	188	#+NAME: funcdefs
	189	#+BEGIN_SRC c
	190	void sk_bigverb_size(sk_bigverb *bv, SKFLT size);
	191	#+END_SRC
	192
	193	Size is a variable between 0-1, which
	194	controls the feedback level for the delay line.
	195
	196	The =size= parameter is stored as a variable called
	197	=size=, and is set to be a pretty sounding value of
	198	=0.93=.
	199
	200	#+NAME: sk_bigverb
	201	#+BEGIN_SRC c
	202	SKFLT size;
	203	#+END_SRC
	204
	205	#+NAME: init_variables
	206	#+BEGIN_SRC c
	207	sk_bigverb_size(bv, 0.93);
	208	#+END_SRC
	209
	210	#+NAME: funcs
	211	#+BEGIN_SRC c
	212	void sk_bigverb_size(sk_bigverb *bv, SKFLT size)
	213	{
	214	bv->size = size;
	215	}
	216	#+END_SRC
	217	** Cutoff
	218	The tone of bigverb can be set with =sk_bigverb_cutoff=.
	219
	220	#+NAME: funcdefs
	221	#+BEGIN_SRC c
	222	void sk_bigverb_cutoff(sk_bigverb *bv, SKFLT cutoff);
	223	#+END_SRC
	224
	225	=cutoff= is a parameter in Hz that determines the overall
	226	timbre of the reverb. This controls the cutoff frequency of
	227	the one pole lowpass filter applied to the reverb.
	228
	229	It is set to be a default value of 10kHz, or 10,000 hz.
	230
	231	#+NAME: init_variables
	232	#+BEGIN_SRC c
	233	sk_bigverb_cutoff(bv, 10000.0);
	234	#+END_SRC
	235
	236	Cutoff uses caching in order to monitor if the parameter
	237	has changed. It does this in order to prevent needing to
	238	compute filter coefficients every sample. The main variable
	239	to be set is =cutoff=, and the cached variable is =pcutoff=.
	240	At the beginning, =pcutoff= is set to be a negative value,
	241	which will cause bigverb to calculate coefficients in the
	242	first call to the tick function after initialization.
	243
	244	#+NAME: sk_bigverb
	245	#+BEGIN_SRC c
	246	SKFLT cutoff;
	247	SKFLT pcutoff;
	248	#+END_SRC
	249
	250	#+NAME: init_variables
	251	#+BEGIN_SRC c
	252	bv->pcutoff = -1;
	253	#+END_SRC
	254
	255	#+NAME: funcs
	256	#+BEGIN_SRC c
	257	void sk_bigverb_cutoff(sk_bigverb *bv, SKFLT cutoff)
	258	{
	259	bv->cutoff = cutoff;
	260	}
	261	#+END_SRC
	262	* Filter
	263	State in a constant called =filt=.
	264
	265	#+NAME: sk_bigverb
	266	#+BEGIN_SRC c
	267	SKFLT filt;
	268	#+END_SRC
	269
	270	#+NAME: init_variables
	271	#+BEGIN_SRC c
	272	bv->filt = 1.0;
	273	#+END_SRC
	274	* Computing Audio
	275	After bigverb has been initialized, it is ready to
	276	process audio. This implementation uses what is known
	277	as a =tick= function, or a function that computes audio
	278	one sample at a time instead of one block at a time. This
	279	simplifies the implementation at the cost of a little bit
	280	of performance overhead, depending on the compiler and
	281	optimization settings.
	282	** Top-Level Tick Function
	283	The function to tick one sample unit of audio is done with
	284	=sk_bigverb_tick=. It takes in two stereo input
	285	values, and returns two stereo output values.
	286
	287	#+NAME: funcdefs
	288	#+BEGIN_SRC c
	289	void sk_bigverb_tick(sk_bigverb *bv,
	290	SKFLT inL, SKFLT inR,
	291	SKFLT outL, SKFLT outR);
	292	#+END_SRC
	293
	294	#+NAME: funcs
	295	#+BEGIN_SRC c
	296	void sk_bigverb_tick(sk_bigverb *bv,
	297	SKFLT inL, SKFLT inR,
	298	SKFLT outL, SKFLT outR)
	299	{
	300	/* TODO: implement */
	301	SKFLT lsum, rsum;
	302
	303	lsum = 0;
	304	rsum = 0;
	305
	306	<<update_filter_coefficients>>
	307	<<calculate_junction_pressure>>
	308	<<compute_delay_bank>>
	309
	310	*outL = lsum;
	311	*outR = rsum;
	312	}
	313	#+END_SRC
	314	** Updating filter coefficients
	315	Bigverb uses parameter caching for the =cutoff= parameter in
	316	order to save on computation time.
	317
	318	Any time =cutoff= changes, the filter coefficients must be
	319	updated. This happens in the tick function, before any
	320	computation happens.
	321
	322	The filter is a simple 1-pole IIR lowpass filter whose
	323	difference equation been reduced to only require a single
	324	parameter. This in turn then gets used in each filter delay
	325	line.
	326
	327	#+NAME: update_filter_coefficients
	328	#+BEGIN_SRC c
	329	if (bv->pcutoff != bv->cutoff) {
	330	bv->pcutoff = bv->cutoff;
	331	bv->filt = 2.0 - cos(bv->pcutoff * 2 * M_PI / bv->sr);
	332	bv->filt = bv->filt - sqrt(bv->filt * bv->filt - 1.0);
	333	}
	334	#+END_SRC
	335	** Calculating Resultant Junction Pressure Amount
	336	The resultant junction pressure amount is calculated from
	337	the delay bank, and then factored into the input signals.
	338
	339	Sum of all the delay line signals, and scaled by 0.25, or
	340	2/N, where N is the number of delay lines (8).
	341
	342	#+NAME: calculate_junction_pressure
	343	#+BEGIN_SRC c
	344	{
	345	int i;
	346	SKFLT jp;
	347
	348	jp = 0;
	349
	350	for (i = 0; i < 8; i++) {
	351	jp += bv->delay[i].y;
	352	}
	353
	354	jp *= 0.25;
	355
	356	inL = jp + inL;
	357	inR = jp + inR;
	358	}
	359	#+END_SRC
	360	** Computing the delay bank
	361	The delay bank is then computed. Each delay line is computed
	362	and summed with either the left or right input signal, and
	363	then sent to a corresponding left or right channel.
	364
	365	At the end, a final scaling out of the output happens. This
	366	is hard coded to be 35 percent.
	367
	368	#+NAME: compute_delay_bank
	369	#+BEGIN_SRC c
	370	{
	371	int i;
	372	for (i = 0; i < 8; i++) {
	373	if (i & 1) {
	374	rsum += delay_compute(&bv->delay[i],
	375	inR,
	376	bv->size,
	377	bv->filt,
	378	bv->sr);
	379	} else {
	380	lsum += delay_compute(&bv->delay[i],
	381	inL,
	382	bv->size,
	383	bv->filt,
	384	bv->sr);
	385	}
	386	}
	387	}
	388	rsum *= 0.35f;
	389	lsum *= 0.35f;
	390	#+END_SRC
	391	* The Feedback Delay Line Bank
	392	8 delay units come together to make the delay line bank.
	393	Each is initialized using one of the parameter sets.
	394	** Memory Allocation + Setup
	395	#+NAME: sk_bigverb
	396	#+BEGIN_SRC c
	397	SKFLT *buf;
	398	#+END_SRC
	399
	400	#+NAME: init_variables
	401	#+BEGIN_SRC c
	402	bv->buf = NULL;
	403	#+END_SRC
	404
	405	#+NAME: setup_delay_lines
	406	#+BEGIN_SRC c
	407	{
	408	unsigned long total_size;
	409	int i;
	410	SKFLT *buf;
	411
	412	total_size = 0;
	413	buf = NULL;
	414	<<calculate_pool_size>>
	415	<<allocate_memory>>
	416	<<initialize_delay_banks>>
	417	}
	418	#+END_SRC
	419	The delay bank is the abstraction in charge of properly
	420	allocating all the memory needed for the buffers.
	421
	422	Memory is allocated in one giant chunk, and then divied up
	423	to each delay line.
	424
	425	The total memory size is obtained by summing all the delay
	426	times. These times are stored as fixed delay times in units
	427	of samples. These parameters assume a sampling rate of
	428	44.1kHz. If this is not the case, this value must be scaled
	429	accordingly, and then truncated to be an integer. This value
	430	is used again to properly slice up the big memory chunk.
	431
	432	#+NAME: static_funcdefs
	433	#+BEGIN_SRC c
	434	static int get_delay_size(const struct bigverb_paramset *p, int sr);
	435	#+END_SRC
	436
	437	#+NAME: funcs
	438	#+BEGIN_SRC c
	439	static int get_delay_size(const struct bigverb_paramset *p, int sr)
	440	{
	441	SKFLT sz;
	442	sz = (SKFLT)p->delay/44100 + (p->drift * 0.0001) * 1.125;
	443	return floor(16 + sz*sr);
	444	}
	445	#+END_SRC
	446
	447	#+NAME: calculate_pool_size
	448	#+BEGIN_SRC c
	449	for (i = 0; i < 8; i++) {
	450	total_size += get_delay_size(&params[i], sr);
	451	}
	452	#+END_SRC
	453
	454	Allocation is done with =calloc=, which zeros out the memory
	455	as well. This memory will eventually be freed in
	456	=sk_bigverb_del=.
	457
	458	#+NAME: allocate_memory
	459	#+BEGIN_SRC c
	460	buf = calloc(1, sizeof(SKFLT) * total_size);
	461	bv->buf = buf;
	462	#+END_SRC
	463
	464	#+NAME: cleanup
	465	#+BEGIN_SRC c
	466	free(bv->buf);
	467	#+END_SRC
	468
	469	#+NAME: sk_bigverb
	470	#+BEGIN_SRC c
	471	sk_bigverb_delay delay[8];
	472	#+END_SRC
	473
	474	#+NAME: initialize_delay_banks
	475	#+BEGIN_SRC c
	476	{
	477	unsigned long bufpos;
	478	bufpos = 0;
	479	for (i = 0; i < 8; i++) {
	480	unsigned int sz;
	481	sz = get_delay_size(&params[i], sr);
	482
	483	delay_init(&bv->delay[i], &params[i],
	484	&buf[bufpos], sz, sr);
	485	bufpos += sz;
	486	}
	487	}
	488	#+END_SRC
	489	** A Single Delay Line Unit
	490	A delay unit in a bank consists of variable delay line with
	491	cubic interpolation with a 1 pole low-pass filter for tone
	492	control, whose frequency is determined using a master
	493	parameter, as well as a jitter generator. Feedback as well.
	494	*** Struct Declaration
	495	A delay unit is known as a struct called =sk_bigverb_delay=.
	496
	497	#+NAME: typedefs
	498	#+BEGIN_SRC c
	499	typedef struct sk_bigverb_delay sk_bigverb_delay;
	500	#+END_SRC
	501
	502	#+NAME: delay_struct
	503	#+BEGIN_SRC c
	504	struct sk_bigverb_delay {
	505	<<bigverb_delay>>
	506	};
	507	#+END_SRC
	508	*** Initialization
	509	#+NAME: static_funcdefs
	510	#+BEGIN_SRC c
	511	static void delay_init(sk_bigverb_delay *d,
	512	const struct bigverb_paramset *p,
	513	SKFLT *buf,
	514	size_t sz,
	515	int sr);
	516	#+END_SRC
	517
	518	#+NAME: funcs
	519	#+BEGIN_SRC c
	520	static void delay_init(sk_bigverb_delay *d,
	521	const struct bigverb_paramset *p,
	522	SKFLT *buf,
	523	size_t sz,
	524	int sr)
	525	{
	526	SKFLT readpos;
	527	<<delay_init>>
	528	}
	529	#+END_SRC
	530
	531	Set up buffer + sz
	532
	533	#+NAME: bigverb_delay
	534	#+BEGIN_SRC c
	535	SKFLT *buf;
	536	size_t sz;
	537	#+END_SRC
	538
	539	#+NAME: delay_init
	540	#+BEGIN_SRC c
	541	d->buf = buf;
	542	d->sz = sz;
	543	#+END_SRC
	544
	545	Initialize write position (0), abbreviated as =wpos=.
	546
	547	#+NAME: bigverb_delay
	548	#+BEGIN_SRC c
	549	int wpos;
	550	#+END_SRC
	551
	552	#+NAME: delay_init
	553	#+BEGIN_SRC c
	554	d->wpos = 0;
	555	#+END_SRC
	556
	557	Initialize read position. Based on delay time, =drift= and
	558	initial seed. Read position has to components, an integer
	559	read position, and a floating point read position. These
	560	will be abbreviated =irpos= and =frpos=.
	561
	562	#+NAME: bigverb_delay
	563	#+BEGIN_SRC c
	564	int irpos;
	565	int frpos;
	566	#+END_SRC
	567
	568	Seed value is multiplied by the initial drift value, and then
	569	divided by 32767.
	570
	571	#+NAME: bigverb_delay
	572	#+BEGIN_SRC c
	573	int rng;
	574	#+END_SRC
	575
	576	#+NAME: delay_init
	577	#+BEGIN_SRC c
	578	d->rng = p->seed;
	579	<<setup_readpos>>
	580	#+END_SRC
	581
	582
	583	The initial time is added to this.
	584
	585	=bufsize - (readpos * sr)= <-- this puts the read position at
	586	the end of the buffer.
	587
	588	Truncate (using integer cast).'
	589
	590	#+NAME: setup_readpos
	591	#+BEGIN_SRC c
	592	readpos = ((SKFLT)p->delay / 44100);
	593	readpos += d->rng * (p->drift * 0.0001) / 32768.0;
	594	readpos = sz - (readpos * sr);
	595	d->irpos = floor(readpos);
	596	d->frpos = floor((readpos - d->irpos) * FRACSCALE);
	597	#+END_SRC
	598
	599	Create first random segments.
	600
	601	#+NAME: delay_init
	602	#+BEGIN_SRC c
	603	<<init_jitter>>
	604	generate_next_line(d, sr);
	605	#+END_SRC
	606	*** Top-Level Compute
	607	The delay line computation is done in a tick function. It
	608	takes in an input sample, returns an output sample. In
	609	addition to delay, filtering, feedback, and jittering
	610	happens as well.
	611
	612	Because feedback + filtering are global options, these
	613	are passed in as parameters on the stack. What is required
	614	is the feedback amount, and the calculated filter
	615	coeffecient used in the filter.
	616
	617	#+NAME: static_funcdefs
	618	#+BEGIN_SRC c
	619	static SKFLT delay_compute(sk_bigverb_delay *d,
	620	SKFLT in,
	621	SKFLT fdbk,
	622	SKFLT filt,
	623	int sr);
	624	#+END_SRC
	625
	626	#+NAME: funcs
	627	#+BEGIN_SRC c
	628	static SKFLT delay_compute(sk_bigverb_delay *del,
	629	SKFLT in,
	630	SKFLT fdbk,
	631	SKFLT filt,
	632	int sr)
	633	{
	634	SKFLT out;
	635	SKFLT frac_norm;
	636	SKFLT a, b, c, d;
	637	SKFLT s[4];
	638	out = 0;
	639	<<write_to_buffer>>
	640	<<increment_write_position>>
	641	<<update_fractional_read_position>>
	642	<<update_integer_read_position>>
	643	<<normalize_fractional_component>>
	644	<<calculate_interpolation_coefficients>>
	645	<<read_from_buffer>>
	646	<<compute_interpolation>>
	647	<<increment_fractional_read_position>>
	648	<<apply_feedback_and_filter>>
	649	<<update_jitter>>
	650	return out;
	651	}
	652	#+END_SRC
	653
	654	# since filter is the same cutoff... probably
	655	# only need to filter once to get same result.
	656	# this is the cummutative property I think?
	657	# May do this later.
	658
	659	The following things happen:
	660
	661	Write the to delay buffer and pre-filter the input by
	662	subtracting the filter state =y=.
	663
	664	#+NAME: write_to_buffer
	665	#+BEGIN_SRC c
	666	del->buf[del->wpos] = in - del->y;
	667	#+END_SRC
	668
	669	Increment the write position. If this is greater than the
	670	buffer size, wrap around.
	671
	672	#+NAME: increment_write_position
	673	#+BEGIN_SRC c
	674	del->wpos++;
	675	if (del->wpos >= del->sz) del->wpos -= del->sz;
	676	#+END_SRC
	677
	678	Update the fractional read position. If the read position
	679	exceeds the maximum fractional scale amount, it means it has
	680	bits that must carry over to the integer read position.
	681	After these bits have been carried over, mask out the upper
	682	bits to keep the range in bounds.
	683
	684	#+NAME: update_fractional_read_position
	685	#+BEGIN_SRC c
	686	if (del->frpos >= FRACSCALE) {
	687	del->irpos += del->frpos >> FRACNBITS;
	688	del->frpos &= FRACMASK;
	689	}
	690	#+END_SRC
	691
	692	If needed, update the read position with wrap-around.
	693
	694	#+NAME: update_integer_read_position
	695	#+BEGIN_SRC c
	696	if (del->irpos >= del->sz) del->irpos -= del->sz;
	697	#+END_SRC
	698
	699	Normalize the fractional component so that it is in range
	700	0 and 1. This is done by dividing the amount by fractional
	701	scaling factor =FRACSCALE=.
	702
	703	#+NAME: normalize_fractional_component
	704	#+BEGIN_SRC c
	705	frac_norm = del->frpos / (SKFLT)FRACSCALE;
	706	#+END_SRC
	707
	708	Calculate interpolation coefficients. These are 4
	709	pre-derived coefficents used to compute third-order
	710	lagrangian interpolation.
	711	Derivation of these is currently beyond the scope of this
	712	document. These will be called =a=, =b=, =c=, and =d=,
	713	respectively, and will correspond to =x(n - 1)=, =x(n)=,
	714	=x(n + 1)=, and =x(n + 2)=, respectively.
	715
	716	#+NAME: calculate_interpolation_coefficients
	717	#+BEGIN_SRC c
	718	{
	719	SKFLT tmp[2];
	720	d = ((frac_norm * frac_norm) - 1) / 6.0;
	721	tmp[0] = ((frac_norm + 1.0) * 0.5);
	722	tmp[1] = 3.0 * d;
	723	a = tmp[0] - 1.0 - d;
	724	c = tmp[0] - tmp[1];
	725	b = tmp[1] - frac_norm;
	726	}
	727	#+END_SRC
	728
	729	Read the samples needed, based on the current playhead
	730	position. When the read position is in regular bounds, this
	731	means reading the previous, current, two next samples.
	732	Otherwise, this means the same thing, but with wrap-around
	733	and bounds checks.
	734
	735	#+NAME: read_from_buffer
	736	#+BEGIN_SRC c
	737	{
	738	int n;
	739	SKFLT *x;
	740	n = del->irpos;
	741	x = del->buf;
	742
	743	if (n > 0 && n < (del->sz - 2)) {
	744	s[0] = x[n - 1];
	745	s[1] = x[n];
	746	s[2] = x[n + 1];
	747	s[3] = x[n + 2];
	748	} else {
	749	int k;
	750	n--;
	751	if (n < 0) n += del->sz;
	752	s[0] = x[n];
	753	for (k = 0; k < 3; k++) {
	754	n++;
	755	if (n >= del->sz) n -= del->sz;
	756	s[k + 1] = x[n];
	757	}
	758	}
	759	}
	760	#+END_SRC
	761
	762	Calculate interpolation. Using the coefficents
	763	described above and the fractional component =f=, one can compute
	764	cubic interpolation with the following expression:
	765
	766	$$
	767	y(n) = (a x(n - 1) + b x(n) + c x(n + 1) + d x(n + 2)) \cdot f + x(n)
	768	$$
	769
	770	#+NAME: compute_interpolation
	771	#+BEGIN_SRC c
	772	out = (as[0] + bs[1] + cs[2] + ds[3]) * frac_norm + s[1];
	773	#+END_SRC
	774
	775	Increment fractional read position, as determined by the
	776	jitter.
	777
	778	#+NAME: increment_fractional_read_position
	779	#+BEGIN_SRC c
	780	del->frpos += del->inc;
	781	#+END_SRC
	782
	783	Apply feedback and filter. The feedback will scale the delay
	784	output. The filtering is a difference equation, optimized
	785	and factored to use a minimum number of arithmetic
	786	operations.
	787
	788	#+NAME: apply_feedback_and_filter
	789	#+BEGIN_SRC c
	790	out *= fdbk;
	791	out += (del->y - out) * filt;
	792	del->y = out;
	793	#+END_SRC
	794
	795	# Original difference equation
	796	#
	797	# y(n) = c1 * x(n) + c2 * y(n - 1)
	798	#
	799	# c2 = b - sqrt(b^2 - 1.0)
	800	# c1 = (1 - c2)
	801	# b = 2 - cos(\omega))
	802	# \omega = freq * 2 * PI / sr
	803	#
	804	# can be derived from tone in the following way:
	805	# y(n) = c1 * x(n) + c2 * Y(n - 1)
	806	# = (1 - c2)x(n) + c2 y(n - 1)
	807	# = x(n) - c2x(n) + c2y(n-1)
	808	# = c2y(n - 1) - c2x(n) + x(n)
	809	# = c2*(y(n - 1) - x(n)) + x(n)
	810
	811	Update jitter, if needed. When the counter zeros out
	812	(or worse), it is time to find a new random target to lerp
	813	to.
	814
	815	#+NAME: update_jitter
	816	#+BEGIN_SRC c
	817	del->counter--;
	818	if (del->counter <= 0) {
	819	generate_next_line(del, sr);
	820	}
	821	#+END_SRC
	822	*** Feedback Fractional Delay Line
	823	A delay line is initialized with a pre-allocated zeroed
	824	buffer its size. Memory will be managed outside of this
	825	abstraction.
	826
	827	Being a fractional delay line means the read position has
	828	two components: an integer component and a fractional
	829	component. The integer component is the current position
	830	in the delay buffer. The fractional component tells how much
	831	it goes over into the next discrete sample position.
	832	In a way, interpolation can be thought of as the process of
	833	using these two values to make a really good guess of what
	834	lies in between the samples.
	835
	836	The fractional delay component is maintained as a 28 bit
	837	integer. This is done to avoid some of the weirdness found
	838	in floating point operations. The remaining upper bits
	839	are "carry-over" samples, that get accumulated in integer
	840	component of the read position.
	841
	842	A few constants are used to conveniently work with this
	843	fractional delay component.
	844
	845	=FRACSCALE= is the fractional scaling amount,
	846	which is $1^28$, or =0x10000000=. Multiplied with a uniform
	847	scalar, this is used to calculate the increment.
	848
	849	#+NAME: fdelay_constants
	850	#+BEGIN_SRC c
	851	#define FRACSCALE 0x10000000
	852	#+END_SRC
	853
	854	=FRACMASK= is the bitmask used to keep the fractional
	855	position in 28-bit range. It is $1^28 - 1$, or =0xFFFFFFF=.
	856	This is in particular is used to filter out upper bits that
	857	get carried over to the integer read position.
	858
	859	#+NAME: fdelay_constants
	860	#+BEGIN_SRC c
	861	#define FRACMASK 0xFFFFFFF
	862	#+END_SRC
	863
	864	=FRACNBITS= is the number of bits in the number. set to be
	865	28.
	866
	867	#+NAME: fdelay_constants
	868	#+BEGIN_SRC c
	869	#define FRACNBITS 28
	870	#+END_SRC
	871	*** Jitter
	872	Jitter in this context, is a random line segment
	873	generator. It linearly interpolates between random values in
	874	a given range, using random durations in a given range.
	875
	876	A line generator stores a counter and increment amount.
	877
	878	#+NAME: bigverb_delay
	879	#+BEGIN_SRC c
	880	int inc;
	881	int counter;
	882	#+END_SRC
	883
	884	#+NAME: init_jitter
	885	#+BEGIN_SRC c
	886	d->inc = 0;
	887	d->counter = 0;
	888	#+END_SRC
	889
	890	The most significant thing to happen in the jitter signal is
	891	calculating the next random segment. This is done in a
	892	static function called =generate_next_line=.
	893
	894	#+NAME: static_funcdefs
	895	#+BEGIN_SRC c
	896	static void generate_next_line(sk_bigverb_delay *d, int sr);
	897	#+END_SRC
	898
	899	To begin, another random value is created based on the
	900	previous random value.
	901
	902	The RNG algorithm used is quite simple, and is used to
	903	produce a 16-bit value.
	904
	905	$$x(n) = (5^{6}x(n - 1) + 1) \& 0xFFFF$$
	906
	907	Before and after this equation, the value is balanced so that
	908	is a 16-bit bipolar signal.
	909	Before, it adds 0x10000 if the
	910	value is less than 0. After, it substracts 0x10000 if the
	911	seed value is greater than 0x8000.
	912
	913	#+NAME: generate_random_number
	914	#+BEGIN_SRC c
	915	if (d->rng < 0) d->rng += 0x10000;
	916	/* 5^6 = 15625 */
	917	d->rng = (1 + d->rng * 0x3d09);
	918	d->rng &= 0xFFFF;
	919	if (d->rng >= 0x8000) d->rng -= 0x10000;
	920	#+END_SRC
	921
	922	This new random value is used to produce the next random
	923	value in seconds.
	924
	925	The line counter is reset. This value comes from the
	926	high-level parameter.
	927
	928	#+NAME: bigverb_delay
	929	#+BEGIN_SRC c
	930	int maxcount;
	931	#+END_SRC
	932
	933	#+NAME: init_jitter
	934	#+BEGIN_SRC c
	935	d->maxcount = round((sr / ((SKFLT)p->randfreq * 0.001)));
	936	#+END_SRC
	937
	938	#+NAME: reset_counter
	939	#+BEGIN_SRC c
	940	d->counter = d->maxcount;
	941	#+END_SRC
	942
	943	Compute delay time values. The current delay time,
	944	=curdel=, is obtained by subtracting the write + integer
	945	read positions, then adding in the fractional component.
	946	Wraparound is applied.
	947
	948	#+NAME: compute_delay_values
	949	#+BEGIN_SRC c
	950	curdel = d->wpos -
	951	(d->irpos + (d->frpos/(SKFLT)FRACSCALE));
	952	while (curdel < 0) curdel += d->sz;
	953	curdel /= sr;
	954	#+END_SRC
	955
	956	The next delay time to lerp to is derived from the RNG
	957	and drift amount.
	958
	959	#+NAME: compute_delay_values
	960	#+BEGIN_SRC c
	961	nxtdel = (d->rng * (d->drift * 0.0001) / 32768.0) + d->dels;
	962	#+END_SRC
	963
	964	The delay time, in seconds (=dels=).
	965
	966	#+NAME: bigverb_delay
	967	#+BEGIN_SRC c
	968	SKFLT dels;
	969	#+END_SRC
	970
	971	#+NAME: init_jitter
	972	#+BEGIN_SRC c
	973	d->dels = p->delay / 44100.0;
	974	#+END_SRC
	975
	976	#+NAME: bigverb_delay
	977	#+BEGIN_SRC c
	978	SKFLT drift;
	979	#+END_SRC
	980
	981	#+NAME: init_jitter
	982	#+BEGIN_SRC c
	983	d->drift = p->drift;
	984	#+END_SRC
	985
	986	The linear increment value is the difference between the current
	987	and next delay times, divided by the number of steps needed
	988	to draw a line between them (=counter=). This value is then
	989	converted into samples. An extra sample is tacked on to
	990	prevent nil values.
	991
	992	#+NAME: compute_increment
	993	#+BEGIN_SRC c
	994	inc = ((curdel - nxtdel) / (SKFLT)d->counter)*sr;
	995	inc += 1;
	996	#+END_SRC
	997
	998	This increment value is truncated and converted to the
	999	fractional read.
	1000
	1001	#+NAME: set_fractional_read
	1002	#+BEGIN_SRC c
	1003	d->inc = floor(inc * FRACSCALE);
	1004	#+END_SRC
	1005
	1006	#+NAME: funcs
	1007	#+BEGIN_SRC c
	1008	static void generate_next_line(sk_bigverb_delay *d, int sr)
	1009	{
	1010	SKFLT curdel;
	1011	SKFLT nxtdel;
	1012	SKFLT inc;
	1013	<<generate_random_number>>
	1014	<<reset_counter>>
	1015	<<compute_delay_values>>
	1016	<<compute_increment>>
	1017	<<set_fractional_read>>
	1018	}
	1019	#+END_SRC
	1020	*** Filter Memory
	1021	A one-pole lowpass filter such as the one used in the delay
	1022	line requires one sample of memory, which stores the output
	1023	of the previous filter. In a difference equation, this would
	1024	be known as $y(n - 1)$. In C code, we abbreviate this as
	1025	=y=.
	1026
	1027	#+NAME: bigverb_delay
	1028	#+BEGIN_SRC c
	1029	SKFLT y;
	1030	#+END_SRC
	1031
	1032	#+NAME: delay_init
	1033	#+BEGIN_SRC c
	1034	d->y = 0.0;
	1035	#+END_SRC

dcblocker.org+118−0

@@ -0,0 +1,118 @@
	1	#+TITLE: DC Blocker
	2	* Introduction
	3	This document outlines an implementation of a DC Blocker, a
	4	typical utility used to block out what is known as DC offset.
	5	DC Blockers are occasionally used to clean up recorded
	6	audio, but they are more often seen handling digital
	7	algorithms that tend to naturally introduce a lot of DC
	8	into their signal. This is common in the physical modelling
	9	family of DSP algorithms (Karplus-Strongs, waveguides,
	10	FDN reverberators, etc).
	11
	12	Without knowing how a DC blocker works, it would seem to have
	13	magical characteristics. A DC blocker has the ability to
	14	"drag down" a waveform of audio with DC offset so that it is
	15	centered. But how does it know how much to drag?
	16
	17	In truth, a DC blocker is typically nothing more than
	18	a recursive highpass filter whose cutoff set to be a
	19	sub-audio frequency. This has the effect of removing any
	20	slow moving waveforms that often are the culrpits of
	21	unwanted DC biasing.
	22	* Difference Equation
	23	The difference equation for a
	24	[[https://ccrma.stanford.edu/~jos/filters/DC_Blocker.html][DC blocker]] looks like this:
	25
	26	$$
	27	y(n) = x(n) - x(n - 1) + Ry(n - 1)
	28	$$
	29
	30	Where $R$ is usually a constant between 0.9 and 1.
	31
	32	This difference equation outlines the whole darn algorithm!
	33	What will be needed is a function that takes in an input,
	34	and stores the previous input, output, as well as the
	35	constant $R$. And that's it.
	36	* Tangled Files
	37	This program outputs two files: =dcblocker.c= and
	38	=dcblocker.h=. Using the macro =SK_DCBLOCKER_PRIV= will
	39	expose the main struct. Otherwise, it is opaque.
	40
	41	#+NAME: dcblocker.c
	42	#+BEGIN_SRC c :tangle dcblocker.c
	43	#define SK_DCBLOCKER_PRIV
	44	#include "dcblocker.h"
	45	<<funcs>>
	46	#+END_SRC
	47
	48	#+NAME: dcblocker.h
	49	#+BEGIN_SRC c :tangle dcblocker.h
	50	#ifndef SK_DCBLOCKER_H
	51	#define SK_DCBLOCKER_H
	52	#ifndef SKFLT
	53	#define SKFLT float
	54	#endif
	55	<<typedefs>>
	56	#ifdef SK_DCBLOCKER_PRIV
	57	<<structs>>
	58	#endif
	59	<<funcdefs>>
	60	#endif
	61	#+END_SRC
	62	* Struct Definition
	63	The main struct for this DC blocker will be called
	64	=sk_dcblocker=. It only needs three parameters: an
	65	=x= value storing the previous input value, a =y= value
	66	storing the previous output value, and an =R= value
	67	storing the DC filter constant.
	68
	69	#+NAME: typedefs
	70	#+BEGIN_SRC c
	71	typedef struct sk_dcblocker sk_dcblocker;
	72	#+END_SRC
	73	#+NAME: structs
	74	#+BEGIN_SRC c
	75	struct sk_dcblocker {
	76	SKFLT x, y, R;
	77	};
	78	#+END_SRC
	79	* Initialization
	80	The =sk_dcblocker= can be initialized with
	81	=sk_dcblocker_init=. It zeros out =x= and =y=, and sets
	82	=R= to be a value of =0.99=, which is a reasonable DC
	83	coefficent value.
	84
	85	#+NAME: funcdefs
	86	#+BEGIN_SRC c
	87	void sk_dcblocker_init(sk_dcblocker *dcblk);
	88	#+END_SRC
	89
	90	#+NAME: funcs
	91	#+BEGIN_SRC c
	92	void sk_dcblocker_init(sk_dcblocker *dcblk)
	93	{
	94	dcblk->x = 0;
	95	dcblk->y = 0;
	96	dcblk->R = 0.99; /* quite reasonable, indeed! */
	97	}
	98	#+END_SRC
	99	* Computation
	100	A single sample of audio is processed with the function
	101	=sk_dcblocker_tick=. It performs the arithmetic for the
	102	difference equation, updates the internal values, then
	103	returns the output.
	104
	105	#+NAME: funcdefs
	106	#+BEGIN_SRC c
	107	SKFLT sk_dcblocker_tick(sk_dcblocker *dcblk, SKFLT in);
	108	#+END_SRC
	109
	110	#+NAME: funcs
	111	#+BEGIN_SRC c
	112	SKFLT sk_dcblocker_tick(sk_dcblocker *dcblk, SKFLT in)
	113	{
	114	dcblk->y = in - dcblk->x + dcblk->R*dcblk->y;
	115	dcblk->x = in;
	116	return dcblk->y;
	117	}
	118	#+END_SRC

fmpair.org+576−0

@@ -0,0 +1,576 @@
	1	#+TITLE: FM Pair Oscillator
	2	* Introduction
	3	This document outlines an implementation of a classic
	4	FM oscillator with frequency, C:M ratio, and modulation
	5	index control.
	6
	7	There is quite a lot to be said about FM synthesis, but a
	8	classic FM oscillator implementation is poetic tribute to
	9	the technique. Requiring only two oscillators (usually sine
	10	waves), an FM operator pair is capable of generating very
	11	rich and dynamic spectrums. For how little computation is
	12	required, these suckers can pack a punch!
	13
	14	The FM pair is the Bechamel sauce of the FM world. Creating
	15	one from scratch provides great sonic intution into the
	16	mostly unintuitive world of FM synthesis programming. It can
	17	be the building block for grokking more complex modulation
	18	configurations, such as those found on the DX7.
	19	* Some words on FM synthesis
	20	(I can't help but ramble a little bit about the topic of
	21	FM synthesis. Feel free to scroll/skip/page ahead.)
	22
	23	FM synthesis, or Frequently Modulation synthesis, was first
	24	implemented by John Chowning while doing research at the
	25	Stanford Artificial Intelligence Laboratory (SAIL, whose
	26	spritual successor would become CCRMA, or the Center for
	27	Computer Research in Music and Acoustics). As the legend
	28	goes, Chowning stumbled upon the technique by accident.
	29	Experimenting with different vibrato speeds, he found that
	30	when vibrato got fast enough, it started changing timbre
	31	and spectrum.
	32
	33	The humble beginings of FM synthesis started out with an FM
	34	pair configuration similar to this one. But can you imagine
	35	the surprise he got when he first heard that new FM sound?
	36
	37	To paint a picture:
	38
	39	At the time, computers were about the
	40	size of fridges. They were expensive devices owned by the
	41	University, and time spent using the computer was precious.
	42	Chowing no doubt had to fight to get his time slot. Being
	43	a composer amongst scientists and engineers, he probably
	44	could only work in late hours.
	45
	46	Typing away at a mainframe 3am, Chowning would have to had
	47	to wait at least a day for the computer to run his FM
	48	synthesis program. Assuming it successfully ran without any
	49	errors, Chowning would then have to ship the output of that
	50	program off somewhere to get it converted to magnetic
	51	reel-to-reel tape. Back then,
	52	digital-to-analog converers (DACs) were things that you had
	53	to get in a car and drive to. Now-a-days, audio DACs are
	54	things found on pretty much every electronic device you own.
	55
	56	So, if you're keeping track, that's at least a week from
	57	code to sound. That's an enormous amount of round-trip
	58	latency! To go that long, and to finally hear these new
	59	sounds must have been such a rewarding experience.
	60
	61	Chowing's original FM synthesis algorithm would eventually
	62	become one of Stanford University's highest grossing patents
	63	of all time. To this day, the royalties of this patent help
	64	to sustain the CCRMA facilities.
	65
	66	Famously, the patent was sold to Yamaha, and that gave us
	67	things like the DX7 synthesizer (and all the cheesy music
	68	that came with it). Due to the efficiency of
	69	the algorithm, it also found itself in a whole generation of
	70	sound cards at the time, used in video game consoles and
	71	arcade machines. The most famous of these chips was the
	72	[[https://en.wikipedia.org/wiki/Yamaha_YM2151][YM2151]], also known as the OPM chip.
	73	* Overview of Algorithm and Implementation
	74	An FM Pair uses two oscillators, sometimes referred two as
	75	=operators=. One of these operators synthesizes audble sound.
	76	This known as the =carrier=, and can be remembered because
	77	it carries the sound signal to the speaker. The other
	78	operator modulates the frequency of the carrier, and it
	79	is known as the =modulator=.
	80
	81	If you set the modulator to any frequency that is audio rate
	82	(give or take 20hz or higher), you will technically get some
	83	kind of FM synthesis. However, depending on the frequency of
	84	the carrier, the spectrum can range from being very pitched
	85	with harmonically related overtones, to something quite
	86	clangorous and not at all something you could easily play
	87	"Hot Cross Buns" on. So, to make an FM pair easier to use,
	88	three unit-less parameters are used to control the
	89	FM characteristics of the sound that work evenly across
	90	the spectrum: a carrier-to-modulator frequency ratio known
	91	as the C:M ratio, and something known as the modulation
	92	index, which can kind of be thought of as a thing that
	93	controls the amount of modulation.
	94
	95	The C and M parameters are values that primarily multiply
	96	the frequency value of the carrier and modulator operators,
	97	respectively.
	98
	99	The modulation index can be defined as having the following
	100	relationship (taken from the Computer Music Tutorial's
	101	chapter on Modulation synthesis, page 229):
	102
	103	$$ I = D/M $$
	104
	105	Where $I$ is the modulation index, $M$ is the modulation
	106	frequency (which would be the oscillator frequency multiplied
	107	by the modulator ratio), and $D$ is the amount of modulation
	108	=depth=. The equation can be reworked in terms of $D$:
	109
	110	$$ D = MI $$
	111
	112	This depth value can be thought of as the "amplitude" amount
	113	the modulator operator (in Hz). This will produce a signal
	114	that modulates the carrier frequency up the depth amount,
	115	and then down that same amount.
	116
	117	In this implemetation, two table-lookup oscillators will be
	118	used for the operators, using an algorithm similar to the
	119	@!(ref "osc" "table lookup oscillator")!@, which employ a
	120	fixed-point phasor with linear interpolation. Details for
	121	the gory details of these oscillators are beyond the scope
	122	of this article, but are described in the implementation
	123	referenced above.
	124	* Tangled Files
	125	=fmpair.c= contains the ANSI C code required for the
	126	implementation.
	127
	128	#+NAME: fmpair.c
	129	#+BEGIN_SRC c :tangle fmpair.c
	130	#include <math.h>
	131	#define SK_FMPAIR_PRIV
	132	#include "fmpair.h"
	133	<<constants>>
	134	<<funcs>>
	135	#+END_SRC
	136
	137	=fmpair.h= is the corresponding header file for the C file.
	138	It contains all the forward declarations. If
	139	=SK_FMPAIR_PRIV= is defined, it exposes the structs.
	140
	141	#+NAME: fmpair.h
	142	#+BEGIN_SRC c :tangle fmpair.h
	143	#ifndef SK_FMPAIR_H
	144	#define SK_FMPAIR_H
	145	#ifndef SKFLT
	146	#define SKFLT float
	147	#endif
	148	<<typedefs>>
	149	<<funcdefs>>
	150	#ifdef SK_FMPAIR_PRIV
	151	<<structs>>
	152	#endif
	153	#endif
	154	#+END_SRC
	155	* Struct
	156	#+NAME: typedefs
	157	#+BEGIN_SRC c
	158	typedef struct sk_fmpair sk_fmpair;
	159	#+END_SRC
	160
	161	#+NAME: structs
	162	#+BEGIN_SRC c
	163	struct sk_fmpair {
	164	<<sk_fmpair>>
	165	};
	166	#+END_SRC
	167	* Table-lookup oscillator data and constants
	168	** Constants
	169	The constant =SK_FMPAIR_MAXLEN= is the maximum length a
	170	lookup table can be.
	171
	172	The constant =SK_FMPAIR_PHASEMASK= is a bitmask used by
	173	the fixed point phasor.
	174	#+NAME: constants
	175	#+BEGIN_SRC c
	176	#define SK_FMPAIR_MAXLEN 0x1000000L
	177	#define SK_FMPAIR_PHASEMASK 0x0FFFFFFL
	178	#+END_SRC
	179	** Struct Data
	180	The FM oscillator implements two independent table-[ookup
	181	oscillators. Each one needs their own set of wavetables,
	182	incrementors, sizes, and cached phase positions.
	183
	184	The convention used will this: any variables used for the
	185	carrier oscillator will begin with 'c':
	186
	187	#+NAME: sk_fmpair
	188	#+BEGIN_SRC c
	189	SKFLT *ctab;
	190	int csz;
	191	int clphs;
	192	#+END_SRC
	193
	194	#+NAME: init
	195	#+BEGIN_SRC c
	196	fmp->ctab = ctab;
	197	fmp->csz = msz;
	198	#+END_SRC
	199
	200	Any variables used for the modulator oscillator will begin
	201	with 'm':
	202
	203	#+NAME: sk_fmpair
	204	#+BEGIN_SRC c
	205	SKFLT *mtab;
	206	int msz;
	207	int mlphs;
	208	#+END_SRC
	209
	210	#+NAME: init
	211	#+BEGIN_SRC c
	212	fmp->mtab = mtab;
	213	fmp->msz = msz;
	214	#+END_SRC
	215
	216	In order to be used by the fixed point phasor,
	217	the initial phases need to be rescaled by =SK_FMPAIR_MAXLEN=
	218	and then converted to an integer value. These are
	219	stored in each respective variable keeping track
	220	of phase.
	221
	222	#+NAME: init
	223	#+BEGIN_SRC c
	224	fmp->clphs = floor(ciphs * SK_FMPAIR_MAXLEN);
	225	fmp->mlphs = floor(miphs * SK_FMPAIR_MAXLEN);
	226	#+END_SRC
	227
	228	A number of constants are used with these tables: =nlb=,
	229	=inlb=, =mask=, and =maxlens=. These are beyond the scope
	230	of this document, but are touched upon in @!(ref "osc")!@.
	231	Since these are all based on table size, two sets of
	232	constants are derived for the carrier and modulator.
	233
	234	#+NAME: sk_fmpair
	235	#+BEGIN_SRC c
	236	/* carrier constants */
	237	int cnlb;
	238	SKFLT cinlb;
	239	unsigned long cmask;
	240
	241	/* modulator constants */
	242	int mnlb;
	243	SKFLT minlb;
	244	unsigned long mmask;
	245
	246	SKFLT maxlens;
	247	#+END_SRC
	248
	249	#+NAME: init
	250	#+BEGIN_SRC c
	251	{
	252	int tmp;
	253
	254	/* carrier */
	255	tmp = SK_FMPAIR_MAXLEN / csz;
	256	fmp->cnlb = 0;
	257	while (tmp >>= 1) fmp->cnlb++;
	258
	259	/* modulator */
	260	tmp = SK_FMPAIR_MAXLEN / msz;
	261	fmp->mnlb = 0;
	262	while (tmp >>= 1) fmp->mnlb++;
	263	}
	264
	265	/* phase mask for dividing lower/upper bits */
	266
	267	fmp->cmask = (1<<fmp->cnlb) - 1;
	268	fmp->mmask = (1<<fmp->mnlb) - 1;
	269
	270	/* constant used to convert to floating point */
	271
	272	fmp->cinlb = 1.0 / (1<<fmp->cnlb);
	273	fmp->minlb = 1.0 / (1<<fmp->mnlb);
	274
	275	/* max table length in seconds */
	276	/* used to convert cycles-per-second units to cycles */
	277
	278	fmp->maxlens = 1.0 * SK_FMPAIR_MAXLEN / sr;
	279	#+END_SRC
	280	* Initialization
	281	The FM oscillator is initialize with =sk_fmpair_init=.
	282	It needs the sampling rate =sr=, and two sets of
	283	wavetables (=cwt= and =mwt=),
	284	their sizes (=csz=) and (=msz=), and
	285	initial phases (=ciphs= and =miphs=) for the carrier (c)
	286	and modulator (m) oscillators.
	287
	288	In a classic FM pair, these would be identical,
	289	and the wavetable would contain a sampled period of a
	290	sine wave.
	291
	292	#+NAME: funcdefs
	293	#+BEGIN_SRC c
	294	void sk_fmpair_init(sk_fmpair *fmp, int sr,
	295	SKFLT *ctab, int csz, SKFLT ciphs,
	296	SKFLT *mtab, int msz, SKFLT miphs);
	297	#+END_SRC
	298
	299	#+NAME: funcs
	300	#+BEGIN_SRC c
	301	void sk_fmpair_init(sk_fmpair *fmp, int sr,
	302	SKFLT *ctab, int csz, SKFLT ciphs,
	303	SKFLT *mtab, int msz, SKFLT miphs)
	304	{
	305	<<init>>
	306	}
	307	#+END_SRC
	308	* Parameter Control
	309	** Frequency
	310	Set with =sk_fmpair_freq=.
	311
	312	#+NAME: funcdefs
	313	#+BEGIN_SRC c
	314	void sk_fmpair_freq(sk_fmpair *fmp, SKFLT freq);
	315	#+END_SRC
	316
	317	#+NAME: funcs
	318	#+BEGIN_SRC c
	319	void sk_fmpair_freq(sk_fmpair *fmp, SKFLT freq)
	320	{
	321	fmp->freq = freq;
	322	}
	323	#+END_SRC
	324
	325	#+NAME: sk_fmpair
	326	#+BEGIN_SRC c
	327	SKFLT freq;
	328	#+END_SRC
	329
	330	A sensible default of A440.
	331
	332	#+NAME: init
	333	#+BEGIN_SRC c
	334	sk_fmpair_freq(fmp, 440);
	335	#+END_SRC
	336	** Carrier/Modulator Ratios
	337	Set with =sk_fmpair_modulator= and =sk_fmpair_carrier=.
	338
	339	#+NAME: funcdefs
	340	#+BEGIN_SRC c
	341	void sk_fmpair_modulator(sk_fmpair *fmp, SKFLT mod);
	342	void sk_fmpair_carrier(sk_fmpair *fmp, SKFLT car);
	343	#+END_SRC
	344
	345	#+NAME: funcs
	346	#+BEGIN_SRC c
	347	void sk_fmpair_modulator(sk_fmpair *fmp, SKFLT mod)
	348	{
	349	fmp->mod = mod;
	350	}
	351
	352	void sk_fmpair_carrier(sk_fmpair *fmp, SKFLT car)
	353	{
	354	fmp->car = car;
	355	}
	356	#+END_SRC
	357
	358	#+NAME: sk_fmpair
	359	#+BEGIN_SRC c
	360	SKFLT car;
	361	SKFLT mod;
	362	#+END_SRC
	363
	364	A 1:1 ratio is a good strong and sensible default.
	365
	366	#+NAME: init
	367	#+BEGIN_SRC c
	368	sk_fmpair_carrier(fmp, 1);
	369	sk_fmpair_modulator(fmp, 1);
	370	#+END_SRC
	371	** Modulation Index
	372	Set with =sk_fmpair_modindex=.
	373
	374	#+NAME: funcdefs
	375	#+BEGIN_SRC c
	376	void sk_fmpair_modindex(sk_fmpair *fmp, SKFLT index);
	377	#+END_SRC
	378
	379	#+NAME: funcs
	380	#+BEGIN_SRC c
	381	void sk_fmpair_modindex(sk_fmpair *fmp, SKFLT index)
	382	{
	383	fmp->index = index;
	384	}
	385	#+END_SRC
	386
	387	#+NAME: sk_fmpair
	388	#+BEGIN_SRC c
	389	SKFLT index;
	390	#+END_SRC
	391
	392	A modulation index of 1 on with 1:1 C-to-M ratio produces
	393	a very warm and versatile initial sound.
	394
	395	#+NAME: init
	396	#+BEGIN_SRC c
	397	sk_fmpair_modindex(fmp, 1);
	398	#+END_SRC
	399	* Computing a Sample
	400	A single sample of audio is computed with =sk_fmpair_tick=.
	401
	402	#+NAME: funcdefs
	403	#+BEGIN_SRC c
	404	SKFLT sk_fmpair_tick(sk_fmpair *fmp);
	405	#+END_SRC
	406
	407	#+NAME: funcs
	408	#+BEGIN_SRC c
	409	SKFLT sk_fmpair_tick(sk_fmpair *fmp)
	410	{
	411	SKFLT out;
	412	SKFLT cfreq, mfreq;
	413	SKFLT modout;
	414	int ipos;
	415	SKFLT frac;
	416	SKFLT x[2];
	417	out = 0;
	418	<<calculate_frequencies>>
	419	<<table_lookup_for_modulator>>
	420	<<scale_modulator_output>>
	421	<<modulate_carrier_frequency>>
	422	<<table_lookup_for_carrier>>
	423	<<update_phase_positions>>
	424	return out;
	425	}
	426	#+END_SRC
	427
	428	The carrier and modulator frequencies are calculated by
	429	multiplying their respective ratio values with
	430	the oscillator frequency.
	431
	432	#+NAME: calculate_frequencies
	433	#+BEGIN_SRC c
	434	cfreq = fmp->freq * fmp->car;
	435	mfreq = fmp->freq * fmp->mod;
	436	#+END_SRC
	437
	438	The modulator oscillator does its table-lookup and
	439	linear interpolation with some bitwise magic. The details of
	440	this can be found in @!(ref "osc")!@.
	441
	442	#+NAME: table_lookup_for_modulator
	443	#+BEGIN_SRC c
	444	fmp->mlphs &= SK_FMPAIR_PHASEMASK;
	445	ipos = fmp->mlphs >> fmp->mnlb;
	446	x[0] = fmp->mtab[ipos];
	447
	448	if (ipos == fmp->msz - 1) {
	449	x[1] = fmp->mtab[0];
	450	} else {
	451	x[1] = fmp->mtab[ipos + 1];
	452	}
	453
	454	frac = (fmp->mlphs & fmp->mmask) * fmp->minlb;
	455	modout = (x[0] + (x[1] - x[0]) * frac);
	456	#+END_SRC
	457
	458	Before it can modulate the carrier, the output of the
	459	modulator needs to be appropriately scaled. This amplitude
	460	amount is the modulator frequency and modulation
	461	index multiplied together.
	462
	463	#+NAME: scale_modulator_output
	464	#+BEGIN_SRC c
	465	modout = mfreq fmp->index;
	466	#+END_SRC
	467
	468	The frequency of the carrier oscillator is now ready to be
	469	modulated. Modulation is an addition operation; The output
	470	of the scaled modulator is tacked on the carrier frequency.
	471
	472	#+NAME: modulate_carrier_frequency
	473	#+BEGIN_SRC c
	474	cfreq += modout;
	475	#+END_SRC
	476
	477	A table-lookup operation computation happens for the
	478	carrier oscillator, similar to what happened with the
	479	modulator oscillator. This output is what gets returned.
	480
	481	#+NAME: table_lookup_for_carrier
	482	#+BEGIN_SRC c
	483	fmp->clphs &= SK_FMPAIR_PHASEMASK;
	484	ipos = (fmp->clphs) >> fmp->cnlb;
	485	x[0] = fmp->ctab[ipos];
	486
	487	if (ipos == fmp->csz - 1) {
	488	x[1] = fmp->ctab[0];
	489	} else {
	490	x[1] = fmp->ctab[ipos + 1];
	491	}
	492
	493	frac = (fmp->clphs & fmp->cmask) * fmp->cinlb;
	494	out = (x[0] + (x[1] - x[0]) * frac);
	495	#+END_SRC
	496
	497	To wrap things up, the phase increments and positions
	498	of both oscillators are updated, based on their respective
	499	frequencies.
	500
	501	increment value is derived by multiplying oscillator
	502	frequency (in units of cycles-per-second) to the maximum
	503	length of the phasor (in units of seconds). The seconds unit
	504	cancels, and the resulting output is cycles. This value,
	505	truncated to an integer value, becomes the phase increment
	506	value.
	507
	508	#+NAME: update_phase_positions
	509	#+BEGIN_SRC c
	510	fmp->clphs += floor(cfreq * fmp->maxlens);
	511	fmp->mlphs += floor(mfreq * fmp->maxlens);
	512	#+END_SRC
	513	* C:M Ratio tips
	514	Some suggestions to get started with picking out good C,M
	515	and I parameters. For those starting out, these should help
	516	build some sonic intuition in FM.
	517
	518	Generally speaking, the more rational C:M ratios are
	519	(1:1, 1:2, 3:2, etc), the more harmonic and pitched the
	520	spectrum will be. The more irrational they are
	521	(1:1.6180339..., 1:3.1415926..., etc) the more clangorous
	522	and unpitched the spectrum will be.
	523
	524	Fractional or strange looking ratios that look irrational
	525	can sometimes be simple ratios in disguise. For example,
	526	1:0.5 is harmonically similar to 2:1.
	527
	528	The carrier ratio value determines what the pitch will be.
	529	Whole integer values correspond with the harmonic series.
	530	A value of 1 is the base frequency, 2 is an octave, 3 is
	531	an octave and a fifth, etc.
	532
	533	The modulator ratio can be thought of as how spread-out
	534	the harmonic spectrum will be. A 1:1 ratio will give you
	535	a denser spectrum than a 1:7 ratio. These spread out
	536	spectrums were what made classic DX7 keyboard tine sounds
	537	work so well. Using ratios with high M values, they were
	538	able to cut through a mix while also leaving room for other
	539	instruments.
	540
	541	Adding small fractional values to the M ratio (ex: 1:1.007) can
	542	add some very unique spectral beating. FM pairs can be
	543	layered together with different fractional amounts to create
	544	thicker sounds.
	545
	546	In an FM pair, the behavior of the modulation index can be
	547	thought of as a kind of tone control, similar to a filter
	548	cutoff control that you'd find in a subtractive synthesizer.
	549	It's a gross simplification for all the amazing things
	550	happening to the spectrum, but it's close enough.
	551
	552	Controlling the modulation index of a 1:1 FM pair kind of
	553	feels like controlling a subtractive synthesizer, especially
	554	for lower modulation index values.
	555
	556	FM sounds can alias like crazy, especially when frequency,
	557	modulator, and modulation index values are high. When you
	558	crank up the modulation index crazy high, you can get some
	559	very interesting alias noises, but finding the sweet spots
	560	can be a bit of a treasure hunt. The plugin
	561	[[https://www.kvraudio.com/product/mr-alias-by-insert-piz-here][Mr.Alias]] is built around this concept.
	562
	563	Keyboard/frequency scaling the modulation index is a Good
	564	Idea. As frequencies get lower, the modulation index can
	565	get larger without risk of aliasing. More modulation index
	566	means more harmonic content, means more sonic beef for
	567	lower notes.
	568
	569	The modulation index spectrum behavior is determined by
	570	[[https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_first_kind][Bessel Functions of the First kind]].
	571	The wiki page on this has a good chart of this that
	572	plots the amplitudes first 3 harmonics. You can see how
	573	ahey dip in and out. One should at least be passively aware
	574	that this happens. Changing the modulation index on a bass
	575	sound, for example, can sometimes cause the fundamental
	576	to drop out, which can produces thin patches of sound.

index.org+21−0

@@ -0,0 +1,21 @@
	1	#+TITLE: sndkit
	2	sndkit is a sonic toolkit for everywhere.
	3
	4	It is a collection of highly portable ANSI-C implementations
	5	of audio DSP algorithms commonly used in computer music, all
	6	wrapped up inside of a statically generated wiki.
	7
	8	Each DSP algorithm is written in a literate programming
	9	style, a technique pioneered by Donald Knuth. Surrounding
	10	the implementation itself are words and equations used to
	11	help explain and contextualize what every single line of
	12	code is doing.
	13
	14	The text in this wiki is licensed under the
	15	[[https://creativecommons.org/share-your-work/public-domain/cc0/][Creatives Commons Zero]] license, otherwise
	16	known as "CC0".
	17
	18	The tangled output is dual-licensed under MIT or
	19	the [[https://unlicense.org/][Unlicense]].
	20
	21	Happy Hacking Humans.

modalres.org+306−0

@@ -0,0 +1,306 @@
	1	#+TITLE: ModalRes
	2	* Overview
	3	ModalRes is an implementation of a 2nd-order resonator
	4	filter suitable for use in modal synthesis.
	5
	6	This filter is particularly useful in constructing bell or
	7	other impact-based sounds. In order to do this, you'd need
	8	to know the modal ratios pre-measured of time. This is what
	9	is known as black-box modelling. A similar example of this
	10	is formant-based speech synthesis.
	11
	12	Some good modal ratios and measurements can be found [[https://csound.com/docs/manual/MiscModalFreq.html][here]].
	13	* Tangled Files
	14	=modalres.c= and =modalres.h= are the two files that are
	15	tangled.
	16
	17	#+NAME: modalres.c
	18	#+BEGIN_SRC c :tangle modalres.c
	19	#include <math.h>
	20	#define SK_MODALRES_PRIV
	21	#include "modalres.h"
	22
	23	#ifndef M_PI
	24	#define M_PI 3.14159265358979323846
	25	#endif
	26
	27	<<funcs>>
	28	#+END_SRC
	29
	30	#+NAME: modalres.h
	31	#+BEGIN_SRC c :tangle modalres.h
	32	#ifndef SK_MODALRES_H
	33	#define SK_MODALRES_H
	34
	35	#ifndef SKFLT
	36	#define SKFLT float
	37	#endif
	38
	39	<<typedefs>>
	40	<<funcdefs>>
	41
	42	#ifdef SK_MODALRES_PRIV
	43	<<structs>>
	44	#endif
	45	#endif
	46	#+END_SRC
	47	* Initiazation
	48	Initialized with =sk_modalres_init=. Only the sampling rate
	49	needs to be supplied.
	50
	51	#+NAME: funcdefs
	52	#+BEGIN_SRC c
	53	void sk_modalres_init(sk_modalres *mr, int sr);
	54	#+END_SRC
	55
	56	#+NAME: funcs
	57	#+BEGIN_SRC c
	58	void sk_modalres_init(sk_modalres *mr, int sr)
	59	{
	60	<<init>>
	61	}
	62	#+END_SRC
	63	* Struct
	64	The encapsulating struct is called =sk_modalres=.
	65
	66	#+NAME: typedefs
	67	#+BEGIN_SRC c
	68	typedef struct sk_modalres sk_modalres;
	69	#+END_SRC
	70
	71	#+NAME: structs
	72	#+BEGIN_SRC c
	73	struct sk_modalres {
	74	<<sk_modalres>>
	75	};
	76	#+END_SRC
	77	* Filter State Variables and Coefficients
	78	The filter state parameters are the typical things you'd
	79	need in a biquadratic filter difference equation.
	80
	81	Variables are created to store the previous x value and
	82	2 y values.
	83
	84	#+NAME: sk_modalres
	85	#+BEGIN_SRC c
	86	SKFLT x;
	87	SKFLT y[2];
	88	#+END_SRC
	89
	90	#+NAME: init
	91	#+BEGIN_SRC c
	92	mr->x = 0;
	93	mr->y[0] = 0;
	94	mr->y[1] = 0;
	95	#+END_SRC
	96
	97	Coefficients for the biquad are =b1= for the x value, and
	98	=a1= and =a2= for the y values. These coefficents get
	99	computed in the computation section.
	100
	101	#+NAME: sk_modalres
	102	#+BEGIN_SRC c
	103	SKFLT b1;
	104	SKFLT a1;
	105	SKFLT a2;
	106	#+END_SRC
	107
	108	#+NAME: init
	109	#+BEGIN_SRC c
	110	mr->b1 = 0;
	111	mr->a1 = 0;
	112	mr->a2 = 0;
	113	#+END_SRC
	114
	115	As per usual, a copy of the sampling rate is stored, as it
	116	is needed to calculate filter coefficients.
	117
	118	#+NAME: sk_modalres
	119	#+BEGIN_SRC c
	120	int sr;
	121	#+END_SRC
	122
	123	#+NAME: init
	124	#+BEGIN_SRC c
	125	mr->sr = sr;
	126	#+END_SRC
	127
	128	An unusual part of this filter design is the scaler
	129	variable, which is used to scale the audio output, but
	130	is not part of the actual transer function. We'll call
	131	this one =s=.
	132
	133	#+NAME: sk_modalres
	134	#+BEGIN_SRC c
	135	SKFLT s;
	136	#+END_SRC
	137
	138	#+NAME: init
	139	#+BEGIN_SRC c
	140	mr->s = 0;
	141	#+END_SRC
	142	* Parameters
	143	** Frequency
	144	Set with =sk_modalres_freq=.
	145
	146	#+NAME: funcdefs
	147	#+BEGIN_SRC c
	148	void sk_modalres_freq(sk_modalres *mr, SKFLT freq);
	149	#+END_SRC
	150
	151	#+NAME: funcs
	152	#+BEGIN_SRC c
	153	void sk_modalres_freq(sk_modalres *mr, SKFLT freq)
	154	{
	155	mr->freq = freq;
	156	}
	157	#+END_SRC
	158
	159	#+NAME: sk_modalres
	160	#+BEGIN_SRC c
	161	SKFLT freq;
	162	SKFLT pfreq;
	163	#+END_SRC
	164
	165	The frequency parameter uses caching to check if filter
	166	coefficients needs to be updated. The cached value is set
	167	to be negative in order to force computation on the first
	168	tick.
	169
	170	#+NAME: init
	171	#+BEGIN_SRC c
	172	sk_modalres_freq(mr, 1000);
	173	mr->pfreq = -1;
	174	#+END_SRC
	175	** Q value
	176	Set with =sk_modalres_q=.
	177
	178	#+NAME: funcdefs
	179	#+BEGIN_SRC c
	180	void sk_modalres_q(sk_modalres *mr, SKFLT q);
	181	#+END_SRC
	182
	183	#+NAME: funcs
	184	#+BEGIN_SRC c
	185	void sk_modalres_q(sk_modalres *mr, SKFLT q)
	186	{
	187	mr->q = q;
	188	}
	189	#+END_SRC
	190
	191	The Q parameter uses caching to check if filter
	192	coefficients needs to be updated. The cached value
	193	is set to be negative in order to force computation on
	194	the initial tick.
	195
	196	#+NAME: sk_modalres
	197	#+BEGIN_SRC c
	198	SKFLT q;
	199	SKFLT pq;
	200	#+END_SRC
	201
	202	#+NAME: init
	203	#+BEGIN_SRC c
	204	sk_modalres_q(mr, 1);
	205	mr->pq = -1;
	206	#+END_SRC
	207	* Computation
	208	A single sample is computed with =sk_modalres_tick=.
	209
	210	#+NAME: funcdefs
	211	#+BEGIN_SRC c
	212	SKFLT sk_modalres_tick(sk_modalres *mr, SKFLT in);
	213	#+END_SRC
	214
	215	#+NAME: funcs
	216	#+BEGIN_SRC c
	217	SKFLT sk_modalres_tick(sk_modalres *mr, SKFLT in)
	218	{
	219	SKFLT out;
	220
	221	out = 0;
	222	<<update_coefficients>>
	223	<<compute_difference_equation>>
	224	<<update_filter_state>>
	225	<<scale_output>>
	226	return out;
	227	}
	228	#+END_SRC
	229
	230	Update the coefficients, if needed.
	231	First, frequency is converted to radians.
	232	While traditionally represented as the
	233	lowercase $\omega$ in mathematical notation, the variable
	234	=w= will be used in C as crude approximation.
	235
	236	From there the following coefficients are computed:
	237
	238	#+BEGIN_SRC tex
	239	b1 = {1 \over \omega^2 + {\omega \over 2Q}}
	240	a1 = {1 - 2\omega^2}b1
	241	a2 = {\omega^2 - {\omega \over 2Q}}
	242	#+END_SRC
	243
	244	#+NAME: update_coefficients
	245	#+BEGIN_SRC c
	246	if (mr->freq != mr->pfreq \|\| mr->q != mr->pq) {
	247	SKFLT w;
	248	SKFLT a, b, d;
	249
	250	w = mr->freq * 2.0 * M_PI;
	251
	252	a = mr->sr / w;
	253	b = a*a;
	254	d = 0.5*a;
	255
	256	mr->pfreq = mr->freq;
	257	mr->pq = mr->q;
	258
	259	mr->b1 = 1.0 / (b + d/mr->q);
	260	mr->a1 = (1.0 - 2.0b) mr->b1;
	261	mr->a2 = (b - d/mr->q) * mr->b1;
	262	mr->s = d;
	263	}
	264	#+END_SRC
	265
	266	Compute difference equation. This computes the difference
	267	equation:
	268
	269	#+BEGIN_SRC tex
	270	$$
	271	y(n) = b0 x(n) + b1 x(n - 1) - a1 y(n - 1) - a2 y(n - 2)
	272	$$
	273	#+END_SRC
	274
	275	Where $b0$ is 0, effectively cancelling out this term.
	276
	277	What's interesting about this is that it adds an implicit
	278	one-sample delay to the filter.
	279	This is done because more accurately simulates the impulse
	280	response in an ideal mass-spring-damper system. The produces
	281	is, after all, a response to the impulse.
	282
	283	#+NAME: compute_difference_equation
	284	#+BEGIN_SRC c
	285	out = mr->b1mr->x - mr->a1mr->y[0] - mr->a2*mr->y[1];
	286	#+END_SRC
	287
	288	Update filter state. Shift things forward in time one
	289	sample. What was once =y[0]= is now =y[1]=. What is now
	290	the current output becomes =y[0]=, what is now the current
	291	input value because the previous input value =x=.
	292
	293	#+NAME: update_filter_state
	294	#+BEGIN_SRC c
	295	mr->y[1] = mr->y[0];
	296	mr->y[0] = out;
	297	mr->x = in;
	298	#+END_SRC
	299
	300	The output is scaled by the variable =s= before it is
	301	returned. It is not part of the recursive filter function.
	302
	303	#+NAME: scale_output
	304	#+BEGIN_SRC c
	305	out *= mr->s;
	306	#+END_SRC

osc.org+409−0

@@ -0,0 +1,409 @@
	1	#+TITLE: The Oscillator
	2	* Introduction.
	3	This document will describe an implementation of a classic
	4	table-lookup oscillator with linear interpolation.
	5
	6	The algorithm for this is an interesting mix of numerical
	7	processing, with the phasor being implemented in fixed
	8	point, and everything else being done in floating point.
	9	A big advantage to doing it this way is numerical stability:
	10	there is no risk of any phase acculumation or truncation
	11	like you'd get with floating-point. This
	12	is an important characteristic needed for phasor to avoid
	13	things like drift over time. As will be seen later, the
	14	fixed-point approach can be a little bit harder to
	15	understand, especially for those unfamiliar with fixed-point
	16	DSP. The author has put extra effort to explain these
	17	portions as clearly as possible.
	18
	19	The implementation also has a famous limitation of only
	20	being able to take in tables with sizes that are a power of
	21	2.
	22	* Generated Files
	23	Header:
	24
	25	#+NAME: osc.h
	26	#+BEGIN_SRC c :tangle osc.h
	27	#ifndef SK_OSC_H
	28	#define SK_OSC_H
	29	#ifndef SKFLT
	30	#define SKFLT float
	31	#endif
	32	<<typedefs>>
	33	#ifdef SK_OSC_PRIV
	34	<<structs>>
	35	#endif
	36	<<funcdefs>>
	37	#endif
	38	#+END_SRC
	39
	40	C file:
	41
	42	#+NAME: osc.c
	43	#+BEGIN_SRC c :tangle osc.c
	44	#include <stdint.h>
	45	#include <stdlib.h>
	46	#include <math.h>
	47	#define SK_OSC_PRIV
	48	#include "osc.h"
	49	<<constants>>
	50	<<funcs>>
	51	#+END_SRC
	52	* Top-level functions
	53	The table-lookup oscillator is initialized with
	54	=sk_osc_init=. The following arguments must be provided:
	55
	56	=osc= is a pre-allocated struct of =sk_osc=.
	57
	58	=sr= is the sampling rate.
	59
	60	=tab= is a pre-allocated wavetable, an array of =SKFLT=
	61	floating-point values.
	62
	63	=sz= is the array size of the wavetable =wt=.
	64
	65	=iphs= provides the initial phase of the oscillator. It is
	66	a value between 0-1.
	67
	68	#+NAME: funcdefs
	69	#+BEGIN_SRC c
	70	void sk_osc_init(sk_osc osc, int sr, SKFLT wt, int sz, SKFLT iphs);
	71	#+END_SRC
	72
	73	To compute a sample of audio, use =sk_osc_tick=.
	74
	75	#+NAME: funcdefs
	76	#+BEGIN_SRC c
	77	SKFLT sk_osc_tick(sk_osc *osc);
	78	#+END_SRC
	79
	80	This oscillator has 2 main parameters: frequency (freq),
	81	and amplitude (amp). They can be
	82	set with the following functions.
	83
	84	#+NAME: funcdefs
	85	#+BEGIN_SRC c
	86	void sk_osc_freq(sk_osc *osc, SKFLT freq);
	87	void sk_osc_amp(sk_osc *osc, SKFLT amp);
	88	#+END_SRC
	89
	90	#+NAME: funcs
	91	#+BEGIN_SRC c
	92	void sk_osc_freq(sk_osc *osc, SKFLT freq)
	93	{
	94	osc->freq = freq;
	95	}
	96
	97	void sk_osc_amp(sk_osc *osc, SKFLT amp)
	98	{
	99	osc->amp = amp;
	100	}
	101	#+END_SRC
	102	* Constants
	103	#+NAME: constants
	104	#+BEGIN_SRC c
	105	#define SK_OSC_MAXLEN 0x1000000L
	106	#define SK_OSC_PHASEMASK 0x0FFFFFFL
	107	#+END_SRC
	108	* Struct and Constants
	109	The main struct of this oscillator is called =sk_osc=.
	110	#+NAME: typedefs
	111	#+BEGIN_SRC c
	112	typedef struct sk_osc sk_osc;
	113	#+END_SRC
	114	#+NAME: structs
	115	#+BEGIN_SRC c
	116	struct sk_osc {
	117	<<sk_osc>>
	118	};
	119	#+END_SRC
	120
	121	#+NAME: sk_osc
	122	#+BEGIN_SRC c
	123	SKFLT freq, amp;
	124	SKFLT *tab;
	125	int inc;
	126	size_t sz;
	127	uint32_t nlb;
	128	SKFLT inlb;
	129	uint32_t mask;
	130	SKFLT maxlens;
	131	int32_t lphs;
	132	#+END_SRC
	133
	134	The oscillator stores it's main parameters =freq= and =amp=
	135	as floating point parameters. They are set to be values
	136	440 and 0.2 by default.
	137
	138	#+NAME: osc_init
	139	#+BEGIN_SRC c
	140	osc->freq = 440;
	141	osc->amp = 0.2;
	142	#+END_SRC
	143
	144	A reference to the table is stored in the variable =tab=,
	145	with its size =sz=.
	146
	147	#+NAME: osc_init
	148	#+BEGIN_SRC c
	149	osc->tab = wt;
	150	osc->sz = sz;
	151	#+END_SRC
	152
	153	A table lookup oscillator indexes through the table using
	154	the increment rate stored in the integer value =inc=. This
	155	value can be positive or negative. Is is zeroed out at
	156	init-time.
	157
	158	#+NAME: osc_init
	159	#+BEGIN_SRC c
	160	osc->inc = 0;
	161	#+END_SRC
	162
	163	The variable =lphs= stores the phase position of the
	164	previous sample. The initial phase value =iphs= is
	165	multiplied with the the maximum table value, and then
	166	masked to keep values in range.
	167
	168	#+NAME: osc_init
	169	#+BEGIN_SRC c
	170	osc->lphs = ((int32_t)(iphs * SK_OSC_MAXLEN)) & SK_OSC_PHASEMASK;
	171	#+END_SRC
	172
	173	For the fixed point table-lookup, some constants are derived
	174	and stored.
	175
	176	Phasor position is stored by splitting the bits of an
	177	N-bit integer number into two parts. The upper bits store the
	178	integer portion, while the lower bits store fractional
	179	portion. The maximum number of bits is arbitrary, but the
	180	underlying architecture must be able to accomodate for the
	181	width. In this implementation, the phasor uses 28 bits
	182	inside of a 32-bit number. This implicitely means the largest
	183	value can be
	184
	185	Split in the phasor position is measured by counting the
	186	number of lower bits. This value is stored in the variable
	187	=nlb=. This value is calculated with the equation
	188
	189	@!
	190	#$$n = log2(M/s)$$
	191	(img "/res/wiki/osc/osc_log2.png")
	192	!@
	193
	194	Where =n= is the number of lower bits, =M= is the maximum
	195	wavetable size, and =s= is the size of the wavetable.
	196
	197	To calculate =nlb=, and hand-rolled =log2= function is
	198	created.
	199
	200	Smaller values of =s= mean more bits in the fractional
	201	component of the number.
	202
	203	#+NAME: osc_init
	204	#+BEGIN_SRC c
	205	{
	206	uint32_t tmp;
	207	tmp = SK_OSC_MAXLEN / sz;
	208	osc->nlb = 0;
	209	while (tmp >>= 1) osc->nlb++;
	210	}
	211	<<calculate_mask>>
	212	<<calculate_inlb>>
	213	<<calculate_maxlens>>
	214	#+END_SRC
	215
	216	The =mask= is the lower-bits masking variable. When an AND
	217	operation is used against this constant, it filters out all
	218	the upper bits, so only the lower bits can pass through.
	219	This constant is necessary for being able to extract the
	220	lower bits from the fixed-point phase value representation.
	221	This sort of value is known in the bit-twiddling world as a
	222	=mask=. In binary, all the lower bits up to the number of
	223	lower bits are set to be on, with the remaining bits set
	224	to be 0.
	225
	226	#+NAME: calculate_mask
	227	#+BEGIN_SRC c
	228	osc->mask = (1<<osc->nlb) - 1;
	229	#+END_SRC
	230
	231	The inverse of max lower bits value, or $1/2^{nlb}$, is
	232	stored as a constant. This cached value is used to replace
	233	and divide operation with a multiply operation, which has
	234	traditionally been a cheaper operation to do on a computer.
	235
	236	#+NAME: calculate_inlb
	237	#+BEGIN_SRC c
	238	osc->inlb = 1.0 / (1<<osc->nlb);
	239	#+END_SRC
	240
	241	The constant =maxlens= is the maximum table length in units
	242	of seconds. This is a value used to efficiently convert the
	243	frequency parameter to sample increment value.
	244
	245	#+NAME: calculate_maxlens
	246	#+BEGIN_SRC c
	247	osc->maxlens = 1.0 * SK_OSC_MAXLEN / sr;
	248	#+END_SRC
	249	* Initialization
	250	In addition to setting variables, the init function will
	251	also set the starting phase of oscillator.
	252	#+NAME: funcs
	253	#+BEGIN_SRC c
	254	void sk_osc_init(sk_osc osc, int sr, SKFLT wt, int sz, SKFLT iphs)
	255	{
	256	<<osc_init>>
	257	}
	258	#+END_SRC
	259	* Computation
	260	The meat of the algorithm is here. Here outlines the
	261	=tick= function, where a single sample of an oscillator
	262	is computer.
	263
	264	#+NAME: funcs
	265	#+BEGIN_SRC c
	266	SKFLT sk_osc_tick(sk_osc *osc)
	267	{
	268	SKFLT out;
	269	SKFLT fract;
	270	SKFLT x1, x2;
	271	int32_t phs;
	272	int pos;
	273
	274	out = 0;
	275	<<update_increment_amount>>
	276	<<lookup_values>>
	277	<<obtain_fractional_component>>
	278	<<interpolate_values>>
	279	<<update_the_state>>
	280	return out;
	281	}
	282	#+END_SRC
	283
	284	First, the increment amount =inc= is updated.
	285	The increment amount tells how much further to move the read
	286	pointer in the table. This increment amount is based on the
	287	current oscillator frequency =freq= and the variable
	288	=maxlens=. =lrintf= is used to round to the nearest integer.
	289
	290	This is kind of a baffling operation. How could multiplying
	291	the frequency by some arbitrary duration yield a phasor
	292	increment amount? And where is the sampling rate in all
	293	this?
	294
	295	The thing that throws everything off is this fixed point
	296	business. Things make a lot more sense if you wanted to
	297	do this the floating-point way.
	298
	299	Recall that a phasor is a repeating line that ramps from
	300	0 to 1 over a given period of time. If we call this period
	301	of time =t=, the increment value =I= needs to work so that
	302	=t * sr * I = 1=. In other words, it's the slope of a line
	303	discretised.
	304
	305	Linear slope is found using good-ol rise over run, change in
	306	value over change in time. Inverting the frequency =1/F=
	307	gives it's period length in seconds. Multiplying by the
	308	sampling rate =sr= converts that value to samples. This
	309	gives us =sr/f=, our change in time. Because of the
	310	normalized range, change in value is just 1. Putting it all
	311	together This gives us a slope of =1/(sr/f)=, or =f/sr=.
	312	That is the normalized increment value.
	313
	314	If we wanted this to work with our fixed point phasor,
	315	we'd need to scale it by the maximum length of the phasor.
	316	That looks like =(f/sr)maxlen)=, or =(f maxlen)/sr)=.
	317	The frequency =f= can be pulled out of the numerator to
	318	get =f * (maxlen/sr)=, which can be reduced to the operation
	319	similarly seen below of =f * maxlens=.
	320
	321	Before frequency units were measured in Hertz, they were
	322	measured in cycles-per-second. If you ever took a highschool
	323	chemistry or physics class, you may recall that units can
	324	"cancel out" one another like a fraction. When
	325	cycles-per-second (cycles / seconds) gets multiplied by
	326	a value in seconds, the seconds cancel. What you are left
	327	with is a unit called cycles.
	328
	329	#+NAME: update_increment_amount
	330	#+BEGIN_SRC c
	331	osc->inc = floor(lrintf(osc->freq * osc->maxlens));
	332	#+END_SRC
	333
	334	Look up values A =n= and B =n + 1= samples from wavetable.
	335	Perform table lookup. Both the current position, and it's
	336	neighor are needed.
	337	This position is found by looking at the upper bits of
	338	the current phase.
	339
	340	#+NAME: lookup_values
	341	#+BEGIN_SRC c
	342	phs = osc->lphs;
	343	pos = phs >> osc->nlb;
	344	x1 = osc->tab[pos];
	345	x2 = osc->tab[(pos + 1) % osc->sz];
	346	#+END_SRC
	347
	348	Now, it's time to interpolate between points A and B. This
	349	oscillator uses linear interpolation, which can be thought
	350	of as a crossfade between two values. The equation for
	351	linear interpolation is commonly shown as:
	352
	353	@!
	354	#$$
	355	#y = (1 - \alpha) x_1 + \alpha x_2
	356	#$$
	357	(img "/res/wiki/osc/osc_lerp.png")
	358	!@
	359
	360	Where $x_1$ and $x_2$ are two values, and $\alpha$ is a
	361	fractional value between 0 and 1. The $\alpha$ values
	362	determines the distribution balance of
	363	the two values. When $\alpha = 0$, it is entirely $x_1$, and
	364	when $\alpha = 1$, the value is $x_2$.
	365
	366	The fractional (alpha) value is obtained by taking the lower
	367	bits portion of the current fixed-point phase position, and
	368	normalizing it to be a floating-point value between 0 and 1.
	369
	370	#+NAME: obtain_fractional_component
	371	#+BEGIN_SRC c
	372	fract = (phs & osc->mask) * osc->inlb;
	373	#+END_SRC
	374
	375	There are now all the parts to do the interpolation. It
	376	turns out the equation above can be simplified further to
	377	shave off a multiply operation.
	378
	379	@!
	380	# $$
	381	# \eqalign{
	382	# y &= (1 - \alpha)x_1 + \alpha x_2 \cr
	383	# &= x_1 - \alpha x_1 + \alpha x_2 \cr
	384	# &= x_1 + \alpha x_2 - \alpha x_1 \cr
	385	# &= x_1 + (x_2 - x_1)\alpha
	386	# }
	387	# $$
	388	(img "/res/wiki/osc/osc_lerp_simp.png")
	389	!@
	390
	391	Which then gets translated to the following C code below.
	392	In this step, the output is also being scaled by the
	393	amplitude =amp=.
	394
	395	#+NAME: interpolate_values
	396	#+BEGIN_SRC c
	397	out = (x1 + (x2 - x1) * fract) * osc->amp;
	398	#+END_SRC
	399
	400	And now the sample has been computed! To wrap up, the
	401	internal phase amount =lphs= is updated and masked to
	402	prevent overflow.
	403
	404	#+NAME: update_the_state
	405	#+BEGIN_SRC c
	406	phs += osc->inc;
	407	phs &= SK_OSC_PHASEMASK;
	408	osc->lphs = phs;
	409	#+END_SRC

peakeq.org+358−0

@@ -0,0 +1,358 @@
	1	#+TITLE: PeakEQ
	2	* Overview
	3	PeakEQ is an implementation of a single
	4	[[https://ccrma.stanford.edu/~jos//filters/Peaking_Equalizers.html][peaking equalizer filter]], with sweepable
	5	frequency, bandwidth, and gain.
	6	When used in series, the PeakEQ filter can be used
	7	to construct parametric and graphic equalizers commonly
	8	found in DAWs.
	9	* Algorithm Details
	10	The filter design based on the 2nd-order equalization filter
	11	design by Reglia and Mitra, from the paper
	12	"Tunable Digital Frequency Response Equalization Filters",
	13	IEEE Trans. on Ac., Sp. and Sig Proc., 35 (1), 1987
	14
	15	This implementation takes the original filter design and
	16	converts it to a direct-form 2 filter configuration, which
	17	has better numerical stability.
	18	* Conversion to Direct Form 2
	19	We begin with the transfer function, described in the 1987
	20	paper:
	21
	22	#+BEGIN_SRC tex
	23	$$
	24	{
	25	\alpha + \beta(1 + \alpha)z^{-1} + z^{-2}
	26	\over
	27	1 + \beta(1 + \alpha)z^{-1} + \alpha z^{-2}
	28	}
	29	$$
	30	#+END_SRC
	31
	32	This has the following difference equation:
	33
	34	#+BEGIN_SRC tex
	35	$$
	36	y(n) = \alpha x(n) + \beta(1 + \alpha)x(n - 1) + x(n - 2)
	37	- \beta(1 + \alpha)z^{-1} y(n - 1) - \alpha y(n - 2)
	38	$$
	39	#+END_SRC
	40
	41	The
	42	[[https://ccrma.stanford.edu/~jos/filters/Direct_Form_II.html][Direct Form II Structure]]. can be
	43	represented in this generalized difference equation:
	44
	45	#+BEGIN_SRC tex
	46	$$
	47	v(n) = x(n) - a_1 v(n - 1) - a_2 v(n - 2)
	48
	49	y(n) = b_0 v(n) + b_1 v(n - 1) + b_2 v(n - 2)
	50	$$
	51	#+END_SRC
	52
	53	Converting to direct form II is a matter of putting
	54	coefficients in their correct places. Placing in the
	55	coefficients from the Reglia and Mitra design yields the
	56	following difference equation:
	57
	58	#+BEGIN_SRC tex
	59	$$
	60	v(n) = x(n) - \beta(1 + \alpha)v(n - 1) - \alpha v(n - 2)
	61
	62	y(n) = \alpha v(n) + \beta(1 + \alpha) v(n - 1) + v(n - 2)
	63	$$
	64	#+END_SRC
	65
	66	Which can be translated to C code.
	67	* Tangled Files
	68	Tangles to files =peakeq.c= and =peakeq.h=.
	69
	70	#+NAME: peakeq.c
	71	#+BEGIN_SRC c :tangle peakeq.c
	72	#include <math.h>
	73	#define SK_PEAKEQ_PRIV
	74	#include "peakeq.h"
	75
	76	#ifndef M_PI
	77	#define M_PI 3.14159265358979323846
	78	#endif
	79
	80	<<funcs>>
	81	#+END_SRC
	82
	83	When =SK_PEAKEQ_PRIV= is defined, it exposes the main
	84	struct.
	85
	86	#+NAME: peakeq.h
	87	#+BEGIN_SRC c :tangle peakeq.h
	88	#ifndef SK_PEAKEQ_H
	89	#define SK_PEAKEQ_H
	90
	91	#ifndef SKFLT
	92	#define SKFLT float
	93	#endif
	94
	95	<<typedefs>>
	96	<<funcdefs>>
	97
	98	#ifdef SK_PEAKEQ_PRIV
	99	<<structs>>
	100	#endif
	101	#endif
	102	#+END_SRC
	103	* Struct
	104	Everything is stored in a struct called =sk_peakeq=.
	105
	106	#+NAME: typedefs
	107	#+BEGIN_SRC c
	108	typedef struct sk_peakeq sk_peakeq;
	109	#+END_SRC
	110
	111	#+NAME: structs
	112	#+BEGIN_SRC c
	113	struct sk_peakeq {
	114	<<sk_peakeq>>
	115	};
	116	#+END_SRC
	117	* Init
	118	PeakEQ is initialized with =sk_peakeq_init=. The sampling
	119	rate is required here.
	120
	121	#+NAME: funcdefs
	122	#+BEGIN_SRC c
	123	void sk_peakeq_init(sk_peakeq *eq, int sr);
	124	#+END_SRC
	125
	126	#+NAME: funcs
	127	#+BEGIN_SRC c
	128	void sk_peakeq_init(sk_peakeq *eq, int sr)
	129	{
	130	<<init>>
	131	}
	132	#+END_SRC
	133	* Filter State Memory and Constants
	134	The direct-form 2 form can be thought of as a 2-pole filter,
	135	followed by a 2-zero filter.
	136
	137	The memory for the 2-pole filter is stored in a variable
	138	called =v=.
	139
	140	#+NAME: sk_peakeq
	141	#+BEGIN_SRC c
	142	SKFLT v[2];
	143	#+END_SRC
	144
	145	#+NAME: init
	146	#+BEGIN_SRC c
	147	eq->v[0] = 0;
	148	eq->v[1] = 0;
	149	#+END_SRC
	150
	151	The alpha and beta coefficients are stored as variables
	152	called =a= and =b=.
	153
	154	#+NAME: sk_peakeq
	155	#+BEGIN_SRC c
	156	SKFLT a;
	157	SKFLT b;
	158	#+END_SRC
	159
	160	#+NAME: init
	161	#+BEGIN_SRC c
	162	eq->a = 0;
	163	eq->b = 0;
	164	#+END_SRC
	165
	166	The sampling rate is required to calculate new coefficients
	167	on-the-fly. A copy is stored in this struct.
	168
	169	#+NAME: sk_peakeq
	170	#+BEGIN_SRC c
	171	int sr;
	172	#+END_SRC
	173
	174	#+NAME: init
	175	#+BEGIN_SRC c
	176	eq->sr = sr;
	177	#+END_SRC
	178	* Parameters
	179	** Frequency
	180	Set with =sk_peakeq_freq=. Used to set the center frequency
	181	of the filter, in units of Hz.
	182
	183	#+NAME: funcdefs
	184	#+BEGIN_SRC c
	185	void sk_peakeq_freq(sk_peakeq *eq, SKFLT freq);
	186	#+END_SRC
	187
	188	#+NAME: funcs
	189	#+BEGIN_SRC c
	190	void sk_peakeq_freq(sk_peakeq *eq, SKFLT freq)
	191	{
	192	eq->freq = freq;
	193	}
	194	#+END_SRC
	195
	196	This parameter uses caching to detect if the parameter
	197	changes.
	198
	199	#+NAME: sk_peakeq
	200	#+BEGIN_SRC c
	201	SKFLT freq;
	202	SKFLT pfreq;
	203	#+END_SRC
	204
	205	=pfreq= is set to be a negative value to force initial
	206	coefficient calculation.
	207
	208	#+NAME: init
	209	#+BEGIN_SRC c
	210	sk_peakeq_freq(eq, 1000);
	211	eq->pfreq = -1;
	212	#+END_SRC
	213	** Bandwidth
	214	Set with =sk_peakeq_bandwidth=. This parameter sets the
	215	bandwidth of the EQ filter, in units of Hz.
	216
	217	#+NAME: funcdefs
	218	#+BEGIN_SRC c
	219	void sk_peakeq_bandwidth(sk_peakeq *eq, SKFLT bw);
	220	#+END_SRC
	221
	222	#+NAME: funcs
	223	#+BEGIN_SRC c
	224	void sk_peakeq_bandwidth(sk_peakeq *eq, SKFLT bw)
	225	{
	226	eq->bw = bw;
	227	}
	228	#+END_SRC
	229
	230	This parameter uses caching in order to detect if the
	231	parameter changes.
	232
	233	#+NAME: sk_peakeq
	234	#+BEGIN_SRC c
	235	SKFLT bw;
	236	SKFLT pbw;
	237	#+END_SRC
	238
	239	Like with freq, =pbw= is set to be a negative value to force
	240	coefficient calculation at the initial computation.
	241
	242	#+NAME: init
	243	#+BEGIN_SRC c
	244	sk_peakeq_bandwidth(eq, 1000);
	245	eq->pbw = -1;
	246	#+END_SRC
	247	** Gain
	248	This sets the gain of the filter. Positive values will cause
	249	a boost. Negative values will create a cut.
	250
	251	#+NAME: funcdefs
	252	#+BEGIN_SRC c
	253	void sk_peakeq_gain(sk_peakeq *eq, SKFLT gain);
	254	#+END_SRC
	255
	256	#+NAME: funcs
	257	#+BEGIN_SRC c
	258	void sk_peakeq_gain(sk_peakeq *eq, SKFLT gain)
	259	{
	260	eq->gain = gain;
	261	}
	262	#+END_SRC
	263
	264	#+NAME: sk_peakeq
	265	#+BEGIN_SRC c
	266	SKFLT gain;
	267	#+END_SRC
	268
	269	#+NAME: init
	270	#+BEGIN_SRC c
	271	sk_peakeq_gain(eq, 1.0);
	272	#+END_SRC
	273	* Compute
	274	=sk_peakeq_tick=.
	275
	276	#+NAME: funcdefs
	277	#+BEGIN_SRC c
	278	SKFLT sk_peakeq_tick(sk_peakeq *eq, SKFLT in);
	279	#+END_SRC
	280
	281	#+NAME: funcs
	282	#+BEGIN_SRC c
	283	SKFLT sk_peakeq_tick(sk_peakeq *eq, SKFLT in)
	284	{
	285	SKFLT out;
	286	SKFLT v;
	287	SKFLT y;
	288	out = 0;
	289
	290	<<update_coefficients>>
	291	<<compute_difference_equations>>
	292	<<compute_gain>>
	293	<<update_filter_state>>
	294
	295	return out;
	296	}
	297	#+END_SRC
	298
	299	Update coefficients, if needed. This happens at init time
	300	or when any of the parameters change.
	301
	302	The coefficents computed are alpha and beta. The beta
	303	coefficient is the negative cosine of the center frequency
	304	in units of radians. The alpha coefficient is the expression
	305	=(1 - c) / (1 + c)=, where =c= is the tangent of the
	306	bandwidth, times =PI=, divided by the sampling rate. More
	307	details on the both coefficient derivations can be found in
	308	the original Reglia and Mitra paper.
	309
	310	#+NAME: update_coefficients
	311	#+BEGIN_SRC c
	312	if (eq->bw != eq->pbw \|\| eq->freq != eq->pfreq) {
	313	SKFLT c;
	314	eq->b = -cos(2 * M_PI * eq->freq / eq->sr);
	315	c = tan(M_PI * eq->bw / eq->sr);
	316	eq->a = (1.0 - c) / (1.0 + c);
	317
	318	eq->pbw = eq->bw;
	319	eq->pfreq = eq->freq;
	320	}
	321	#+END_SRC
	322
	323	Compute difference equations. First the 2-pole filter,
	324	followed by the 2-zero filter. The output of the 2-pole
	325	filter goes into the 2-zero filter. See the previous section
	326	on conversion to direct-form 2 for more information.
	327
	328	#+NAME: compute_difference_equations
	329	#+BEGIN_SRC c
	330	v = in - eq->b(1.0 + eq->a)eq->v[0] - eq->a*eq->v[1];
	331	y = eq->av + eq->b(1.0 + eq->a)*eq->v[0] + eq->v[1];
	332	#+END_SRC
	333
	334	Compute the gain. This is done with the following equation:
	335
	336	#+BEGIN_SRC tex
	337	$$
	338	y = {((x + f) + g (x - f) \over 2}
	339	$$
	340	#+END_SRC
	341
	342	#+NAME: compute_gain
	343	#+BEGIN_SRC c
	344	out = ((in + y) + eq->gain(in - y)) 0.5;
	345	#+END_SRC
	346
	347	Where $y$ is the output, $x$ is the input signal, $g$ is
	348	the gain amount, and =f= is the filtered version of =x=.
	349
	350	Update filter state. The output of the 2-pole filter becomes
	351	=v[0]=, or $v(n - 1)$, and the previous =v[0]= becomes
	352	=v[1]=, or $v(n - 2)$.
	353
	354	#+NAME: update_filter_state
	355	#+BEGIN_SRC c
	356	eq->v[1] = eq->v[0];
	357	eq->v[0] = v;
	358	#+END_SRC

phasewarp.org+104−0

@@ -0,0 +1,104 @@
	1	#+TITLE: PhaseWarp
	2	* Overview
	3	PhaseWarp can be thought of as a special kind of filter for
	4	normalized phasor signals. When used with a phasor and a
	5	table-lookup oscillator this produces what is known as
	6	[[https://en.wikipedia.org/wiki/Phase_distortion_synthesis][phase distortion synthesis]], famously used on Casio
	7	synthesizers like the VZ-1.
	8	* Algorithm
	9	The algorithm for PhaseWarp works by adjusting the midpoint
	10	position where the phasor reaches 0.5. The adjustment amount
	11	is usually a normalized bipolar signal in the range -1 to 1.
	12	When the position value is positive, it warps the midpoint
	13	towards the upper bounds, causing a greater slope in the
	14	upper half, while decreasing the slope in the lower half.
	15	A negative position does opposite, creating a higher slope
	16	in the lower half, and a lower slope in the upper half.
	17
	18	(The whole thing makes a lot more sense when you graph it.
	19	Someday, I'll do that here. For now, you'll just have to
	20	take my word for it.)
	21	* Tangled Files
	22	=phasewarp.c= and =phasewarp.h=.
	23
	24	#+NAME: phasewarp.c
	25	#+BEGIN_SRC c :tangle phasewarp.c
	26	#include "phasewarp.h"
	27	<<funcs>>
	28	#+END_SRC
	29
	30	#+NAME: phasewarp.h
	31	#+BEGIN_SRC c :tangle phasewarp.h
	32	#ifndef SK_PHASEWARP_H
	33	#define SK_PHASEWARP_H
	34
	35	#ifndef SKFLT
	36	#define SKFLT float
	37	#endif
	38
	39	<<funcdefs>>
	40	#endif
	41	#+END_SRC
	42	* Compute
	43	This algorithm is so simple, it can be done in a stateless
	44	function called =sk_phasewarp_tick=. Its arguments are
	45	the input phasor signal and the warp factor.
	46
	47	#+NAME: funcdefs
	48	#+BEGIN_SRC c
	49	SKFLT sk_phasewarp_tick(SKFLT in, SKFLT warp);
	50	#+END_SRC
	51
	52	#+NAME: funcs
	53	#+BEGIN_SRC c
	54	SKFLT sk_phasewarp_tick(SKFLT in, SKFLT warp)
	55	{
	56	SKFLT out;
	57	SKFLT wmp;
	58
	59	out = 0;
	60
	61	<<calculate_warped_midpoint>>
	62	<<determine_side_and_warp>>
	63	return out;
	64	}
	65	#+END_SRC
	66
	67	Calculate warped midpoint. Call it =wmp=.
	68
	69	#+NAME: calculate_warped_midpoint
	70	#+BEGIN_SRC c
	71	wmp = (warp + 1.0) * 0.5;
	72	#+END_SRC
	73
	74	Determine side and compute. If the phasor position is
	75	less than the warped midpoint, than compute
	76	with lefthand slope. Otherwise, compute righthand slope.
	77
	78	#+NAME: determine_side_and_warp
	79	#+BEGIN_SRC c
	80	if (in < wmp) {
	81	<<compute_with_lefthand_slope>>
	82	} else {
	83	<<compute_with_righthand_slope>>
	84	}
	85	#+END_SRC
	86
	87	Compute With Lefthand Slope. This slope is equal to
	88	=0.5/wmp=. This value scales the phasor input =slope * in=.
	89
	90	#+NAME: compute_with_lefthand_slope
	91	#+BEGIN_SRC c
	92	if (wmp != 0) out = ((SKFLT)0.5 / wmp) * in;
	93	#+END_SRC
	94
	95	Compute Righthand Slope. The righthand slope is computed as
	96	=0.5/(1 - wmp)=. The output is computed as
	97	=slope * (in - wmp) + 0.5=.
	98
	99	#+NAME: compute_with_righthand_slope
	100	#+BEGIN_SRC c
	101	if (wmp != 1.0) {
	102	out = ((SKFLT)0.5 / (SKFLT)(1.0 - wmp)) * (in - wmp) + 0.5;
	103	}
	104	#+END_SRC

rline.org+403−0

@@ -0,0 +1,403 @@
	1	#+TITLE: R-Line
	2	* Introduction
	3	R-Line is a random line segment generator. It
	4	produces random values at regular rate, and then
	5	linearly interpolates between them. The algorithm has
	6	audio-rate parametric control of the range of values, as
	7	well as the frequency that values change.
	8
	9	This algorithm is a great little control signal generator
	10	that can add a bit of zest to a otherwise
	11	static parameter in a patch. It's a very efficient way to
	12	add a warm organic feel to a patch.
	13	* Notable Implementation Details
	14	This particular algorithm uses a fixed-point phasor similar
	15	to the ones seen in @!(ref "osc")!@ and @!(ref "fmpair")!@.
	16	Fixed point has very good numerical stability, which makes
	17	it ideal for a phasor algorithm.
	18
	19	To keep things more portable, R-Line implements its own
	20	pseudo-Random Number Generator (pRNG), using a common
	21	algorithm generating random values known as
	22	[[https://en.wikipedia.org/wiki/Linear_congruential_generator][Linear Congruential Generator]], or an LCG. This is
	23	what most =rand= function
	24	implementations use under the hood. The modulus, multiplier,
	25	and and increment values come from the Microsoft
	26	Visual/Quick C/C++ library (as found on Wikipedia).
	27	* Tangled Files
	28	R-Line tangles to a C file and a header file.
	29
	30	The C file =rline.c= contains the core implementation
	31	code.
	32
	33	#+NAME: rline.c
	34	#+BEGIN_SRC c :tangle rline.c
	35	#include <math.h>
	36	#define SK_RLINE_PRIV
	37	#include "rline.h"
	38	<<macros>>
	39	<<funcs>>
	40	#+END_SRC
	41
	42	The header file =rline.h= contains forward function
	43	declarations. When =SK_RLINE_PRIV= is enabled, the structs
	44	are exposed. Otherwise, they are left as opaque.
	45
	46	#+NAME: rline.h
	47	#+BEGIN_SRC c :tangle rline.h
	48	#ifndef SK_RLINE_H
	49	#define SK_RLINE_H
	50
	51	#ifndef SKFLT
	52	#define SKFLT float
	53	#endif
	54
	55	<<typedefs>>
	56	<<funcdefs>>
	57
	58	#ifdef SK_RLINE_PRIV
	59	<<structs>>
	60	#endif
	61	#endif
	62	#+END_SRC
	63	* Core Struct
	64	The core struct is called =sk_rline=.
	65
	66	#+NAME: typedefs
	67	#+BEGIN_SRC c
	68	typedef struct sk_rline sk_rline;
	69	#+END_SRC
	70
	71	#+NAME: structs
	72	#+BEGIN_SRC c
	73	struct sk_rline {
	74	<<sk_rline>>
	75	};
	76	#+END_SRC
	77	* Constants and Variables
	78	The value from the LCG is often scaled to be normalized
	79	between 0 and 1. The constant required for this is
	80	$1/(2^31 - 1$, or =1/((1L<<31) - 1)=. This value is
	81	stored in a constant.
	82
	83	#+NAME: sk_rline
	84	#+BEGIN_SRC c
	85	SKFLT rngscale;
	86	#+END_SRC
	87
	88	#+NAME: set_random_scaler
	89	#+BEGIN_SRC c
	90	rl->rngscale = 1.0 / ((1L<<31) - 1);
	91	#+END_SRC
	92
	93	The current state of the random number generator is stored
	94	as an integer =rng=. It is initialized to be =seed=.
	95
	96	#+NAME: sk_rline
	97	#+BEGIN_SRC c
	98	int rng;
	99	#+END_SRC
	100
	101	#+NAME: init
	102	#+BEGIN_SRC c
	103	rl->rng = seed;
	104	<<set_random_scaler>>
	105	<<generate_initial_values>>
	106	#+END_SRC
	107
	108	The fixed point phasor has two major constants: the maximum
	109	length, as well as the masking value. Detailed explanation
	110	of these are beyond the scope of this document.
	111
	112	#+NAME: macros
	113	#+BEGIN_SRC c
	114	#define SK_RLINE_PHSLEN 0x1000000L
	115	#define SK_RLINE_PHSMSK 0xFFFFFFL
	116	#+END_SRC
	117
	118	The variable =maxlens= is maximum phasor length converted to
	119	units of seconds by dviding it by the sampling rate.
	120	Multipyling this value with a frequency will return an
	121	increment value in units of phasor, which is used
	122	as the phasor increment value.
	123
	124	#+NAME: sk_rline
	125	#+BEGIN_SRC c
	126	SKFLT maxlens;
	127	#+END_SRC
	128
	129	#+NAME: init
	130	#+BEGIN_SRC c
	131	rl->maxlens = (SKFLT)SK_RLINE_PHSLEN / sr;
	132	#+END_SRC
	133
	134	Phase position is stored in an unsigned long called
	135	=phspos=. When it needs to be converted for use in floating
	136	point operations , the constant =scale= is used.
	137
	138	#+NAME: sk_rline
	139	#+BEGIN_SRC c
	140	unsigned long phasepos;
	141	SKFLT scale;
	142	#+END_SRC
	143
	144	#+NAME: init
	145	#+BEGIN_SRC c
	146	rl->phasepos = 0;
	147	#+END_SRC
	148
	149	The scale value works by normalizing to be in the range 0-1
	150	by dividing by the max phase length =SK_RLINE_PHSLEN=, then
	151	scaling it to be the difference between the start/end point
	152	values. Doing this shaves off a divide operation later.
	153
	154	#+NAME: calculate_initial_scale
	155	#+BEGIN_SRC c
	156	rl->scale = (rl->end - rl->start) / SK_RLINE_PHSLEN;
	157	#+END_SRC
	158
	159	A line has two points: a start point, and an end point.
	160	These are stored as normalized floating point variables
	161	=start= and =end=, and then are dynamically scaled to the
	162	=min= and =max= values during computation.
	163
	164	#+NAME: sk_rline
	165	#+BEGIN_SRC c
	166	SKFLT start;
	167	SKFLT end;
	168	#+END_SRC
	169
	170	#+NAME: generate_initial_values
	171	#+BEGIN_SRC c
	172	rl->rng = LCG(rl->rng);
	173	rl->start = RNG(rl->rng) * rl->rngscale;
	174	rl->rng = LCG(rl->rng);
	175	rl->end = RNG(rl->rng) * rl->rngscale;
	176	<<calculate_initial_scale>>
	177	#+END_SRC
	178	* Linear Congruential Generator
	179	An internal
	180	[[https://en.wikipedia.org/wiki/Linear_congruential_generator][Linear Congruential Generator]] is
	181	used to generate sequences of pseudo-random numbers.
	182
	183	It is the following equation
	184
	185	$$
	186	y(n) = ((y(n - 1) * a + c) \gg 1)
	187	\& m
	188	$$
	189
	190	Where $m$ is the masking value =0x7ffffff=, $a$ is a
	191	=multiplier=, and $c$ is the =increment=. In this
	192	implementation, $a$ will be =0x343fd=, and $c$ will be
	193	=0x2693ec3=. These constants come from the wikipedia page
	194	on LCGs, and originate the Microsoft Visual C/C++ standard
	195	library.
	196
	197	The LCG here can be implemented as a stateless function or
	198	macro. In this case, we will go with the macro.
	199
	200	#+NAME: macros
	201	#+BEGIN_SRC c
	202	#define LCG(y) (y * 0x343fd + 0x269ec3)
	203	#+END_SRC
	204
	205	The =LCG= operation only computes the next state state of
	206	the random-number generator. To actually get it within the
	207	correct bounds for this sytem, it has to be right-shifted
	208	to knock it down 1 bit, then masked by =0x7ffffff= as a kind
	209	of modulo operation.
	210
	211	This macro operation =RNG= assumes that =y= is the current
	212	state of the LCG.
	213	#+NAME: macros
	214	#+BEGIN_SRC c
	215	#define RNG(y) ((y >> 1) & 0x7fffffff)
	216	#+END_SRC
	217	* Initialization
	218	Initialization is done with =sk_rline_init=.
	219
	220	The main things needed for initialization are the sampling
	221	rate =sr=, as well as the initial seed value for the random
	222	number generator.
	223
	224	#+NAME: funcdefs
	225	#+BEGIN_SRC c
	226	void sk_rline_init(sk_rline *rl, int sr, int seed);
	227	#+END_SRC
	228
	229	#+NAME: funcs
	230	#+BEGIN_SRC c
	231	void sk_rline_init(sk_rline *rl, int sr, int seed)
	232	{
	233	<<init>>
	234	}
	235	#+END_SRC
	236	* Parameters
	237	** Minimum value
	238	Set with =sk_rline_min=.
	239
	240	#+NAME: funcdefs
	241	#+BEGIN_SRC c
	242	void sk_rline_min(sk_rline *rl, SKFLT min);
	243	#+END_SRC
	244
	245	#+NAME: funcs
	246	#+BEGIN_SRC c
	247	void sk_rline_min(sk_rline *rl, SKFLT min)
	248	{
	249	rl->min = min;
	250	}
	251	#+END_SRC
	252
	253	#+NAME: sk_rline
	254	#+BEGIN_SRC c
	255	SKFLT min;
	256	#+END_SRC
	257
	258	Initialized to be 0.
	259
	260	#+NAME: init
	261	#+BEGIN_SRC c
	262	sk_rline_min(rl, 0);
	263	#+END_SRC
	264	** Maximum value
	265	Set with =sk_rline_max=.
	266
	267	#+NAME: funcdefs
	268	#+BEGIN_SRC c
	269	void sk_rline_max(sk_rline *rl, SKFLT max);
	270	#+END_SRC
	271
	272	#+NAME: funcs
	273	#+BEGIN_SRC c
	274	void sk_rline_max(sk_rline *rl, SKFLT max)
	275	{
	276	rl->max= max;
	277	}
	278	#+END_SRC
	279
	280	#+NAME: sk_rline
	281	#+BEGIN_SRC c
	282	SKFLT max;
	283	#+END_SRC
	284
	285	Initialized to be 1.
	286
	287	#+NAME: init
	288	#+BEGIN_SRC c
	289	sk_rline_max(rl, 1);
	290	#+END_SRC
	291	** Rate
	292	Set with =sk_rline_rate=.
	293
	294	#+NAME: funcdefs
	295	#+BEGIN_SRC c
	296	void sk_rline_rate(sk_rline *rl, SKFLT rate);
	297	#+END_SRC
	298
	299	#+NAME: funcs
	300	#+BEGIN_SRC c
	301	void sk_rline_rate(sk_rline *rl, SKFLT rate)
	302	{
	303	rl->rate= rate;
	304	}
	305	#+END_SRC
	306
	307	#+NAME: sk_rline
	308	#+BEGIN_SRC c
	309	SKFLT rate;
	310	#+END_SRC
	311
	312	Initialized to be 1.
	313
	314	#+NAME: init
	315	#+BEGIN_SRC c
	316	sk_rline_rate(rl, 1);
	317	#+END_SRC
	318	* Computing a sample
	319	A single sample is computed with =sk_rline_tick=.
	320
	321	#+NAME: funcdefs
	322	#+BEGIN_SRC c
	323	SKFLT sk_rline_tick(sk_rline *rl);
	324	#+END_SRC
	325
	326	#+NAME: funcs
	327	#+BEGIN_SRC c
	328	SKFLT sk_rline_tick(sk_rline *rl)
	329	{
	330	SKFLT out;
	331
	332	out = 0;
	333
	334	<<compute_current_sample>>
	335	<<update_phase>>
	336	<<generate_next_line_segment>>
	337
	338	return out;
	339	}
	340	#+END_SRC
	341
	342	Compute the current sample. The line interpolation is
	343	calculated in a normalized space, then scaled to be
	344	in the min/max range. Doing it this way allows the min/max
	345	values to be dynamically changed over time without having
	346	to wait for the next line.
	347
	348	The normalized output can be computed with the expression:
	349
	350	$$
	351	y = x_1 + pc
	352	$$
	353
	354	Where $x_1$ is the starting point of the line, $p$ is
	355	the current phase increment, represented in fixed point, and
	356	$c$ is the constant that normalizes and scales the phase to
	357	be the amount of progress to value $x_2$.
	358
	359	#+NAME: compute_current_sample
	360	#+BEGIN_SRC c
	361	out = rl->start + rl->phasepos*rl->scale;
	362	out = out * (rl->max - rl->min) + rl->min;
	363	#+END_SRC
	364
	365
	366	Update phase position. The phase is updated by incrementing
	367	it by a amount obtained by multilying the frequency by the
	368	maximum phase length in units of seconds. How this works is
	369	beyond the scope of this document, but is explained in @!(ref
	370	"osc")!@.
	371
	372	#+NAME: update_phase
	373	#+BEGIN_SRC c
	374	rl->phasepos += floor(rl->rate * rl->maxlens);
	375	#+END_SRC
	376
	377	Generate next line segment. Preparation for a new line
	378	segment happens when the phase of the phasor reaches the
	379	end, and is greater than or equal to the max length.
	380	The phasor is masked in order to filter out upper
	381	bits and allow the lower bits to roll over. The the starting
	382	value is set to be the current end value, and a new end
	383	value is obtained using the random number generator.
	384
	385	After the new points have been obtained, the constant used
	386	to normalize + scale the phasor value is computed. Dividing
	387	by the maximum phasor length normalizes the phasor to be in
	388	range 0 and 1. Multiplying by the difference of the two
	389	segment values scales this value to be in the correct range.
	390	This constant is useful because it shaves off a divide
	391	operation, which has traditionally been costly arithmetic
	392	computatoin compared to a multiply.
	393
	394	#+NAME: generate_next_line_segment
	395	#+BEGIN_SRC c
	396	if (rl->phasepos >= SK_RLINE_PHSLEN) {
	397	rl->phasepos &= SK_RLINE_PHSMSK;
	398	rl->start = rl->end;
	399	rl->rng = LCG(rl->rng);
	400	rl->end = RNG(rl->rng) * rl->rngscale;
	401	rl->scale = (rl->end - rl->start) / SK_RLINE_PHSLEN;
	402	}
	403	#+END_SRC

scrubber.org+56−0

@@ -0,0 +1,56 @@
	1	#+TITLE: Scrubber
	2	* Overview
	3	An FFT based in-memory sample scrubber with pitch control.
	4	* Tangled Files
	5	* FFT procedures, and Soundpipe
	6	One could spend a whole lifetime building FFT algorithms,
	7	so an existing FFT implementation is used. For FFT and
	8	inverse-FFT operations, the Soundpipe fft utilities are
	9	used. So yes, Soundpipe is a hard dependency for this
	10	algorithm. Or at least, bits of Soundpipe are. The FFT
	11	portions can easily be extracted from the Soundpipe
	12	codebase.
	13
	14	The soundpipe FFT utility is a wrapper around a public
	15	domain FFT implementation by John Green.
	16	* Initialization
	17	Done with =sk_scrubber_init=.
	18	* Cleanup
	19	Done with =sk_scrubber_cleanup=. Mainly done for the sake
	20	of the =sp_fft=. Future iterations will try to remove any
	21	mandatory use of =malloc= and =free=.
	22	* Constants
	23	Scrubber uses a fixed FFT size, which greatly simplifies
	24	things. The FFT size chosen is 2048, or =2^11=. It's a good
	25	value that equal parts effeciency and fidelity.
	26
	27	The LOG2 of this FFT size is 11, and this is what is passed
	28	into the soundpipe FFT interface.
	29
	30	The hopsize determines how quickly to slide along the FFT
	31	window. In this case it's 512, or the FFT size divided up
	32	into 4 bits.
	33
	34	This divisor is a significant constant that will be
	35	stored as a constant as well.
	36	* Parameters
	37	** Position
	38	Sets the current read position of the buffer, in samples.
	39	** Playback
	40	Sets the relative playback speed. 1 is normal, 2 is 2x
	41	speed, and 0.5 is 2x slower.
	42	* Computation
	43	Read interpolated samples from buffer forward/backward
	44	frames.
	45
	46	Take FFT of frames
	47
	48	Phase Vocode.
	49
	50	Inverse FFT operation.
	51
	52	Overlap-add to output frames.
	53
	54	Write to output and scale;
	55
	56	Update counter and frame position.

vardelay.org+376−0

@@ -0,0 +1,376 @@
	1	#+TITLE: Vardelay
	2	* Overview
	3	Vardelay is an implementation of a variable-length delay
	4	line, using third order [[https://ccrma.stanford.edu/~jos/pasp/Lagrange_Interpolation.html][Langrange Interpolation]].
	5	This delay time of this delay line can be adjusted over
	6	time. Because it uses interpolation, any adjustments can
	7	cause pitch shifting to occur. This itself is often a
	8	desirable artifact, and is the basis for audio effects
	9	like flanging and chorusing.
	10
	11	In addition to delay time, a feedback parameter will also be
	12	supplied, as it is so often used with delay lines.
	13	* Tangled Files
	14	=vardelay.c= and =vardelay.h= are the tangled files.
	15
	16	#+NAME: vardelay.c
	17	#+BEGIN_SRC c :tangle vardelay.c
	18	#include <math.h>
	19	#include <stdlib.h>
	20	#define SK_VARDELAY_PRIV
	21	#include "vardelay.h"
	22	<<funcs>>
	23	#+END_SRC
	24
	25	If =SK_VARDELAY_PRIV= is defined, private structs are made
	26	public.
	27
	28	#+NAME: vardelay.h
	29	#+BEGIN_SRC c :tangle vardelay.h
	30	#ifndef SK_VARDELAY_H
	31	#define SK_VARDELAY_H
	32
	33	#ifndef SKFLT
	34	#define SKFLT float
	35	#endif
	36
	37	<<typedefs>>
	38	<<funcdefs>>
	39	#ifdef SK_VARDELAY_PRIV
	40	<<structs>>
	41	#endif
	42	#endif
	43	#+END_SRC
	44	* Initialization
	45	Initialized with =sk_vardelay_init=. must supply your own
	46	buffer, with buffer size. This way, no internal memory
	47	allocation is required.
	48
	49	#+NAME: funcdefs
	50	#+BEGIN_SRC c
	51	void sk_vardelay_init(sk_vardelay *vd, int sr,
	52	SKFLT *buf, unsigned long sz);
	53	#+END_SRC
	54
	55	#+NAME: funcs
	56	#+BEGIN_SRC c
	57	void sk_vardelay_init(sk_vardelay *vd, int sr,
	58	SKFLT *buf, unsigned long sz)
	59	{
	60	<<init>>
	61	}
	62	#+END_SRC
	63
	64	Anything buffer size less than 4 samples is not enough
	65	to do third-order interpolation, so it returns an error.
	66	* Struct
	67	Called =sk_vardelay=.
	68
	69	#+NAME: typedefs
	70	#+BEGIN_SRC c
	71	typedef struct sk_vardelay sk_vardelay;
	72	#+END_SRC
	73
	74	#+NAME: structs
	75	#+BEGIN_SRC c
	76	struct sk_vardelay {
	77	<<sk_vardelay>>
	78	};
	79	#+END_SRC
	80	* Constants and Variables
	81	** Sample rate
	82	a copy of the sampling rate is stored in the struct.
	83	It is used to convert delay time seconds to samples.
	84
	85	#+NAME: sk_vardelay
	86	#+BEGIN_SRC c
	87	int sr;
	88	#+END_SRC
	89
	90	#+NAME: init
	91	#+BEGIN_SRC c
	92	vd->sr = sr;
	93	#+END_SRC
	94	** Buffer, Buffer Size
	95	The delay buffer is the chunk of memory where samples will
	96	be read/written. It is assumed to be neatly pre-allocated
	97	and zeroed out externally. The size is also assumed to be
	98	known ahead of time and correct.
	99
	100	If the buffer size is less than 4, the delay line is
	101	disabled, as there aren't enough samples to do the
	102	interpolation.
	103
	104	#+NAME: sk_vardelay
	105	#+BEGIN_SRC c
	106	SKFLT *buf;
	107	unsigned long sz;
	108	#+END_SRC
	109
	110	#+NAME: init
	111	#+BEGIN_SRC c
	112	if (sz < 4) {
	113	vd->buf = NULL;
	114	vd->buf = 0;
	115	} else {
	116	vd->buf = buf;
	117	vd->sz = sz;
	118	}
	119	#+END_SRC
	120	** Previous output
	121	A variable is used to store the output of the previous delay
	122	in what is known as a single-sample delay or a unit delay.
	123	This delay is used to introduce feedback into the system.
	124	In discrete digital systems like this one, unit delays are
	125	inevitable for feedback, which is why they are referred to
	126	as implicit delays.
	127
	128	#+NAME: sk_vardelay
	129	#+BEGIN_SRC c
	130	SKFLT prev;
	131	#+END_SRC
	132
	133	#+NAME: init
	134	#+BEGIN_SRC c
	135	vd->prev = 0;
	136	#+END_SRC
	137	** Write Position
	138	The write position is used to store the current index
	139	position in the buffer being written to. The write position
	140	starts at zero, moves forward until it reaches the end of
	141	the buffer, and then goes back to the beginning again. The
	142	delayed samples start the write position, and move backwards
	143	in time. If that position is a negative, it wraps around.
	144
	145	Notably, the write position is a signed =long= to compensate
	146	for negative values. This is needed for the delay to read
	147	stuff back in time.
	148
	149	#+NAME: sk_vardelay
	150	#+BEGIN_SRC c
	151	long writepos;
	152	#+END_SRC
	153
	154	#+NAME: init
	155	#+BEGIN_SRC c
	156	vd->writepos = 0;
	157	#+END_SRC
	158	* Parameters
	159	** Delay Time
	160	Set with =sk_vardelay_delay=, In units of seconds.
	161
	162	#+NAME: funcdefs
	163	#+BEGIN_SRC c
	164	void sk_vardelay_delay(sk_vardelay *vd, SKFLT delay);
	165	#+END_SRC
	166
	167	#+NAME: funcs
	168	#+BEGIN_SRC c
	169	void sk_vardelay_delay(sk_vardelay *vd, SKFLT delay)
	170	{
	171	vd->delay = delay;
	172	}
	173	#+END_SRC
	174
	175	#+NAME: sk_vardelay
	176	#+BEGIN_SRC c
	177	SKFLT delay;
	178	#+END_SRC
	179
	180	The initial delay time is set to be half the delay time, in
	181	units of seconds.
	182
	183	#+NAME: init
	184	#+BEGIN_SRC c
	185	sk_vardelay_delay(vd, ((SKFLT)sz / sr) * 0.5);
	186	#+END_SRC
	187	** Feedback
	188	Set with =sk_vardelay_feedback=. Should be a value between
	189	0 and 1.
	190
	191	#+NAME: funcdefs
	192	#+BEGIN_SRC c
	193	void sk_vardelay_feedback(sk_vardelay *vd, SKFLT feedback);
	194	#+END_SRC
	195
	196	#+NAME: funcs
	197	#+BEGIN_SRC c
	198	void sk_vardelay_feedback(sk_vardelay *vd, SKFLT feedback)
	199	{
	200	vd->feedback = feedback;
	201	}
	202	#+END_SRC
	203
	204	#+NAME: sk_vardelay
	205	#+BEGIN_SRC c
	206	SKFLT feedback;
	207	#+END_SRC
	208
	209	No feedback by default.
	210
	211	#+NAME: init
	212	#+BEGIN_SRC c
	213	sk_vardelay_feedback(vd, 0);
	214	#+END_SRC
	215	* Computation
	216	Done with =sk_vardelay_tick=.
	217
	218	#+NAME: funcdefs
	219	#+BEGIN_SRC c
	220	SKFLT sk_vardelay_tick(sk_vardelay *vd, SKFLT in);
	221	#+END_SRC
	222
	223	#+NAME: funcs
	224	#+BEGIN_SRC c
	225	SKFLT sk_vardelay_tick(sk_vardelay *vd, SKFLT in)
	226	{
	227	SKFLT out;
	228	SKFLT dels;
	229	SKFLT f;
	230	long i;
	231	SKFLT s[4];
	232	unsigned long n[4];
	233	SKFLT a, b, c, d;
	234
	235	out = 0;
	236	<<return_if_empty>>
	237	<<write_to_buffer>>
	238	<<calculate_read_position>>
	239	<<wrap_and_flip>>
	240	<<read_samples>>
	241	<<calculate_interpolation_coefficients>>
	242	<<interpolate>>
	243	<<update_position>>
	244	<<update_feedback>>
	245
	246	return out;
	247	}
	248	#+END_SRC
	249
	250	If buffer is NULL or size is 0, return 0.
	251
	252	#+NAME: return_if_empty
	253	#+BEGIN_SRC c
	254	if (vd->buf == NULL \|\| vd->sz == 0) return 0;
	255	#+END_SRC
	256
	257	Write to buffer with feedback.
	258
	259	#+NAME: write_to_buffer
	260	#+BEGIN_SRC c
	261	vd->buf[vd->writepos] = in + vd->prev * vd->feedback;
	262	#+END_SRC
	263
	264	Calculate the read position. This is two parts: a fractional
	265	component, and an integer component.
	266
	267	Get floating point + integer components of delay time.
	268	The delay time position is set to be the write position
	269	=dels= samples in the past, where =dels= is the delay time
	270	in samples.
	271
	272	#+NAME: calculate_read_position
	273	#+BEGIN_SRC c
	274	dels = vd->delay * vd->sr;
	275	i = floor(dels);
	276	f = i - dels;
	277	i = vd->writepos - i;
	278	#+END_SRC
	279
	280	Because we're looking backwards, the fractional component is
	281	backwards too. A fractional value here is set to be
	282	negative (reaching back in time).
	283
	284	If there is a fractional component that is negative, or the
	285	integer component is negative, these must be corrected.
	286
	287	The fractional value is flipped to be positive by adding 1
	288	to itself. The integer position is set back in time one
	289	sample. This sets the interpolation up so that instead of
	290	taking a sample and interpolating backwards, you start with
	291	the previous sample and move forwards. The integer position
	292	is then set to be in bounds.
	293
	294	By default, the position is wrapped around using addition
	295	until the sample is in bounds.
	296
	297	#+NAME: wrap_and_flip
	298	#+BEGIN_SRC c
	299	if ((f < 0.0) \|\| (i < 0)) {
	300	/* flip fractional component */
	301	f = f + 1.0;
	302	/* go backwards one sample */
	303	i = i - 1;
	304	while (i < 0) i += vd->sz;
	305	} else while(i >= vd->sz) i -= vd->sz;
	306	#+END_SRC
	307
	308	Read samples. This includes the current sample, the previous
	309	sample, and two samples in the future.
	310
	311	#+NAME: read_samples
	312	#+BEGIN_SRC c
	313	/* x(n) */
	314	n[1] = i;
	315
	316	/* x(n + 1) */
	317	if (i == (vd->sz - 1)) n[2] = 0;
	318	else n[2] = n[1] + 1;
	319
	320	/* x(n - 1) */
	321	if (i == 0) n[0] = vd->sz - 1;
	322	else n[0] = i - 1;
	323
	324	if (n[2] == vd->sz - 1) n[3] = 0;
	325	else n[3] = n[2] + 1;
	326
	327	{
	328	int j;
	329	for (j = 0; j < 4; j++) s[j] = vd->buf[n[j]];
	330	}
	331	#+END_SRC
	332
	333	Calculate interpolation coefficients. These four
	334	coefficients correspond with the four samples read.
	335
	336	#+NAME: calculate_interpolation_coefficients
	337	#+BEGIN_SRC c
	338	{
	339	SKFLT tmp[2];
	340
	341	d = ((f * f) - 1) * 0.1666666667;
	342	tmp[0] = (f + 1.0) * 0.5;
	343	tmp[1] = 3.0 * d;
	344	a = tmp[0] - 1.0 - d;
	345	c = tmp[0] - tmp[1];
	346	b = tmp[1] - f;
	347	}
	348	#+END_SRC
	349
	350	Interpolate. This follows the following equation:
	351
	352	$$
	353	y(n) = (a x(n - 1) + b x(n) + c x(n + 1) + d x(n + 2)) \cdot f + x(n)
	354	$$
	355
	356	#+NAME: interpolate
	357	#+BEGIN_SRC c
	358	out = (as[0] + bs[1] + cs[2] + ds[3]) * f + s[1];
	359	#+END_SRC
	360
	361	Update position. Increment the write position, and wrap
	362	back to zero if it reaches the end of the delay buffer.
	363
	364	#+NAME: update_position
	365	#+BEGIN_SRC c
	366	vd->writepos++;
	367	if (vd->writepos == vd->sz) vd->writepos = 0;
	368	#+END_SRC
	369
	370	Update feedback. The current output is set to be the =prev=
	371	value in the delay line.
	372
	373	#+NAME: update_feedback
	374	#+BEGIN_SRC c
	375	vd->prev = out;
	376	#+END_SRC

importing some org files

Keyboard shortcuts

General

Issues and pull requests

Issues, pull requests and Inbox

Repository code