/* Spline wave synthesis on Linux. One instrument (exponentially * decaying quartic spline waves), using some of the syntax of the * GW-BASIC PLAY command * . * * At present, according to Cachegrind, this needs about 37 * instructions on amd64 per output sample. And it’s about 900 * instructions. It does not use floating point. It needs about * 16+6N bytes of RAM for a * So spline wave synthesis should probably be doable on a * smallish microcontroller. * * Spline waves are spline approximations of sine waves with two * spline nodes per cycle, either at the crest and trough or at the * zero-crossings, depending on whether the order of the spline is * even or odd. You can compute them efficiently by repeatedly * integrating square waves. Audacity tells me that the spline waves * coming out of this program have pretty minimal harmonic content; * the first significant harmonic is the third, with the following * amplitude relative to the fundamental: * * quadratic spline: 3rd harmonic is 27dB below fundamental * cubic spline: 3rd harmonic is 37dB below fundamental * quartic spline: 3rd harmonic is 47dB below fundamental * * I sort of lied above: I said “repeatedly integrating” but didn’t * specify where the constants of integration came from. That’s * because I don’t know, so I’m running each stage of integration * through a single-pole high-pass filter (i.e. exponential decay) to * keep things from blowing up too badly. (It still blows up a bit, * and you can hear the glitches.) You might wonder if perhaps * multiple stages of high-pass fiters might affect the 3rd-harmonic * numbers above, making the 3rd harmonic relatively larger. And yes, * with a sufficiently small time constant, this does happen: using * *buf -= *buf >> 3 for the high-pass filter, I do get only 24dB * instead of 47dB for the quartic, >>4 gives us 37dB, >>5 already * gives us 45dB. Pumping it all the way up to >>9 breaks it * completely, but at >>8 we only get 48dB out of the quartic. * * In its current state this barely demonstrates spline wave * synthesis! There are a couple of tricks left to pull: * * - As I mentioned before, I need a better solution for keeping DC * offsets from blowing up. Integrating once, your constant of * integration (the initial value of the integrator) is a DC offset * on the output; at worst, you'll get clipping/overflow where you * wouldn’t. Integrating twice or more, your constant of * integration becomes a linear drift on the output, which will * produce a sawtooth wave modulated by your signal (which, let me * tell you, sounds like screaming demented demon babies). What do * you do about that? Can you choose the constant of integration * correctly? Do you just subtract a little of the integrated * signal from the integrand, thus sort of low-pass filtering the * integrand to eliminate the DC offset? My current solution seems * ad-hoc. * * - The simple low-pass filtering means that every time we switch the * frequency of the driver, we get a big low-frequency spike, which * can easily cause clipping. * * - I think something is a bit dodgy about my amplitude adjustment * code in square_wave_compensate. Why power-1 instead of power? * It seems to work slightly better. */ /* * Remaining to-do: * - add harmonics for different instruments * - add remaining PLAY command bits: <#+ONPT MF MB MN ML MS X? * - clean up; some things are dirty. */ #include #include #include /* C99 */ #include /* See comment on `half_periods` when changing this. */ #define SAMPLING_RATE 48000 typedef struct { uint16_t half_period; /* amplitude initially contains the positive desired amplitude of the square wave. Every time the countdown reaches zero, it gets negated. */ int16_t amplitude; uint16_t countdown; } square_oscillator; static inline int16_t square_wave_output(square_oscillator *o) { return o->amplitude; } static inline void square_wave_next(square_oscillator *o) { if (o->countdown--) return; o->countdown = o->half_period; o->amplitude = -o->amplitude; } /* To set up a square wave that will, after integration some number of * times, produce a spline wave with a desired amplitude, we need to * give a higher amplitude to higher-frequency square waves than to * lower-frequency square waves. How much higher depends on the * degree of the spline. */ static inline void square_wave_compensate(square_oscillator *o, uint8_t power) { for (int ii = 0; ii < power-1; ii++) { o->amplitude = o->amplitude * 64 / o->half_period; } } static inline FILE * open_audio_output() { /* Stringification trick from GCC manual section. What was X3J11 * thinking? */ #define STRING(x) XSTRING(x) #define XSTRING(x) #x return popen("aplay -q -f S16_LE -r " STRING(SAMPLING_RATE), "w"); /* return popen("sox -t raw -r " STRING(SAMPLING_RATE) " -b 16 -e signed-integer -L -c 1 - music.wav", "w"); */ } // Arrangement of Simple Gifts in F major from // ; // should probably consider replacing this with the original from // . char *simple_gifts = "L8 CC F4 FG AFAB- >C4 >CB- A4 GF" "G4 G4 G4 F4 GAGE C4 C4 FEFG A4 GG" "A4 B4- >C4 A4 G4 GG A4 AG F4 FF F2" ">C2 A4. G AB-AG F4. G A4 AB- >C4 A4 G4 GA G4 C4" "F2 F4. G A4 AB- >C4 B-A G4 G4 A4 AG F4 F4 F2"; /* Here we have the square-wave half-periods for the chromatic scale from middle C (C₄) for an octave, inclusive (including C₅): python -c 'print [int(round(48000 / (261.6 * 2**(ii/12.0)) / 2)) for ii in range(13)]' Sadly, the rounding error here is quite large, although my own ear is bad enough that it doesn’t bother me. Also, changing SAMPLING_RATE requires recomputing this. */ const uint8_t half_periods[] = {92, 87, 82, 77, 73, 69, 65, 61, 58, 55, 51, 49, 46}; /* starting here with A: A B C D E F G */ const int8_t n_halfsteps[] = {9, 11, 0, 2, 4, 5, 7}; /* Doing our own buffering makes the whole program use about 12% less * instruction fetches than using stdio for buffering. Note we’re * doing it kind of half-assed and not flushing the buffer at the end. * Using fputc instead of putc or doing our own buffering was costing * 75% of the program’s runtime. */ static inline bool output_sample(FILE *output, int16_t sample) { static char outbuf[256]; static char *p = outbuf; if (p == outbuf + sizeof(outbuf)) { if (1 != fwrite(outbuf, sizeof(outbuf), 1, output)) return false; p = outbuf; } *p++ = sample & 0xff; *p++ = sample >> 8; return true; } typedef int32_t integrator; static inline int32_t integrate(integrator *buf, int32_t sample) { *buf += sample; /* Exponential decay to eliminate DC offsets causing overflow. I * arrived at 5 through trial and error. */ *buf -= *buf >> 5; return *buf; } /* Actually more like a score iterator. */ typedef struct { char *pos; uint8_t length; /* of each note, by default, as reciprocal */ bool pending_octave_up; /* > */ } score; typedef struct { uint8_t half_period; uint8_t duration; bool eof; } note; static inline note next_note(score *buf) { if (!buf->length) buf->length = 8; /* 64th-notes I guess. */ note rv = {0, buf->length, false}; int8_t pitch; /* in half-steps above middle C */ bool got_note = false; /* char note_letter; */ /* Loop over the characters. We have to keep going until we find * the next note, in case we need to find a flat ('-'). */ char *p = buf->pos; for (; *p; p++) { char c = *p; if ('a' <= c && c <= 'z') c -= 'a' - 'A'; if ('A' <= c && c <= 'G') { if (got_note) break; got_note = true; /* note_letter = c; */ pitch = n_halfsteps[c - 'A']; if (buf->pending_octave_up) pitch += 12; buf->pending_octave_up = false; } else if (c == '>') { /* This could be either before or after the note we’re * returning. */ buf->pending_octave_up = true; } else if (c == '.') { rv.duration += rv.duration / 2; } else if ('0' <= c && c <= '9') { rv.duration = 64 / strtol(p, &p, 10); p--; /* Because we’re about to increment it. */ } else if (c == 'L') { buf->length = 64 / strtol(p+1, &p, 10); p--; } else if (c == '-') { pitch--; } } if (!*p && !got_note) { rv.eof = true; return rv; } buf->pos = p; //printf("Playing %c (%d) for %d\n", note_letter, pitch, rv.duration); rv.half_period = half_periods[pitch]; return rv; } static inline void play_music(FILE *output, char *music) { score buf = {music}; int bpm = 120; int beat_samples = 60 * SAMPLING_RATE * 4 / bpm / 64; /* printf("a beat is %d samples\n", beat_samples); */ integrator triangle = {0}, parabola = {0}, cubic = {0}, quartic = {0}; square_oscillator osc = {0}; for (;;) { note next = next_note(&buf); if (next.eof) break; osc.half_period = next.half_period; /* Reboost amplitude. */ osc.amplitude = osc.amplitude > 0 ? 256 : -256; /* Make amplitude independent of frequency. */ square_wave_compensate(&osc, 4); for (int ii = 0; ii < beat_samples * next.duration; ii++) { square_wave_next(&osc); /* Exponential amplitude decay: */ if (!(ii%256)) osc.amplitude -= osc.amplitude >> 6; int16_t sample = integrate(&quartic, integrate(&cubic, integrate(¶bola, integrate(&triangle, square_wave_output(&osc)))) ) >> 11; if (!output_sample(output, sample)) { perror("output sample"); return; } } } } int main(int argc, char **argv) { play_music(open_audio_output(), argc > 1 ? argv[1] : simple_gifts); return 0; }