/*
Copyright 2021 Tyler de Witt (0xDEAFBEEF) www.deafbeef.com

Synth Poems are 1min generative audiovisual pieces using FM synthesis,
and XY oscilloscope visualization. Each entirely generated unsupervised
by a self contained code model, which is stored on chain.

Frequency Modulation(FM) is an encoding method originally used in radio
broadcasting, and later discovered to be useful for audio synthesis allowing a simple
way to generate rich harmonic timbres. The left and right stereo audio channels
are mapped to the XY inputs of a (virtual) oscilloscope to produce a direct
visualization of the audio waveform.

At the musical level, this generative algorithm plays with feedback to produce
an evolving but balanced set of timbres and pitches. Across pieces, a carefully
selected set of rules govern the variation of tempos, time signatures, pitches,
and (in)harmonicity resulting in individual pieces that are unique, while retaining
a common theme.

*/

char *seed="0x2c68ca657f166e1e7a5db2266ee1edca9daf2da5ecc689c39c1d8cfe0b9aad2e";

/*
This self contained C file will compile to run to generate audio and image data.
Replace the example seed value above with with the 'seed' field returned by
the smart contracts getTokenParams(tokenID) function, where tokenID is your unique
ERC721 tokenID.

Code conforms to the ISO c99 standard. Only library dependency is C math library.
Pass '-lm' to gcc. Assumes architecture is little endian.

1. Create directory to store output images:
# mkdir 'frames'
2. Compile and run:
# gcc main.c -lm && ./a.out

This will produce:
a. Numbered image files in BMP file format stored in 'frames' directory. 24 FPS(frames
per second)
b. Audio file named 'output.wav' in WAV format, stereo, 44.1khz, 16bit depth

This is the 'raw' information representing digital audio signals
and image pixel intensities. Based on simple standards. Lists of numbers. Future proof.

This raw information can be encoded perceptually using whatever tools exist at the
time of reconstruction. At present, for example, opensource tools such as imagemagick and ffmpeg
can be used to encode images,video and audio in different format for popular distributions.

Examples:

1. Audio can be encoded to MP3.
# ffmpeg -i output.wav -b 192k output.mp3

2. Images can be assembled into an animation and encoded to MP4, for example:
# ffmpeg -framerate 24 -i frames/frame%04d.bmp -crf 20 video.mp4

3. audio can be encoded with video:
# ffmpeg -i output.wav -i video.mp4 -c:v libx264 -c:a aac -b:a 192k -pix_fmt yuv420p -shortest audiovisual.mp4

4. Screenshots of particular frames(in this example, frame 60), using imagemagick:
# convert frames/frame0060.bmp -scale 720x720 image.png

*/

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string.h>

int myrand();
#define rand() myrand()

#ifndef M_PI
#define M_PI 3.14159265358979323846264338327950288
#endif

#define SAMPLE_RATE 44100

typedef struct {
double freq;
double phase;
double amp;
double pw; //pulse width for pulse generator
int i; // sample counter

// for DSF
double freq2;
double phase2;
double w;
} gen_data;

typedef struct {
int ftype;
double fc; //center freq
double bw; //bandwidth (for resonent filters)

double pole;
double a1;
double a2;
double rA; // for RES2A filter
double rB; // for RES2A filter
double rC; // for RES2A filter
double gain;
double x[3];
double y[3];

} iirfilter_params;
#define FTYPE_LP 1
#define FTYPE_HP 2
#define FTYPE_RES2 3
#define FTYPE_RES2A 4

double *make_bufn(int n);
double *make_buf(double dur,int*n);
void fill(double a, double *env,int n);
void delay(int b, double *data,int n);
void shift(int b, double *data,int n);
double maxsample(double *data,int n);
void normalize(double *data,int n);
void normalize_stereo(double *lchan, double *rchan,int n);
void mult(double *src, double *dst, int n);
void lerp(double a, double *src, double *dst, int n);
void add(double *src, double *dst, int n);
void scale(double a,double b, double *dst, int n);
void softclip(double gain,double offset,double *data,int n);
void gain(double g,double*data,int n);
double drand ( double low, double high );
double nond_drand ( double low, double high );
void dbgain(double db,double*data,int n);

//misc
double limit(double low,double high,double x);
double dmod(double x, double y);

void stereomixdown(double *lchan,double*rchan,int n);
void msmix(double *src,double dbgain, double pos, double *dstl, double *dstr,int n);

//bsets (breakpoint sets)
//breakpoint set
typedef struct bset {
int n; // number of samples used
int max_n; //max allocated samples
int sr; //sample rate
int interp; //default interpolation
double param; //default param
double *t; //pointer to samples
double *y; //pointer to samples
double *iparam; //parameter for interpolation technique (gamma, tau, etc);
int *itype; //type of interpolation for next segment
} bset;

#define INTERP_SAMPLEHOLD 1
#define INTERP_LINEAR 2
#define INTERP_CUBIC 3
#define INTERP_GAMMA 4
#define INTERP_EXP 5

void bset_init(bset *b);
void bset_destroy(bset *b);
int bset_search(bset *b, double t);
double bset_interp(bset *b, double t);
double bset_interp_linear(bset *b, double t);
double bset_interp_cubic(bset *b, double t);
double bset_interp_samplehold(bset *b, double t);
double bset_interp_exp(bset *b, double t,double tau);
double bset_interp_gamma(bset *b, double t,double gamma);
double bset_render(bset *b, double *data);
void bset_dump(bset *b);
void bset_write(bset *b,char *filename);
void bset_read(bset *b,char *filename);
void bset_add(bset *b, double t, double y);
void bset_segment(bset *b,double t1, double t2, double y1, double y2);
void bset_remove(bset *b, int j);
void bset_clear(bset *b);
void bset_filter_glide(bset *b, double dur,double dur2);
void bset_parse(bset *b, char *s);
void bset_parse2(bset *b, char *s);
void bset_parse_file(bset *b, char *filename,int voice);
void bset_stretch(bset *b,double f);
void bset_clip(bset *b,double t1,double dur);
void bset_shift(bset *b,double f);
void bset_delay(bset *b,double f);
double bset_length(bset *b);
double bset_last(bset *b);
void bset_dup(bset *b, double times, double gap);

//generators
void init_gen(gen_data *g);
double dsf_gen(gen_data *g);
double dsf(double fc, double fm, double w, double u, double v);
double sinlu(double w);
double sin_gen(gen_data *g);
void whitenoise(double gain, double *data,int n);

//filters
void lowpass(double fc,double *data,int n);
void highpass(double fc,double *data,int n);
void res2a(double fc,double bw, double *data,int n);
void res2a_bset(bset *fenv,
bset *bwenv,
double fmult,
double *data,int n);

//compression, gating, envelope detection
void envdetect(double attack, double release, double *src,double *dst,int n);
void compress(double ratio, double threshold, double *cs, double *src,double *dst,int n);
void compress_inplace(double attack, double release, double db, double *src,int n);
void compress_inplace_stereo(double attack, double release, double db, double *src0, double *src1, int n);

//effects
void to8bit(int num,double *data, int n);

void write_wav_header(FILE *f,int numSamples);
void write_wav_data(FILE *f,double *lchannel,double *rchannel,int num_samples);


// END HEADER

void gain(double g,double*data,int n) {
for (int i=0;i<n;i++) {
data[i] *= g;
}
}

void dbgain(double db,double*data,int n) {
double g = exp(db/20 * log(10));

for (int i=0;i<n;i++) {
data[i] *= g;
}
}

//TODO: have different streams of PRNG

#define RAND_MAX_P 2147483647

//LCD PRNG with 8 streams
unsigned long rs[8] = {1,2,3,4,5,6,7,8};
int seed_index = 0;
int randp_seed(int i,int entropy) {
rs[i] = entropy;
rs[i] = rs[i] % 2147483647;
if (rs[i] <= 0) rs[i] += 2147483646;
}
int randp_getstream() { return seed_index; }
void randp_setstream(int i) {
if (i>=0 && i< 8) seed_index = i;
else printf("randp_setstream index out of range");
}
long randps(int i) { return rs[i] = rs[i] * 16807 % 2147483647;} //from specific stream
long randp() { return rs[seed_index] = rs[seed_index] * 16807 % 2147483647;}

//sample from uniform distribution
double drand ( double low, double high ) {
return ( (double)randp() * ( high - low ) ) / (double)RAND_MAX_P + low;
}
//random integer
int irand ( int low, int high ) {
return ( randp() % ( high - low ) ) + low;
}

double choose2(double c1,double c2) { return (randp() % 2 ==0 ? c1 : c2); }

// adds src track with dst track
void add(double *src, double *dst, int n) {
for (int i=0;i<n;i++) { dst[i] += src[i]; }
}

//modulo operator for floating point
double dmod(double x, double y) {
return x - (int)(x/y) * y;
}

// multiplies src track with dst track (ring modulator)
void mult(double *src, double *dst, int n) {
for (int i=0;i<n;i++) { dst[i] *= src[i]; }
}

// scale and offset
void scale(double a,double b, double *dst, int n) {
for (int i=0;i<n;i++) { dst[i] = dst[i]*a + b; }
}

void fill(double a, double *env,int n) {
for (int i=0;i<n;i++) {
env[i] = a;
}
}

