#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <math.h>

#define RATE 16000.0

float sinetab[2048];

typedef struct oscillator {
	uint32_t phasor;
	uint32_t step;
	float value;
} osc_t;

void update_osc(osc_t* o) {
	o->phasor += o->step;
	o->value = sinetab[o->phasor >> (32-11)];
}

void set_osc_freq(osc_t* o, float f) {
	o->step = (uint32_t) ((f / RATE) * UINT32_MAX);
}

typedef struct noise {
	float r1, r2, filt, out1, out2;
	uint8_t left;
} noise_t;

float get_noise(noise_t* n) {
	// See https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
	float in;
	if(n->left == 0) {
		n->r1 = (float) random() / UINT32_MAX;
		n->r2 = (float) random() / UINT32_MAX;
		n->left = 2;
	}
	if(n->left == 2) {
		in = (sqrt(-2*log(n->r1))*cos(2*M_PI*n->r2))/6;
	} else if(n->left == 1) {
		in = (sqrt(-2*log(n->r1))*sin(2*M_PI*n->r2))/6;
	}
	n->left--;
	// Two cascaded single order LPFs
	n->out1 += n->filt*(in - n->out1);
	n->out2 += n->filt*(n->out1 - n->out2);
	return n->out2;
}

int main() {
	osc_t lfo;
	noise_t noise;
	int16_t out;

	for(int i=0; i<2048; i++) {
		sinetab[i] = (float) sin(i*2*M_PI/2048);
	}
	set_osc_freq(&lfo, 0.1);
	noise.left = 0;

	while(1) {
		update_osc(&lfo);
		noise.filt = 0.5 + 0.5*lfo.value;
		out = INT16_MAX*get_noise(&noise);
		fwrite(&out, 2, 1, stdout);
	}

	return 1;
}