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

#define SAMPLE_RATE 32000
#define DURATION 57600
#define SMALL_SAMPLES 1

#if SMALL_SAMPLES
#define BITS_SAMPLE 8
#define BYTES_SAMPLE 1
#define SCALE 128
#else
#define BITS_SAMPLE 16
#define BYTES_SAMPLE 2
#define SCALE 32768
#endif

struct HEADER {
    char chunkid[4];
    int chunksize;
    char format[4];
    char subchunk1id[4];
    int subchunk1size;
    short audioformat;
    short numchannels;
    int samplerate;
    int byterate;
    short blockalign;
    short bitspersample;
    char subchunk2id[4];
    int subchunk2size;

    void init(int nsamples);
};

void HEADER::init(int nsamples) {
    strcpy(chunkid, "RIFF");
    chunksize = nsamples*BYTES_SAMPLE + sizeof(HEADER) - 8;
    strcpy(format, "WAVE");
    strcpy(subchunk1id, "fmt ");
    subchunk1size = 16;
    audioformat = 1;
    numchannels = 1;
    samplerate = SAMPLE_RATE;
    byterate = SAMPLE_RATE*BYTES_SAMPLE;
    blockalign = 2;
    bitspersample = BITS_SAMPLE;
    strcpy(subchunk2id, "data");
    subchunk2size = nsamples;
}

struct FILTER {
    double mult0, mult1, mult2, xprev1, xprev2;
    double sample_rate;

    FILTER(double sr) {
        sample_rate = sr;
        xprev1 = xprev2 = 0;
    }

    void apply(double& sample) {
        sample = mult0 * sample + mult1 * xprev1 - mult2 * xprev2;
        xprev2 = xprev1;
        xprev1 = sample;
    }

    void setup(double freq, double band) {
        double PI=3.1415926;
        double gamma, B, R;
        double cos0, angle;

        if (freq > sample_rate) freq = sample_rate;
        else if (freq < 0.) freq = 0.;
        if (band > sample_rate) band = sample_rate;
        else if (band < 0.) band = 0.;

        gamma = freq / sample_rate * PI;
        B = band / sample_rate * PI;

        R  = 1. - B / 2.;
        mult2 = R * R;
        cos0 = 2. * R * cos(gamma) / (1. + mult2);
        angle = acos(cos0);
        mult0 = (1.- mult2) * sin (angle);
        mult1 = 2. * R * cos0;
    }
};

int main() {
    int i;
    HEADER header;
    FILE* f = fopen("sound.wav", "w");
    double pi=3.1415926;
    FILTER filter(SAMPLE_RATE);
    filter.setup(4000.,500.);

    header.init(SAMPLE_RATE*DURATION);
    fwrite(&header, sizeof(header), 1, f);

    int lo=0, hi=0;
    for(i=0; i<SAMPLE_RATE*DURATION; i++) {
        double x = ((double)rand())/RAND_MAX;
        x *= 2;
        x -= 1;
        filter.apply(x);
#if SMALL_SAMPLES
        int y = 64 + (int)(64*x);
#else
        int y = (int)(x*32768);
#endif
        if (y>hi) hi=y;
        if (y<lo) lo=y;
        fwrite(&y, BYTES_SAMPLE, 1, f);
    }
    fclose(f);
}