// shift back in time by 'b', and zero fill the rest
// to account for group delay in filters
void shift(int b, double *data,int n) {
for (int i=0;i<n;i++) data[i] = data[i+b];
for (int i=n-b;i<n;i++) data[i] = 0;
}
//delay
void delay(int b, double *data,int n) {
for (int i=n-1;i>=b;i--) data[i] = data[i-b];
for (int i=0;i<b;i++) { data[i] = 0; }
}


double *make_buf(double dur,int*n) {
int num_samples=dur * SAMPLE_RATE;
if (n!=NULL) { *n = num_samples; }
return (double *)calloc(num_samples, sizeof(double));
}

double *make_bufn(int n) {
return (double *)calloc(n, sizeof(double));
}
double *dup_buf(double *src,int n) {
double *dst= (double *)malloc(n * sizeof(double));
for (int i=0;i<n;i++)
dst[i] = src[i];
return dst;
}
// crossfade src track with dst track
void lerp(double a, double *src, double *dst, int n) {
for (int i=0;i<n;i++) { dst[i] = a*src[i] + (1.0-a)*dst[i]; }
}
// normalize, optionally pass maximum to use
void normalize(double *data,int n) {
double m = maxsample(data,n);
if (m < 1e-009) { printf("Normalize: maxsample < 1e-009\n"); return; }
gain(1.0/m,data,n);
}
double maxsample(double *data,int n) {
double m = 0;
for (int i=0;i<n;i++) {
double x = fabs(data[i]);
if (x > m) m = x;
}
return m;
}

void softclip(double gain,double offset,double *data,int n) {
for (int i=0;i<n;i++) {
double x = data[i] * gain + offset;
if (x < -1.0) {
data[i] = -(2.0/3.0);
} else if (x > 1.0) {
data[i] = 2.0/3.0;
} else {
data[i] = x - (x*x*x)/3.0;
}
data[i] /= gain;
data[i] *= 3.0/2.0;
}
}

void stereomixdown(double *lchan,double*rchan,int n) {

//optional delay
//fdelay(0.8,-10,300,1000,lchan,n);
//fdelay(0.4,-10,300,2000,rchan,n);

//normalize channels

double m = maxsample(lchan,n);
double m2 = maxsample(rchan,n);
if (m < m2) m = m2;
gain(.98/m,lchan,n);
gain(.98/m,rchan,n);

FILE *f = fopen("output.wav","wb");
write_wav_header(f,n);
write_wav_data(f,lchan,rchan,n);
fclose(f);
}

void toBigEndian(int *a) {
int b=0x0000;
for (int i=0;i<4;i++) {
b <<= 8;
b |= (*a & 0xff);
*a >>= 8;
}
*a = b;
}

typedef struct typeWavHeader{
//Riff header
int ChunkID;
int ChunkSize;
int Format;

// WAVE format
int Subchunk1ID;
int Subchunk1Size;
short AudioFormat;
short NumChannels;
int SampleRate;
int ByteRate;
short BlockAlign;
short BitsPerSample;

//Data
int Subchunk2ID;
int Subchunk2Size;
} WavHeader;

void write_wav_header(FILE *f,int numSamples) {
WavHeader h;

h.ChunkID = 0x52494646;
h.ChunkSize = 0;
h.Format = 0x57415645;
h.Subchunk1ID = 0x666d7420;
h.Subchunk1Size = 16;
h.AudioFormat = 1;
h.NumChannels = 2; //Stereo=2,Mono=1
h.SampleRate = 44100;
h.BitsPerSample = 16; //16 bit
h.ByteRate = h.SampleRate * h.NumChannels * h.BitsPerSample/8;
h.BlockAlign = h.NumChannels * h.BitsPerSample/8;

h.Subchunk2ID = 0x64617461;
h.Subchunk2Size = numSamples * h.NumChannels * h.BitsPerSample/8;

h.ChunkSize = 36+h.Subchunk2Size; //now calculate this

//switch big endian
toBigEndian(&h.ChunkID);
toBigEndian(&h.Format);
toBigEndian(&h.Subchunk1ID);
toBigEndian(&h.Subchunk2ID);

fwrite(&h,sizeof(h),1,f);
}

void write_wav_data(FILE *f,double *lchannel,double *rchannel,int num_samples) {
short d;
for (int i =0;i<num_samples;i++) {
d = lchannel[i] * 32767;
fwrite(&d,sizeof(short),1,f);
d = rchannel[i] * 32767;
fwrite(&d,sizeof(short),1,f);
}
}

//1st order IIR lowpass at cutoff freq fc (-3db point)
void lowpass(double fc,double *data,int n) {
double pole = 1 - 4*fc/SAMPLE_RATE; //location of pole in transfer function H(z)
double GAIN = 2.0/(1-pole);

double x[2] = {0,0}; double y[2] = {0,0};
for (int i=0;i<n;i++) {
x[0] = x[1];
x[1] = data[i]/GAIN;
y[0] = y[1];
y[1] = (x[0] + x[1]) + (pole * y[0]);
data[i] = y[1];
}
}

//1st order IIR highpass at cutoff freq fc (-3db point)
void highpass(double fc,double *data,int n) {
double pole = 1 - 4*fc/SAMPLE_RATE; //location of pole in transfer function H(z)
double GAIN = 2.0/(1+pole);

double x[2] = {0,0}; double y[2] = {0,0};
for (int i=0;i<n;i++) {
x[0] = x[1];
x[1] = data[i]/GAIN;
y[0] = y[1];
y[1] = (x[0] - x[1]) + (pole * y[0]);
data[i] = y[1];
}
}
void bandpass(double fl,double fh,double *data,int n) {
highpass(fl,data,n);
lowpass(fh,data,n);
}

void normalize_stereo(double *lchan, double *rchan,int n) {
double m = maxsample(lchan,n);
double m2 = maxsample(rchan,n);
if (m < m2) m = m2;
if (m < 1e-009) { printf("Normalize: maxsample < 1e-009\n"); return; }
gain(1.0/m,lchan,n);
gain(1.0/m,rchan,n);
}

//white noise
void whitenoise(double gain, double *data,int n) {
for (int i=0;i<n;i++) {
//non deterministic drand
//TODO: take this from separate PRNG stream
data[i] = drand(-gain,gain);
}
}

void msmix(double *src,double dbgain, double pos, double *dstl, double *dstr,int n) {
pos = limit(-1.0,1.0,pos);
double g = exp(dbgain/20 * log(10));
for (int i=0;i<n;i++) {
//constant power (-3db) pan law
dstl[i] += g*cos((pos+1.0)*0.25*M_PI)*src[i];
dstr[i] += g*sin((pos+1.0)*0.25*M_PI)*src[i];
}
}

void tomono(double *lchan, double *rchan, double *mono, int n) {
lerp(1.0,lchan,mono,n);
add(rchan,mono,n);
}

//clamp high and low values
double limit(double low,double high,double x) {
return fmin(fmax(x,low),high);
}
//digital resonator described by Platt in speech synthesis
void res2a(double fc,double bw, double *data,int n) {
double C = -1.0 * exp(-2.0*M_PI*bw/SAMPLE_RATE);
double B = 2.0 * exp(-1.0*M_PI*bw/SAMPLE_RATE) * cos(2*M_PI*fc/SAMPLE_RATE);
double A = 1.0 - B - C;

double x[3] = {0,0,0}; double y[3] = {0,0,0};

for (int i=0;i<n;i++) {
x[0] = x[1]; x[1] = x[2]; x[2] = data[i];
y[0] = y[1]; y[1] = y[2];
y[2] = A * x[2] + B * y[1] + C *y[0];
data[i] = y[2];
}
}


// ENVELOPES

//exponential envelope
void expenv(double tc, double *env, int env_n) {
double lambda = 1.0/tc;
for (int i=0;i<env_n;i++) {
double t = (double)i/SAMPLE_RATE;
env[i] = exp(-lambda*t);
}
}

//linear envelope from a to b
void line(double a,double b,double dur,double *env,int env_n) {
for (int i=0;i<dur * SAMPLE_RATE;i++) {
if (i >= env_n) continue;
env[i] = (b-a)/(dur*SAMPLE_RATE) * i + a;
}
}

//naive triangle tone, suitable for lfo
void sin_lfo(double freq,double amp,double offset,double *data,int n) {
for (int i=0;i<n;i++) {
data[i] = offset + amp * sinlu(freq * 2 * M_PI * (double)i/SAMPLE_RATE);
}
}

//exponential adsr envelope
void expadsr(double a,double d, double s, double r, double dur, double *env, int env_n) {
//TODO tweak time constants for different parts of the envelope
int i=0;
double t =0;
double lambda = 1.0/(a/6);
if (a < 1e-9) lambda = 1e9;
int ac = 0;

while (i<a*SAMPLE_RATE && (ac+i) < env_n) {
t = (double)i/SAMPLE_RATE;
env[ac+i] = exp(-lambda*t) * (0-1.0) + 1.0;
i++;
}
// printf("%.3f\n",(double)(ac+i)/SAMPLE_RATE);
ac += i; i=0;
lambda = 1.0/ (d/8);
if (d < 1e-9) lambda = 1e9;

while (i<(d+dur)*SAMPLE_RATE && (ac+i) < env_n) {

t = (double)i/SAMPLE_RATE;
//TODO change time constant to affect the decay val
env[ac+i] = exp(-lambda*0.5*t) * (1.0 - s) + s;
i++;
}
ac += i; i=0;
lambda = 1.0/ (r/6);
if (r < 1e-9) lambda = 1e9;
while (i<r*SAMPLE_RATE && (ac+i) < env_n) {

t = (double)i/SAMPLE_RATE;
env[ac+i] = exp(-lambda*t) * s;
i++;
}

//zero fill rest of envelope
while (ac+i < env_n) {
env[ac+i++] = 0;
}

}

