Notes:

  • This program uses the C99 function nearbyint and the C99 type int16_t. If you use an older compiler which doesn't have them, you have to remove the #include <stdint.h> and define them yourself (e.g., typedef short int int16_t; double nearbyint(double x) { return floor(x + 0.5); }).
  • If you use an operating system treating text files specially (e.g. Microsoft Windows), you'll have to modify the program so that it opens a file rather than writing to stdout.
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#define SRATE 44.1e3
#define FREQ 220.
#define AMP 16384.
#define LEN 5.
#define TWOPI 6.283185307179586476925286766559
#define DITHER(x) (nearbyint((x) + rand()/(RAND_MAX + 1.)	\
				 - rand()/(RAND_MAX + 1.)))

enum { Sine = 0, Triangle = 1, Square = 2, Sawtooth = 3 };
int main(void)
{
        int i, j;
	double a0;		/* (Absolute) amplitude of first harmonic */
        double rel_amp[1000];	/* Amplitude of the j-th harmonic relative to
				 * the first one */
	double tot = 0;		/* Sum of squares of relative amplitudes */
        for (j = 1; j * FREQ < SRATE/2.; j++) {
		switch (WAVEFORM) {
		case Sine:
			rel_amp[j] = (j == 1) ? 1. : 0.;
			break;
		case Triangle:
        		rel_amp[j] = (j % 4 == 1) ?  1./(j * j) :
				     (j % 4 == 3) ? -1./(j * j) : 0.; 
			break;
		case Square:
	        	rel_amp[j] = (j % 2 == 1) ? 1./j : 0.; 
			break;
		case Sawtooth:
			rel_amp[j] = (j % 2 == 1) ? 1./j : -1./j;
			break;
		}
        	tot += rel_amp[j] * rel_amp[j];
        }
        a0 = AMP/sqrt(tot);	/* Normalize absolute amplitude so that
				 * root-mean-square equals AMP */
        for (i = 0; i < LEN * SRATE; i++) {
        	double t = i/SRATE;
        	double val = 0.;
		int16_t sample;
        	for (j = 1; j * FREQ < SRATE/2.; j++)
        		val += a0 * rel_amp[j] * sin(TWOPI * FREQ * j * t);
		sample = DITHER(val);
                fwrite(&sample, sizeof sample, 1, stdout);
        }
        return 0;
}