typedef struct {
double *Xr2;
double *Xi2;
double *cos2pidN;
double *sin2pidN;
int *ind;
int N;
} fft_plan;

void prepare_fft(fft_plan *f,int N) {
if (f->N == N) {
return;
}
if (f->N != 0) {
if (f->Xr2) free(f->Xr2);
if (f->Xi2) free(f->Xi2);

if (f->cos2pidN) free(f->cos2pidN);
if (f->sin2pidN) free(f->sin2pidN);
if (f->ind) free(f->ind);
}
printf("Preparing fft %d\n",N);
f->N = N;
f->Xr2 = (double *)malloc(N * sizeof(double));
f->Xi2 = (double *)malloc(N * sizeof(double));

f->cos2pidN = (double *)malloc(N * sizeof(double));
f->sin2pidN = (double *)malloc(N * sizeof(double));
f->ind = (int *)malloc(N * sizeof(int));

//precompute cos/sin table
for (int i=0;i<N;i++) {
f->cos2pidN[i] = cos(2*M_PI*(double)i/N);
f->sin2pidN[i] = -sin(2*M_PI*(double)i/N);
}

//build index
f->ind[0] = 0;
for (int n= N; n>1; n /=2) {
for (int i=0;i<N;i+=n) {
f->ind[i+n/2] = f->ind[i] + N/n;
}
}
}

void fft(double *x, double *Xr, double *Xi,fft_plan *f,int N) {
double *X = Xr;
if (N != f->N) {
prepare_fft(f,N);
}

double *Xr2 = f->Xr2;
double *Xi2 = f->Xi2;

//compute first layer
for (int i=0;i<N;i+=2) {
Xr[i] = x[f->ind[i]] + x[f->ind[i+1]];
Xr[i+1] = x[f->ind[i]] - x[f->ind[i+1]];

Xi[i] = 0;
Xi[i+1] = 0;
}

double *Er,*Or,*Ei,*Oi;

for (int n = 2;n<N;n*=2) {
int CI = (N/(n*2));
for (int j = 0;j<N;j+=n*2) {
Er = &Xr[j];
Or = &Xr[j+n];

Ei = &Xi[j];
Oi = &Xi[j+n];

double *Xr2ij=&Xr2[j];
double *Xi2ij=&Xi2[j];
double *Xr2ijn=&Xr2[j+n];
double *Xi2ijn=&Xi2[j+n];

double *sint = f->sin2pidN;
double *cost = f->cos2pidN;
for (int i=0;i<n;i++,sint+=CI,cost+=CI) {
//compute n*2 point DFT
//1st half

*Xr2ij++ = *Er + (*Or*(*cost) - *Oi*(*sint));
*Xi2ij++ = *Ei + (*Or*(*sint) + *Oi*(*cost));

//2nd half (TODO: omit in last level)

*Xr2ijn++ = *Er++ - (*Or*(*cost) - *Oi*(*sint));
*Xi2ijn++ = *Ei++ - (*Or++ *(*sint) + *Oi++ *(*cost));
}
}

//now swap pointers (TODO: omit in last level)
double *tmp = Xr;Xr = Xr2;Xr2 = tmp;
tmp = Xi;Xi = Xi2; Xi2 = tmp;
}

// copy to passed array if required
if (X != Xr) {
double *tmp = Xr;Xr = Xr2;Xr2 = tmp;
tmp = Xi;Xi = Xi2; Xi2 = tmp;
//copy Xr2
for (int j=0;j<N;j++) { Xr[j] = Xr2[j]; Xi[j] = Xi2[j];
}
}
}

void ifft(double *x, double *Xr, double *Xi,fft_plan *f,int N) {
double *X = Xr;
if (N != f->N) {
prepare_fft(f,N);
}
double *Xr2 = f->Xr2;
double *Xi2 = f->Xi2;

//compute first layer
for (int i=0;i<N;i+=2) {
Xr2[i] = Xr[f->ind[i]] + Xr[f->ind[i+1]];
Xr2[i+1] = Xr[f->ind[i]] - Xr[f->ind[i+1]];

Xi2[i] = Xi[f->ind[i]] + Xi[f->ind[i+1]];
Xi2[i+1] = Xi[f->ind[i]] - Xi[f->ind[i+1]];
}

//now swap pointers (TODO: omit in last level)
double *tmp = Xr;Xr = Xr2;Xr2 = tmp;
tmp = Xi;Xi = Xi2; Xi2 = tmp;

double *Er,*Or,*Ei,*Oi;

for (int n = 2;n<N;n*=2) {
int CI = (N/(n*2));
for (int j = 0;j<N;j+=n*2) {
Er = &Xr[j];
Or = &Xr[j+n];

Ei = &Xi[j];
Oi = &Xi[j+n];

double *Xr2ij=&Xr2[j];
double *Xi2ij=&Xi2[j];
double *Xr2ijn=&Xr2[j+n];
double *Xi2ijn=&Xi2[j+n];

double *sint = f->sin2pidN;
double *cost = f->cos2pidN;
for (int i=0;i<n;i++,sint+=CI,cost+=CI) {
//compute n*2 point DFT
//1st half

*Xr2ij++ = *Er + (*Or*(*cost) + *Oi*(*sint));
*Xi2ij++ = *Ei + (*Or*(-*sint) + *Oi*(*cost));

//2nd half (TODO: omit in last level)

*Xr2ijn++ = *Er++ - (*Or*(*cost) + *Oi*(*sint));
*Xi2ijn++ = *Ei++ - (*Or++ *(-*sint) + *Oi++ *(*cost));
}
}

if (n == N/2 && X == Xr) {
//don't swap
}else {
//now swap pointers (TODO: omit in last level)
double *tmp = Xr;Xr = Xr2;Xr2 = tmp;
tmp = Xi;Xi = Xi2; Xi2 = tmp;
}
}

// copy to passed array if required
if (X != Xr) {
double *tmp = Xr;Xr = Xr2;Xr2 = tmp;
tmp = Xi;Xi = Xi2; Xi2 = tmp;
//copy Xr2
for (int j=0;j<N;j++) { Xr[j] = Xr2[j]; Xi[j] = Xi2[j]; }
}

for (int j=0;j<N;j++) { x[j] = Xr[j]/(double)N;}
}

//convolution using fft
void convfft(double *src,int src_n,double *dst,int dst_n) {
int N = 1;
while (src_n+dst_n > N) N*=2;
N*=2;
fft_plan f;
f.N = 0;

printf("Performing FFT with N = %d\n",N);
printf("Src: %d, DST: %d\n",src_n,dst_n);

double *srcx = (double*)malloc(N * sizeof(double));
double *SRCr = (double*)malloc(N * sizeof(double));
double *SRCi = (double*)malloc(N * sizeof(double));

fill(0,srcx,N);
add(src,srcx,src_n);

fft(srcx,SRCr,SRCi,&f,N);
printf("Completed FFT\n");
double *dstx = (double*)malloc(N * sizeof(double));
double *DSTr = (double*)malloc(N * sizeof(double));
double *DSTi = (double*)malloc(N * sizeof(double));

fill(0,dstx,N);
add(dst,dstx,dst_n);

fft(dstx,DSTr,DSTi,&f,N);

for (int i=0;i<N;i++) {
double Xr = DSTr[i];
double Xi = DSTi[i];

DSTr[i] = Xr * SRCr[i] - Xi*SRCi[i];
DSTi[i] = Xr * SRCi[i] + Xi*SRCr[i];
}

ifft(dstx,DSTr,DSTi,&f,N);
for (int i=0;i<dst_n;i++) {
dst[i] = dstx[i];
}

free(srcx); free(SRCr); free(SRCi);
free(dstx); free(DSTr); free(DSTi);
}

//filtered white noise exponential trail
//g gain in db
void reverb(double dur, double g, double predelay, double *src,double *dst,int n) {
int dn;
double *buf = make_buf(dur*1.5 + predelay,&dn);
double *env = make_bufn(dn);

expenv(dur/6,env,dn);
if (predelay > 1e-6) delay((int)(predelay*SAMPLE_RATE),env,dn);

whitenoise(1.0,buf,dn);
lowpass(1500,buf,dn);
highpass(100,buf,dn);
highpass(100,buf,dn);
mult(env,buf,dn);

lerp(1.0,src,dst,n);
convfft(buf,dn,dst,n);
dbgain(g,dst,n);

free(buf);
free(env);
}

// BSETS

void bset_init(bset *b) {
b->interp = INTERP_LINEAR; //linear by default
b->n = 0;
b->max_n = 2048*4;
b->sr = SAMPLE_RATE;
b->t = (double *)calloc(b->max_n,sizeof(double));
b->y = (double *)calloc(b->max_n,sizeof(double));
b->iparam = (double *)calloc(b->max_n,sizeof(double)); //param (gamma, tau);
b->itype = (int *)calloc(b->max_n,sizeof(int)); //interpolation type
}

double bset_length(bset *b) {if (b->n==0) return 0;return b->t[b->n-1];
}

double bset_last(bset *b) {if (b->n==0) return 0;return b->y[b->n-1];}

void bset_destroy(bset *b) {
if (b->t) free(b->t);
if (b->y) free(b->y);
if (b->iparam) free(b->iparam);
if (b->itype) free(b->itype);
b->t =0;
b->y =0;
}


//returns index of left bound of interval containing t
int bset_search(bset *b, double t) {
//TODO: optimization, change this to binary search
int i=0;

for (i=0;i<b->n;i++) {
if (t < b->t[i]) break;
}

return i;
}

//interpolate breakpoint set at time t, linear interpolation
double bset_interp_linear(bset *b, double t) {
int i = bset_search(b,t) - 1;
if (i >= 0 && i < b->n-1) {
//linear interpolation
double a = (t - b->t[i])/ (b->t[i+1]-b->t[i]);
return (1-a)*b->y[i] + a*b->y[i+1];
}
return 0;
}
//interpolate breakpoint set at time t, exponential interpolation
double bset_interp_exp(bset *b, double t, double tau) {
int i = bset_search(b,t) - 1;

if (i >= 0 && i < b->n-1) {
//exponential interpolation
double a = (t - b->t[i])/ (b->t[i+1]-b->t[i]);
double trans_val = 0.75;
if (a >= trans_val) {
double exp_v = (b->y[i] - b->y[i+1])*exp(-1.0*a*tau) + b->y[i+1];
double lin1 = (b->y[i] - b->y[i+1])*exp(-1.0*trans_val*tau) + b->y[i+1];
double c = (a-trans_val)/(1-trans_val);
double lin_v = (1-c)*lin1 + c*b->y[i+1];
return (1-c)*exp_v + c*lin_v;
} else {
return (b->y[i] - b->y[i+1])*exp(-1.0*a*tau) + b->y[i+1];
}
}
return 0;
}

//interpolate breakpoint set at time t, power law
double bset_interp_gamma(bset *b, double t,double gamma) {
int i = bset_search(b,t) - 1;

if (i >= 0 && i < b->n-1) {
//power law interpolation
double a = (t - b->t[i])/ (b->t[i+1]-b->t[i]);
a = pow(a,1.0/gamma);
return (1-a)*b->y[i] + a*b->y[i+1];
}
return 0;
}

//interpolate breakpoint set at time t, sample/hold
double bset_interp_cubic(bset *b, double t) {
int i = bset_search(b,t) - 1;
if (i >= 0 && i < b->n-1) {
//cubic interpolation
double a = -2*(b->y[i+1] - b->y[i]);
double c = 3 * (b->y[i+1] - b->y[i]);

double d = (t - b->t[i])/(b->t[i+1] - b->t[i]);
return a*d*d*d + c*d*d + b->y[i];
}
return 0;
}

//interpolate breakpoint set at time t, sample/hold
double bset_interp_samplehold(bset *b, double t) {
int i = bset_search(b,t) - 1;
if (i >= 0 && i < b->n-1) {
//sample and hold
return b->y[i];
}
return 0;
}

//interpolate breakpoint set, using the indicated interpolation method for that segment
double bset_interp(bset *b, double t) {
int i = bset_search(b,t) - 1;
if (i >= 0 && i < b->n-1) {
switch(b->itype[i]) {
case INTERP_SAMPLEHOLD:
return bset_interp_samplehold(b,t);
break;
case INTERP_LINEAR:
return bset_interp_linear(b,t);
break;
case INTERP_CUBIC:
return bset_interp_cubic(b,t);
break;
case INTERP_GAMMA:
return bset_interp_gamma(b,t,b->iparam[i]); //default gamma param 5.0
break;
case INTERP_EXP:
return bset_interp_exp(b,t,b->iparam[i]); //default exp param 6.0
break;
default:
return bset_interp_samplehold(b,t); // sample hold by default
break;
}
}
}

//render breakpoint set to an array of doubles
double bset_render(bset *b, double *data) {
double t = 0;
int k=0;
while (t < b->t[b->n-1]) {
data[k++] = bset_interp(b,t);
t+=1.0/b->sr;
}
}

void bset_add(bset *b, double t, double y) {
if (b->n >= b->max_n) { printf("more than max_n %d entries breakpoint set",b->max_n); return; }

//insertion sort
int i = bset_search(b,t);

for (int j=b->n;j>i;j--) {
b->t[j] = b->t[j-1]; b->y[j] = b->y[j-1];
b->itype[j] = b->itype[j-1];
b->iparam[j] = b->iparam[j-1];
}
b->t[i] = t; b->y[i] = y;
b->itype[i] = b->interp; //default interpolation method for segment
b->iparam[i] = b->param; //default interpolation param
b->n++;
}

//overwrite with linear segment
void bset_segment(bset *b,double t1, double t2, double y1, double y2) {
//first remove any in the way
for (int i=b->n-1;i>=0;i--) {
if (b->t[i] > t1 && b->t[i] <= t2) {
bset_remove(b,i);
}
}
bset_add(b,t1,y1);
bset_add(b,t2,y2);
}
void bset_remove(bset *b, int j) {
if (j >= b->n) return;

for (int i = j;i<b->n-1;i++) {
b->t[i] = b->t[i+1]; b->y[i] = b->y[i+1];
b->itype[i] = b->itype[i+1];
b->iparam[i] = b->iparam[i+1];
}
b->n--;
}

void bset_clear(bset *b) { for (int i=0;i<b->n;i++) { b->t[i] = 0; b->y[i] = 0; b->itype[i] = 0; b->iparam[i] = 0;} b->n =0; }
void bset_sort(bset *b) {
bset b2;
bset_init(&b2);

for (int i=0;i<b->n;i++) {
bset_add(&b2,b->t[i],b->y[i]);
}
bset_clear(b);
for (int i=0;i<b2.n;i++) {
bset_add(b,b2.t[i],b2.y[i]);
}
}
void bset_stretch(bset *b,double f) {
for (int i=0;i<b->n;i++) {
b->t[i] *= f;
}
}

//crops the segment beginning at t1 with duration dur
void bset_crop(bset *b,double t1,double dur) {
//traverse backwards so points can be safely removed
for (int i=b->n-1;i>=0;i--) {
b->t[i] -= t1;
if (b->t[i] < 0 || b->t[i] > dur) bset_remove(b,i);
}
}

void bset_shift(bset *b,double f) {
//traverse backwards so points can be safely removed
for (int i=b->n-1;i>=0;i--) {
b->t[i] -= f;
if (b->t[i] < 0) bset_remove(b,i);
}
}

void bset_delay(bset *b,double f) {
for (int i=b->n-1;i>=0;i--) {
b->t[i] += f;
}
}

double linear_interp_circular(double p, double *data,int n) {
if (p < 0) return 0;

double a1 = p - floor(p);
int i1 = (int)p;
int i2 = i1 + 1;
i1 = i1 % n;
i2 = i2 % n;

return a1*data[i2] + (1.0-a1)*data[i1];

}
double linear_interp(double p, double *data,int n) {
if (p < 0) return 0;

double a1 = p - floor(p);
int i1 = (int)p;
int i2 = i1 + 1;
if (i2 >= n) return 0;

return a1*data[i2] + (1.0-a1)*data[i1];
}


//update filter coefficients based on change in center frequency, bandwidth etc
void iir_update_params(iirfilter_params *f) {
if (f->ftype==FTYPE_LP) {
f->pole = 1 - 4*f->fc/SAMPLE_RATE;
f->gain = 2.0/(1-f->pole);
} else if (f->ftype == FTYPE_HP) {
f->pole = 1 - 4*f->fc/SAMPLE_RATE;
f->gain = 2.0/(1+f->pole);
} else if (f->ftype == FTYPE_RES2) {
double R = exp(-M_PI*f->bw/SAMPLE_RATE); //pole radius in terms of bandwidth
double A = 2*M_PI*f->fc/SAMPLE_RATE; // pole angle

f->a1 = -2*R*cos(A);
f->a2 = R*R;
f->gain = 100 * (1-R*R)/2.0;
} else if (f->ftype == FTYPE_RES2A) {
f->rC = -1.0 * exp(-2.0*M_PI*f->bw/SAMPLE_RATE);
f->rB = 2.0 * exp(-1.0*M_PI*f->bw/SAMPLE_RATE) * cos(2*M_PI*f->fc/SAMPLE_RATE);
f->rA = 1.0 - f->rB - f->rC;
}
}

void init_iirfilter(iirfilter_params *f) {
for (int i=0;i<3;i++) { f->x[i] = 0; f->y[i] = 0; }
iir_update_params(f);
}

//one step of iir filter
double iir_apply(iirfilter_params *f, double s) {
if (f->ftype==FTYPE_LP) {
f->x[0] = f->x[1];
f->x[1] = s/f->gain;

f->y[0] = f->y[1];
f->y[1] = (f->x[0] + f->x[1]) + (f->pole * f->y[0]);
return f->y[1];
} else if (f->ftype==FTYPE_HP) {
f->x[0] = f->x[1];
f->x[1] = s/f->gain;

f->y[0] = f->y[1];
f->y[1] = (f->x[0] - f->x[1]) + (f->pole * f->y[0]);
return f->y[1];
} else if (f->ftype==FTYPE_RES2) {
double *x = f->x; double *y= f->y;
x[0] = x[1];
x[1] = x[2];
x[2] = s * f->gain;

y[0] = y[1];
y[1] = y[2];
y[2] = x[2] - x[0] - f->a1* y[1] -f->a2* y[0];
return y[2];
} else if (f->ftype==FTYPE_RES2A) {
double *x = f->x; double *y= f->y;
x[0] = x[1]; x[1] = x[2]; x[2] = s;
y[0] = y[1]; y[1] = y[2];
y[2] = f->rA * x[2] + f->rB * y[1] + f->rC *y[0];
return y[2];
} else {
//unknown filter, just pass to output
return s;
}

return s;
}

//simple delay using delay buffer, and optional hp and lp filters
void ffdelay(double dur, double feedback, double fl, double fh, double *src, double*dst, int n) {
feedback = exp(feedback/20 * log(10));

int d = (dur * SAMPLE_RATE) + 1;
double *buf = (double *)malloc(d * sizeof(double));
for (int i=0;i<d;i++) buf[i] = 0;

iirfilter_params filt_hp;
filt_hp.ftype = FTYPE_HP;
filt_hp.fc = fl;

iirfilter_params filt_lp;
filt_lp.ftype = FTYPE_LP;
filt_lp.fc = fh;

init_iirfilter(&filt_hp);
init_iirfilter(&filt_lp);

int p = 0;
for (int i=0;i<n;i++) {
//double x = buf[p];

//apply any iir filters to buf[p] here
if (fl > .001) buf[p] = iir_apply(&filt_hp, buf[p]);
if (fh > .001) buf[p] = iir_apply(&filt_lp, buf[p]);
buf[p] = src[i] + feedback*buf[p];
//data[i] += gain * x;
dst[i] += buf[p];

if (++p >= d) p = 0;
}
free(buf);
}

//envelope follower
//attack and release in milliseconds
void envdetect(double attack, double release, double *src,double *dst,int n) {
// double a = exp(-1000*attack*3/SAMPLE_RATE);
// double r = exp(-1000*release*3/SAMPLE_RATE);
double a = pow( 0.01, 1.0 / ( attack * SAMPLE_RATE * 0.001 ) );
double r = pow( 0.01, 1.0 / ( release * SAMPLE_RATE * 0.001 ) );
// printf("attack: %.9f, release: %.9f\n",a,r);

// double a = pow(100,-1000*attack/SAMPLE_RATE);
// double r = pow(100,-1000*release/SAMPLE_RATE);

double x = 0;
for (int i=0;i<n;i++) {
double b = fabs(src[i]);
if (x < b) {
x = a * (x-b) + b;
} else {
x = r * (x-b) + b;
}
dst[i] = x;
}
}

// attenuate signal based on control signal, hard knee threshold/ratio
void compress(double ratio, double threshold, double *cs, double *src,double *dst,int n) {
for (int i=0;i<n;i++) {
if (cs[i] <= threshold) {
dst[i] = src[i];
} else {

double over_db = 20*log(cs[i]/threshold); //db over threshold
double atten = -over_db * (1 - 1.0/ratio); //attenuation in db
// printf("atten: %.3f\n",atten);
dst[i] = src[i] * exp(atten/20* log(10));
}
}
}


void compress_inplace(double attack, double release, double db, double *src,int n) {
double *env = make_bufn(n);
double *dst = make_bufn(n);
envdetect(attack,release,src,env,n);

double thresh = maxsample(env,n) * exp(-1.0*db/20 * log(10));

// printf("top sample: %.3f\n",maxsample(src,n));
// printf("thresh: %.3f\n",thresh);
double ratio = 4;
compress(ratio,thresh,env,src,dst,n);
lerp(1.0,dst,src,n);
free(env);
free(dst);
}

void compress_inplace_stereo(double attack, double release, double db, double *src0, double *src1, int n) {
double *env = make_bufn(n);
double *mono = make_bufn(n);
double *dst0 = make_bufn(n);
double *dst1 = make_bufn(n);

tomono(src0,src1,mono,n);
envdetect(attack,release,mono,env,n);

double thresh = maxsample(env,n) * exp(-1.0*db/20 * log(10));

double ratio = 2;
compress(ratio,thresh,env,src0,dst0,n);
compress(ratio,thresh,env,src1,dst1,n);
lerp(1.0,dst0,src0,n);
lerp(1.0,dst1,src1,n);

free(env);
free(mono);
free(dst0);free(dst1);
}


//generators

void init_gen(gen_data *g) {
g->freq = 440;

g->phase = 0;

g->pw = 0.50;
g->amp = 1.0;
g->i = 0;

//for DSF
g->freq2 = 440;
g->phase2 = 0;
g->w = 0.5;
}

void sin_tone(int freq,double start,double dur,double amp,double *data) {
int s;
for (int i=0;i<dur*SAMPLE_RATE;i++) {
data[i + (int)(start*SAMPLE_RATE)] = amp * sinlu(freq * 2 * M_PI * (double)i/SAMPLE_RATE);
}
}

//sin tone with time varying frequency envelope
void sin_tone_t(double *fenv, double env_n, double start, double dur,double amp, double *data,int n) {
gen_data g;

init_gen(&g);
g.amp = amp;
g.freq = fenv[0];


int offset = start*SAMPLE_RATE;
int ie = 0;
for (int i=0;i<dur*SAMPLE_RATE;i++,ie++) {
if (ie >= env_n) ie = env_n - 1;
g.freq = fenv[ie];
data[offset + i] = sin_gen(&g); //TODO: change this to a function pointer
}
}

double sin_gen(gen_data *g) {
g->phase += 2*M_PI*g->freq/SAMPLE_RATE;
if (g->phase > 2*M_PI) g->phase -= 2*M_PI;
g->phase = dmod(g->phase,2*M_PI);
double x = g->amp * sinlu(g->phase);
return x;
}

// sin lookup table
//int sinl_table_n=2048*256; //4mb lookup table gives about 70db SNR to math sin()
int sinl_table_n=2048*4; //4mb lookup table gives about 70db SNR to math sin()
double *sinl_table=0;
double sinlu(double w) {
if (sinl_table==0) {
printf("Building %d sin lookup table.. ",sinl_table_n);
sinl_table = (double*)malloc(sinl_table_n * sizeof(double));
for (int i=0;i<sinl_table_n;i++) {
sinl_table[i] = sin((double)i*2.0*M_PI/(double)sinl_table_n);
}
printf("done.\n");
}

w = dmod(w, 2.0*M_PI);
if (w < 0) w += 2.0*M_PI;

double p = ((double)w*sinl_table_n/(2.0*M_PI));
//return sinl_table[(int)round(p)];
//optional linear interpolation of sin table

double a1 = p - floor(p);
int i1 = (int)p;
int i2 = i1 + 1;
if (i2 >= sinl_table_n) i2 = 0;

return a1*sinl_table[i2] + (1.0-a1)*sinl_table[i1];
}




//trapezoidal window
void trapz(int a, double *env, int n) {
fill(1.0,env,n);
for (int i=0;i<a;i++) {
env[i] = (double)i/a;
env[n-1-i] = env[i];
}
}

//FM modulation
void modulate_t2(double *A, double *freq, double *src, int src_n, double *dst) {
gen_data g;
init_gen(&g);
g.freq = freq[0] + A[0]*src[0];
int ie = 0;
for (int i=0;i<src_n;i++,ie++) {
g.freq = freq[i] + A[i]*src[i];
dst[i] = sin_gen(&g);
}
}

void fadeout(double dur, double *src,int n) {
int d = dur*SAMPLE_RATE;
if (d > n) d = n;
double t = 0;
for (int i=n-d;i<n;i++) {
double a = (1-t/dur);
src[i] *= a;
t += 1.0/SAMPLE_RATE;
}
}

//pulse envelope made to the rhythm of a beat
//pw in % of dur
void beat_env(double dur, double pw, double *pat, int patn, int n, double *env,int env_n) {
bset b;
bset_init(&b);

double t = 0;
double delta = .001;
for (int i=0;i<n;i++) {
bset_add(&b,t,pat[i % patn]);
bset_add(&b,t+pw*dur-delta,pat[i % patn]);
// bset_add(&b,t+pw*dur-delta,0);
// bset_add(&b,t+delta,pat[i%patn]/2);

t += dur;
}
b.interp = INTERP_CUBIC;
bset_render(&b,env);
lowpass(200,env,env_n);
}

//pulse envelope made to the rhythm of a beat
//pw in % of dur
void beat_env2(double dur,
double *pat, int patn,
double *attack, int attackn,
double *pw, int pwn,
double *release, int releasen,
int n, double *env,int env_n) {
bset b;
bset_init(&b);
int ns = dur*SAMPLE_RATE;
double *e = make_bufn(ns);

int pt = 0;
double delta = .001;
for (int i=0;i<n;i++) {
double amp = pat[i%patn];
double a = attack[i%attackn];
double d = dur * pw[i%pwn];
double r = release[i%releasen];

// if (d-a-r < 0) { r = d - a; }

expadsr(a, //attack
0, //decay
1.0, //sustain
r, //release
d-a-r,// dur
e,
ns);

scale(amp,0,e,ns);
lerp(1.0,e,&env[pt],ns);

pt += ns;
}
free(e);
lowpass(200,env,env_n);
}

//visual code

typedef unsigned int int32;
typedef short int16;
typedef unsigned char byte;

typedef struct img_t {
int h;
int w;
int n; //number of values. 1 for grayscale, 3 for rgb
double *p; //pixel data
} img_t;


typedef struct {
int use_lines; //for '/-|\' type lines
double c; //otherwise the character to use;
} drawlinesettings_t;

int ind(int x, int y,int w);


void val_line(int x0,int y0, int x1,int y1, double v, double *val,int w, int h);
void val_to_img_ds(int w,int h, int os, double *val,img_t *im);
void val_to_img(int w,int h, double *val,img_t *im);
void val_boxblur(double rhox,double rhoy, double *val, int w, int h);
void val_lens(double R1, double *val, int w, int h);

void write_frame(int w,int h,double os,double *val,img_t *img,int *fn,double *t);
void writebmp(img_t *i,char *filename);
void init_img(img_t *i,int n, int w, int h);
double getpixel(img_t *i, int x, int y, int n);
double *getpixelp(img_t *i, int x, int y);


//end header

//index into 2d array
int ind(int x, int y,int w) {
return y*w + x;
}

double getpixel(img_t *i, int x, int y, int n) {
int index = x*i->h*i->n + y*i->n;
double *p = &(i->p[index]);
return p[n];
}

double *getpixelp(img_t *i, int x, int y) {
int index = x*i->h*i->n + y*i->n;
double *p = &(i->p[index]);
return p;
}

void setpixel(img_t *i,int x, int y, double r,double g, double b) {
int index = x*i->h*i->n + y*i->n;
double *p = &(i->p[index]);

p[0] = r; p[1] = g; p[2] = b;
}


void val_to_img(int w,int h, double *val,img_t *im) {
for (int x=0;x<w;x++) { for (int y=0;y<h;y++) {
double *v = &val[ind(x,y,w)];
int index = x*im->h*im->n + y*im->n;
double *p = &(im->p[index]);
double f = *v;
if (f < 0) f = 0; if (f > 1.0) f = 1.0;
p[0] = f; p[1] = f; p[2] = f;
// setpixel(im,x,y,*v,*v,*v);
}
}
}

//val to img with downsampling
void val_to_img_ds(int w,int h, int os, double *val,img_t *im) {
double os_v = 1.0/(double)(os*os);
double *val_ds = make_bufn(((double)(w*h)/(os*os)));
int w2 = w/os;
int h2 = h/os;

for (int x=0;x<w;x++) {
for (int y=0;y<h;y++) {
int x2 = x/os;
int y2 = y/os;
val_ds[ind(x2,y2,w2)] += os_v * val[ind(x,y,w)];
}
}
// normalize(val_ds,w2*h2);
val_to_img(w2,h2,val_ds,im);
free(val_ds);
}

void writebmp(img_t *i,char *filename) {
int WIDTH = i->w;
int HEIGHT= i->h;

unsigned int headers[13];
FILE * outfile;
int extrabytes;
int paddedsize;
int x; int y; int n;
int red, green, blue;

extrabytes = 4 - ((WIDTH * 3) % 4);
if (extrabytes == 4)
extrabytes = 0;

paddedsize = ((WIDTH * 3) + extrabytes) * HEIGHT;

headers[0] = paddedsize + 54;
headers[1] = 0;
headers[2] = 54;
headers[3] = 40;
headers[4] = WIDTH;
headers[5] = HEIGHT;
headers[7] = 0;
headers[8] = paddedsize;
headers[9] = 0;
headers[10] = 0;
headers[11] = 0;
headers[12] = 0;

outfile = fopen(filename, "wb");
fprintf(outfile, "BM");

for (n = 0; n <= 5; n++) {
fprintf(outfile, "%c", headers[n] & 0x000000FF);
fprintf(outfile, "%c", (headers[n] & 0x0000FF00) >> 8);
fprintf(outfile, "%c", (headers[n] & 0x00FF0000) >> 16);
fprintf(outfile, "%c", (headers[n] & (unsigned int) 0xFF000000) >> 24);
}
fprintf(outfile, "%c", 1);
fprintf(outfile, "%c", 0);
fprintf(outfile, "%c", 24);
fprintf(outfile, "%c", 0);

for (n = 7; n <= 12; n++) {
fprintf(outfile, "%c", headers[n] & 0x000000FF);
fprintf(outfile, "%c", (headers[n] & 0x0000FF00) >> 8);
fprintf(outfile, "%c", (headers[n] & 0x00FF0000) >> 16);
fprintf(outfile, "%c", (headers[n] & (unsigned int) 0xFF000000) >> 24);
}

int stride = (i->h*i->n);
unsigned char *data = (unsigned char*)malloc(WIDTH*3);

for (y = HEIGHT - 1; y >= 0; y--) {
double *p = getpixelp(i,0,y);

for (x = 0; x <= WIDTH - 1; x++) {
data[x*3] = p[2]*255;
data[x*3+1] = p[1]*255;
data[x*3+2] = p[0]*255;
p+= stride;
}
fwrite(data,1,WIDTH*3,outfile);
if (extrabytes) {
for (n = 1; n <= extrabytes; n++) {
fprintf(outfile, "%c", 0);
}
}
}

fclose(outfile);
return;
}

void write_frame(int w,int h,double os,double *val,img_t *img,int *fn,double *t) {
val_to_img_ds(w,h,os,val,img);

char filename[128];
sprintf(filename,"frames/frame%04d.bmp",*fn+1);
writebmp(img,filename);

if (*fn % 10==0) {
printf("Frame %d\n",*fn+1);
}

*t += 1/24.0;
*fn = *fn + 1;
}

//DRAWING bresenham line drawing
drawlinesettings_t dls;
void init_drawline() {
dls.use_lines = 1;
dls.c = '.';
}

void init_img(img_t *i,int n, int w, int h) {
i->h = h; i->w = w; i->n = n;
i->p = (double *)malloc(n*w*h*sizeof(double));
}

//bresenham
void val_line(int x0,int y0,
int x1,int y1,
double v,
double *val,int w, int h) {

int dx = abs(x1-x0);
int sx = x0 < x1 ? 1 : -1;
int dy = -1.0*abs(y1-y0);
int sy = y0<y1 ? 1 : -1;
int err = dx+dy;
int x = x0;
int y = y0;
int k=0;

while (1) {
if (x >= 0 && x < w && y >= 0 && y < h)
val[ind(x,y,w)] += v;

if (k >= 4096) { //printf("LINE overflow: %d\n",k);
return; }
if (floor(x)==floor(x1) && floor(y)==floor(y1)) break;
int e2 = 2*err;
if (e2 >= dy) {
err += dy;
x += sx;
}
if (e2 <= dx) {
err += dx;
y += sy;
}
}
}

//non space varying box blur
void val_boxblur(double rhox,double rhoy, double *val, int w, int h) {
double a = 0;
double mx = rhox*2 + 1;
double my = rhoy*2 + 1;
double *dst = make_bufn(w*h);

for (int i=0;i<w;i++) {
a = 0;
for (int j=0;j<h;j++) {
int xx = i;
int yy= j-my;

a += val[ind(i,j,w)];
if (!(xx < 0 || xx >= w || yy < 0 || yy >= h)) {
a -= val[ind(xx,yy,w)];
}

yy = j - rhoy;
if (!(xx < 0 || xx >= w || yy < 0 || yy >= h)) {
dst[ind(xx,yy,w)] = a;
}
}
}
lerp(1.0,dst,val,w*h);
fill(0,dst,w*h);


for (int j=0;j<h;j++) {
a = 0;
for (int i=0;i<w;i++) {
int xx = i - mx;
int yy= j;

a += val[ind(i,j,w)];
if (!(xx < 0 || xx >= w || yy < 0 || yy >= h)) {
a -= val[ind(xx,yy,w)];
}

xx = i-rhox;
if (!(xx < 0 || xx >= w || yy < 0 || yy >= h)) {
dst[ind(xx,yy,w)] = a;
}
}
}
lerp(1.0,dst,val,w*h);
fill(0,dst,w*h);
double fact = 1.0/(mx*my);
for (int i=0;i<w*h;i++) {
val[i] *= fact;
}

free(dst);
}

//radial lens distoration
//caches the nonlinear transform to accelerate applying to multiple frames
void val_lens(double R1, double *val, int w, int h) {
int n = w*h;
int use_cache = 1;
static int *mapping;
static int init =0;
double R0 = 0;
// double R1 = -0.15; //-0.15;
double R3 = -R1*1.05; //0.2
double R5 = 0; //R3*1.25;

double *buf = make_bufn(n);

if (use_cache) {
if (init==1) {
//map is already in place
for (int i=0;i<n;i++) {
if (mapping[i] >= 0) {
buf[i] = val[mapping[i]];
} else {
buf[i] = 0.1;
}
}
lerp(1.0,buf,val,n);
free(buf);
return;
}
mapping = (int *)calloc(n,sizeof(int));
init = 1;
}
int xc = w/2.0;
int yc = h/2.0;

fill(0.3,buf,n);

for (int x=0;x<w;x++) { for (int y = 0;y<h;y++) {
// x =w/2; y=0;
double xr = ((double)x - xc)/xc;
double yr = ((double)y - yc)/yc;

double r = sqrt(xr*xr + yr*yr);
// printf("r: %.3f xr: %.3f yr: %.3f\n",r,xr,yr);
double phi = M_PI/2;
phi = atan2(yr, xr);
//if (fabs(xr) > 1e-14 && fabs(yr) > 1e-14)

// printf("r: %.3f phi: %.3f, x: %d y: %d\n",r,phi/M_PI,x,y);

double dr = R0 + r * R1 + r*r*r*R3 + r*r*r*r*r*R5;
r += dr;

double x2 = r*cos(phi);
double y2 = r*sin(phi);
x2 = x2*xc+xc;
y2 = (y2)*yc+yc;

// printf("r: %.3f x: %d y: %d\n",r,(int)x2,(int)y2);

if (x2 >=0 && x2 < w && y2 >= 0 && y2 < h) {
if (use_cache) mapping[ind(x,y,w)] = ind(x2,y2,w);
buf[ind(x,y,w)] = val[ind(x2,y2,w)];
} else {
buf[ind(x,y,w)] = 0.1;
if (use_cache) mapping[ind(x,y,w)] = -1;
}

} }
lerp(1.0,buf,val,n);
free(buf);
}


//VAL SCALING
void val_radial(int w,int h,double p, double ang2, double *val) {
double SQRT2 = pow(2,0.5);
int xc = w/2.0;
int yc = h/2.0;
for (int x=0;x<w;x++) { for (int y=0;y<h;y++) {
double xr = (double)x - xc;
double yr = (double)y - yc;
double l = xr*xr + yr*yr;
double l2 = ((double)xc*xc + (double)yc*yc);
double r = limit(0,1.0,sqrt(l/l2));

double ang = atan(yr/xr);
if (ang <= 0) ang += M_PI;
if (yr < 0) ang += M_PI;

ang += ang2;

double ang_n = ang/(2*M_PI); //normalized

// double r = limit(0,1.0,sqrt(0.5 + 0.5*sin(0.3*p*l)));
//double r = limit(0,1.0,sqrt(0.5 + 0.5*sin(300*p*l/l2)));

// l = sin(0.1*p*l)+p*p;
// double r = limit(0,1.0,sqrt(l/l2));


// val[ind(x,y,w)] = (sin(ang)*sin(ang)) * (1.0 - pow(r,p));
val[ind(x,y,w)] = (1.0 - pow(r,p));
}
}
}

void cross_hairs(double *v, double os, int w, int h) {
double x = 0.5*w;
double y = 0.5*h;
double f = 0.10;
val_radial(w,h,0.7,0,v);
scale(1.2,0,v,w*h);
val_line(x,0,x,h,f,v,w,h);
val_line(0,y,w,y,f,v,w,h);

f = 0.24;
int ml = 6;
for (int i=0;i<ml;i++) {
x = ((double)i/ml)*w;
val_line(x,0,x,h,f,v,w,h);
y = ((double)i/ml)*h;
val_line(0,y,w,y,f,v,w,h);
}

ml = 30;
f = 0.2;
double t = 0.015;
for (int i=0;i<ml;i++) {
x = ((double)i/ml)*w;
val_line(x,(0.5-t)*h,x,(0.5+t)*h,f,v,w,h);
y = ((double)i/ml)*h;
val_line((0.5-t)*w,y,(0.5+t)*w,y,f,v,w,h);
}
val_lens(-0.15,v,w,h);
val_boxblur(w*os*.002,h*os*0.002,v,w,h);

scale(0.3,0,v,w*h);
}


double *visual_lchan;
double *visual_rchan;
int visual_numsamples;

void visual() {
printf("Generating visual..\n");
if (visual_lchan==0 || visual_rchan==0 || visual_numsamples==0) {
printf("Error: Audio data must be present for visualizer\n");
return;
}

int w2 = 720;
int h2 = 720;

double os=1; //oversampling rate during processing
double dos=1; //sampling rate before rendering (leave at one for downsample from 'os' value)
int w = w2*os;
int h = h2*os;

img_t img;
// init_img(&img,3,w2,h2);
init_img(&img,3,w2*dos,h2*dos);

int n = w*h;

double *v[4]; //buffer storing images
for (int i=0;i<4;i++) {
v[i] = make_bufn(n);
}

int fn=0;
double t=0;

init_drawline();

double dur=(double)visual_numsamples/SAMPLE_RATE;

double x0 = 0.5*w;
double y0 = 0.5*h;
int div = 2000;

cross_hairs(v[3],os,w,h);

while (t <= dur) {
fill(0,v[0],n);

int imax = SAMPLE_RATE/24;
for (int i=0;i<imax;i++) {

int in = (t-1.0/24)*SAMPLE_RATE + i;
if (in < 0) continue;
double x = visual_lchan[in];
double y = visual_rchan[in];
x = (x+1.0)*0.5 * w;
y = (1-(y+1.0)*0.5) * h;

val_line(y0,x0,y,x,(double)i/imax,v[0],w,h);
x0 = x;
y0 = y;
}
double bl = ((double)w/400)*os;

val_boxblur(1*bl,1*bl,v[0],w,h);


add(v[0],v[1],n);


val_boxblur(3*bl,3*bl,v[1],w,h);
scale(0.2,0,v[1],n);

add(v[0],v[1],n);

lerp(1.0,v[1],v[2],n);


for (int j=0;j<3;j++)
val_boxblur(3*bl,3*bl,v[2],w,h);

scale(2.4,0,v[2],n);
add(v[1],v[2],n);


normalize(v[2],n);
softclip(1.8,0,v[2],n);
add(v[3],v[2],n);
normalize(v[2],n);

write_frame(w,h,os/dos,v[2],&img,&fn,&t);
}

}

int only_metadata;
void audio() {
printf("Generating audio..\n");
int tone_n = 20.0 * SAMPLE_RATE;
double *tone[5];
double *bar[5];
for (int i=0;i<5;i++) {
tone[i] = make_bufn(tone_n);
fill(0,tone[i],tone_n);

bar[i] = make_bufn(tone_n);
fill(0,bar[i],tone_n);
}

int env_n = 10.0 * SAMPLE_RATE;
double *env[10];
for (int i=0;i<10;i++) {
env[i] = make_bufn(env_n);
fill(0,env[i],env_n);
}



int chan_n;
double *chan[2];
chan[0] = make_buf(90.0,&chan_n);
chan[1] = make_bufn(chan_n);

double *efx[2];
efx[0] = make_bufn(chan_n);
efx[1] = make_bufn(chan_n);

int p =0;


double bpm = drand(80,300); // PARAM for tempo orig 80
bpm = drand(80,300); //run it twice because first floating point number from PRNG is always close to zero.
printf("PARAM bpm: %.3f\n",bpm);
int beats_per_bar = randp()% 4 + 2; //PARAM orig 5;
printf("PARAM beats_per_bar: %d\n",beats_per_bar);
double beat_dur = 60/bpm;
int bs = beat_dur*SAMPLE_RATE;

double fpat[128];
//make frequency pattern

int fpatn = 10; //PARAM number of notes originaly 10
double dice = drand(0,1.0);
fpatn = 5+randp()%6;
printf("PARAM fpatn: %d\n",fpatn);

/*
if (dice < 0.1) fpatn = 10; //10% are 10
else if (dice < 0.2) fpatn = 5; //10% are 5
else if (dice < 1.0) fpatn = 4; //rest are 4
*/

int inharmonic = 0;
if (randp() % 5==0) inharmonic=1; //PARAM inharmonic 20% are inharmonic

printf("PARAM inharmonic: %d\n",inharmonic);
int note_changes = 0;
dice = drand(0,1.0);
// printf("dice: %.3f\n",dice);
if (dice < 0.1) note_changes = 3; //PARAM note changes
else if (dice < 0.2) note_changes = 2;
else if (dice < 0.4) note_changes = 1;
else note_changes = 0;
printf("PARAM note_changes: %d\n",note_changes);

int scale_index=randp() % 10; //PARAM what type of scale to chose

int scale_maxes[] = {7,4,5,6,5,5,6,5,6,6};
int scale_max = scale_maxes[scale_index];

double freq = drand(60,140); //PARAM 60 - 120
printf("PARAM freq: %.3f\n",freq);
double plow = freq;
double phigh = limit(0,1000,freq*16);
char *scale_names[] = {"major", //0
"add9", //1
"minadd9", //2
"0:2:3:6:7:12", //3
"m7",//4
"diminished", //5
"lydian", //6
"quartal",//7
"harmonic minor", //8
"pentatonic major"}; //9

printf("PARAM scale_index: %d %s\n",scale_index,scale_names[scale_index]);
int scales[][7] = {{0,2,4,5,7,9,11}, //major
{0,2,4,7,0,0,0}, //maj add9
{0,2,3,7,12,0,0}, //minor
{0,2,3,6,7,12,0}, //something
{0,3,7,10,12,0,0}, //m7
{0,3,6,9,12,0,0}, //diminished
{0,2,4,6,7,12,0}, //lydian
{0,2,5,10,12,0,0}, //quartal
{0,3,7,8,11,12,0}, //harm minor
{0,2,4,7,9,12,0}, // pentatonic major
{0,4,7,0,4,7,0}
};
if (only_metadata) {
//only print out parameters
exit(0);
}

for (int i=0;i<fpatn;i++) {

int note = scales[scale_index][randp() % scale_max];
double f= freq * pow(2,(double)note/12.0);
// f = freq;
fpat[i] = f * pow(2,irand(0,3));
printf("fpat[%d] = %.3f\n",i,f);
}

//filter pattern
double filtpat[128];
//make frequency pattern
int filtpatn = 8;
for (int i=0;i<filtpatn;i++) {
filtpat[i] = drand(200,800);
}

double pat[] = {5,3,5,1,
4,0,3,1};
int patn = 8;

normalize(pat,patn);

double apat[] = {.003,.05,.05};
int apatn = 3;

double filtapat[] = {.02,.01,.01};
int filtapatn = 3;

double pwpat[] = {0.9,0.9,0.5};
int pwpatn = 3;

double rpat[] = {.001};
int rpatn = 1;

double pat_beat_dur;

p = bs*4;

int p2 = 0;

int total_samples = bs*beats_per_bar*12;
int kmeasure=1;

double dur = 50.0; //50 seconds
total_samples = dur *SAMPLE_RATE;

while (p2 < total_samples && p2 < chan_n) {
printf("measure: %d\n",kmeasure++);

//screw with tempo and bar length


beat_dur = 60/bpm;
bs = beat_dur*SAMPLE_RATE;


p = bs * beats_per_bar; //duration of measure in samples
double bar_dur = (double)p/SAMPLE_RATE; //duration of measure in seconds
int pat_beat_div;
trapz(200,env[6],p);

pat_beat_div = 3; //irandp(1,4)*2;
pat_beat_dur = beat_dur/pat_beat_div; //fraction of beats, 1/4 is sixteenth note in 4/4

//modify the patterns
for (int j=0;j<4;j++) {
int i = randp() % patn;
// pat[i] = limit(0.2,0.5,pat[i] + irandp(-2,3)*0.25);
}

for (int j=0;j<2;j++) {
int i = randp() % pwpatn;
pwpat[i] = limit(0.5,0.99,pwpat[i] + drand(-0.2,0.2));
}

//PARAM inharmonicity
if (inharmonic) {
//changes fundamental by an incremental amount. this will affect any note changes in the phrases
freq += drand(-200,200);
freq = limit(100,2000,freq);
}
// end inharmonicity

//Each measure, possibly change some notes in the scale
for (int i=0;i<note_changes;i++) { //PARAM note changes per measure
int note = scales[scale_index][randp() % scale_max];
double f= freq * pow(2,(double)note/12.0);
int z = randp() %fpatn;
fpat[z] = f * pow(2,irand(0,3));
}

//change octave of notes
for (int i=0;i<3;i++) {
int z = randp() % fpatn;
fpat[z] *= pow(2,irand(-1,2));
if (fpat[z] < plow) {
fpat[z] *= 2;
}else if (fpat[z] > phigh) {
fpat[z] /= 2;
}
}

//amplitude envelope
fill(0,env[0],env_n);
beat_env(pat_beat_dur, //dur of each entry in 'pat'
1.0, //pw
pat, //pattern
patn, //how many entries in pat to cycle through
beats_per_bar * pat_beat_div, //number of entries to take from 'pat'
env[0],env_n);



//frequency envelope
fill(0,env[2],env_n);
beat_env(pat_beat_dur, //dur of each entry in 'pat'
1.0, //pw
fpat, //pattern
fpatn, //how many entries in pat to cycle through
beats_per_bar * pat_beat_div, //number of entries to take from 'pat'
env[2],env_n);


//amplitude pattern
fill(0,env[1],env_n);
beat_env2(pat_beat_dur, //dur of each entry in 'pat'
pat,patn,
apat,apatn,
pwpat,pwpatn,
rpat,rpatn,
beats_per_bar * pat_beat_div, //number of entries to take from 'pat'
env[1],env_n);

//filter env pattern
fill(0,env[3],env_n);
beat_env2(pat_beat_dur, //dur of each entry in 'pat'
filtpat,filtpatn,
filtapat,filtapatn,
pwpat,pwpatn,
rpat,rpatn,
beats_per_bar * pat_beat_div, //number of entries to take from 'pat'
env[3],env_n);


fill(0,tone[2],p+100);
fill(0,tone[0],p+100);
fill(0,tone[3],p+100);
fill(0,tone[4],p+100);

//make phrase,add new material
if (kmeasure==2 || irand(0,2)==1) {
//reseed it with saw wave
//saw_tone_t(env[2],p, 0, bar_dur, 1.0, tone[1],p); //env[2] is freq env
sin_tone_t(env[2],p, 0, bar_dur, 1.0, tone[1],p); //env[2] is freq env
//reslp(800,0.5,tone[1],p);

mult(env[1],tone[1],p);

// reslp_t(env[3],p,1,tone[1],p);

lerp(0.5,tone[1],bar[0],p);
mult(env[6],bar[0],p);
normalize(bar[0],p+100);
msmix(bar[0],0,0,bar[1],bar[2],p+100);
}

if (1) {
//then keep FM modulating it
for (int i=0;i<1;i++) {
lerp(1.0,env[2],env[4],p); //carrier
scale(1.5,0,env[4],p);

line(0,1000,bar_dur,env[5],p);
sin_lfo(drand(0.5,1.5),200,200,env[5],p);

modulate_t2(env[5],
env[4], //carrier
bar[0],p,
tone[0]); //dst

mult(env[6],tone[0],p);
lerp(0.1,tone[0],bar[1],p);

sin_lfo(drand(0.2,2),200,200,env[5],p);
modulate_t2(env[5],
env[4], //carrier
bar[0],p,
tone[0]); //dst

mult(env[6],tone[0],p);
lerp(0.1,tone[0],bar[2],p);


mult(env[1],bar[1],p);
mult(env[1],bar[2],p);

// reslp_t(env[3],p,1,bar[1],p);
//reslp_t(env[3],p,1,bar[2],p);

bandpass(100,8000,bar[1],p+100); //originally 8000
bandpass(100,8000,bar[2],p+100);

normalize(bar[1],p+100);
normalize(bar[2],p+100);

lerp(1.0,tone[0],bar[0],p);
bandpass(50,10000,bar[0],p+100);

// reslp(8000,0.1,bar[0],p+100);
normalize(bar[0],p+100);

}
}

fill(0,&bar[1][p],100);
fill(0,&bar[2][p],100);
fadeout(.001,bar[1],p);
fadeout(.001,bar[2],p);
//copy the stereo measure into master buffer
add(bar[1],&chan[0][p2],p+100);
add(bar[2],&chan[1][p2],p+100);
p2 += p;

// if (kmeasure > 2) break;
}

p = p2;

normalize_stereo(chan[0],chan[1],p);

lowpass(10000,chan[0],p);
lowpass(10000,chan[1],p);
softclip(1.8,0,chan[0],p);
softclip(1.8,0,chan[1],p);
lowpass(10000,chan[0],p);
lowpass(10000,chan[1],p);

// now mess with this loop. cut it up, distort it, repeat it

int endbuf = 6.0*SAMPLE_RATE;
p += endbuf;
ffdelay(beat_dur,-3,200,400,chan[0],efx[0],p);
ffdelay(beat_dur*2,-3,200,800,chan[1],efx[1],p);

add(efx[0],chan[0],p+endbuf);
add(efx[1],chan[1],p+endbuf);

normalize_stereo(chan[0],chan[1],p);
compress_inplace_stereo(.01,400,3,chan[0],chan[1],p);
softclip(1.1,0,chan[0],p);
softclip(1.1,0,chan[1],p);

printf("Calculating FFT convolution reverb...\n");
reverb(3.0,-43, .02, chan[0],efx[0],p);
reverb(3.0,-43, .02, chan[1],efx[1],p);

add(efx[0],chan[0],p);
add(efx[1],chan[1],p);

stereomixdown(chan[0],chan[1],p);
free(efx[0]);
free(efx[1]);
for (int i=0;i<10;i++) {free(env[i]); }

visual_numsamples = p;
visual_lchan = chan[0];
visual_rchan = chan[1];
}



unsigned long seed2long(char *seed,int i) {
char sbuf[9];
memcpy( sbuf, &seed[2+i*8], 8 );
sbuf[8] = '\0';
unsigned long v = strtol(sbuf,NULL,16);
return v;
}


int defaultSeed() {
if (strcmp(seed,"0x2c68ca657f166e1e7a5db2266ee1edca9daf2da5ecc689c39c1d8cfe0b9aad2e")==0) {
return 1;
}
return 0;
}

#define rand() myrand()

int myrand() {
for (int i=0;i<20;i++) {
printf("WARNING!!!! do not use rand();\n");
}
exit(0);
return 0;
}

int main(int argc, char *argv[]) {
if (defaultSeed(seed)) {
printf("Safety check. You must first replace 'seed' with your tokens seed value.\n");
return 0;
}

only_metadata = 0;
if (argc > 1 && strcmp(argv[1],"-m")==0) {
//pass "-m" switch to print out only meta data
only_metadata = 1;
}
//stream 0 seeds from first 4 bytes
randp_seed(0,seed2long(seed,0));
// printf("%lu\n",seed2long(seed,0));

audio();
visual(); //visual is deterministic, uses audio as input

return 0;
}