in development
2018-09-19

A long time ago (I wanna say around mid-2015) I wrote the following C program:

``````#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#define k double
#define l unsigned short
#define m *((int *) &a[30][5])
#define n rand() %
#define o ((k) rand()) / RAND_MAX
#define q(a) ((char *) &h)[a]
#define r 1870552195.42848
#define s 264600
#define t 573300
#define u return
#define v 7069.358259
#define w a[i]

k a[31][6]={{2.5139212861031446e+88,2.0846937317106930e+112,1.2882297587186030e-231,
6.6113133109285920e+22,2.0602858415292070e-289,-1.3817238424510467e+82}},b;k c(k a,k
b,k d){u(1-d)*a+d*b;}k d(k a,k b,k c){u c<s?(1-c/s)*a/v+b*c/r:b/v;}k e(k a,k b){u b<
s?a*b* (
.1238359
+7e-8):(       b<485100         )?a*(327      *0144+67         ):(1-(b-
485100)/       88200)*a         *32767;}       i,j,f,g,       h; p(int
h){u b=m       <s?(1-m/         s)*h[a][        0]/v+h[a     ][1]*(m)
/r:h[a][       1]/v,g=(         fabs((((         b+a[h][3   ])/(6.28
+.00319)       )-floor(         0.5+(a[h          ][3]+b)/ 6.28319)
)*2)-0.5       )* e(a[h][2],m),h[a] [3]=           fmod(h[a][3]+d(
h[a][0],       h[a][1],m),6.28319),f=((l            )(g*c(a[h][4]
,a[h][5]       ,((k)m/t))/30))*65536+(l)           (g*c(a[h][4],a[
h][5],1-       ((k)m/t)         )/30),(*          ((int*)& a[++h][5
]))++,f;       } main (         ){srand(         time(0))   ;for(i=0
; i <054       ;)fputc(         ((char *        )a)[7-(i     %8)+(i++
&0370)],       stdout);         a:switch       (i){case       44: i=0;
default:       w[0]= (n         200)+200      ,w [1] =         pow(2, n

4+(n 2)*(7.0/12.0))*42,w[1]+=((n 012)/11.0)*(w[1]/0xc8),w[4]=w[5]=o,w[2]=log(w[1]) /
3.5,w[3]=0,i++;goto a;case 30:w[5]=0;break;}i-=30;b:switch(i){default:h=0;for(j=0;j<
30;)h+=p(j++)*2;printf("%c%c%c%c",q(0),q(1),q(2),q(3),i++);goto b; case t: break;};}
``````

I had originally intended this as an IOCCC entry, but it's been long enough since I wrote it (and I think someone else did it better) that I'm going to make something new for the next one.

This code, when run, will produce a WAV file that sounds roughly like the THX Deep Note (hence the shape of the source code). I've mostly just obfuscated the syntax of the program here (another reason to write a new program for the next IOCCC), but I've been pretty thorough with that.

Before we start, a little about the structure of the Deep Note: it consists of about 30 oscillators with randomized pitches that fade in and converge on a E5 chord after six seconds, which is then sustained for an additional five seconds before fading out over two seconds.

Without further ado, let's jump in.

The headers and `#defines` should be fairly self-explanatory; the defines mostly exist to compress the code and make it easier to format.

``````k a[31] = {{ ... }}, b;

double phase_offset;
union {
double d[6];
struct osc {
double f1;
double f2;
double gain;
double phase;
double pan1;
double pan2;
} o;
} oscs[31] = {{...}};
``````

Yeah, I know the initializer for `oscs` isn't valid anymore. What are you gonna do about it?

``````k c(k a,k b,k d) {
u(1-d)*a+d*b;
}

double lerp(double a, double b, double t) {
return (1 - t) * a + t * b;
}
``````

Nothing special here. `c` simply linearly interpolates its arguments. It's pretty hard to obfuscate that further without slowing the code down beyond what's reasonable.

``````k d(k a,k b,k c) {
u c<s
? (1-c/s)*a/v+b*c/r
: b/v;
}

double chirp(double f1, double f2, double sample) {
if(sample < 264600) {
return (1 - sample/264600)*f1/7069.358259 + f2*sample/1870552195.42848;
} else {
return f2/7069.358259;
}
}

double chirp(double f1, double f2, double sample) {
if(sample < 44100 * 6) {
return lerp(f1, f2, sample / (44100 * 6)) / (44100 / 2*pi);
} else {
return f2 / (44100 / 2*pi);
}
}
``````

Lotta magic numbers here. This function takes a start frequency and an end frequency and returns the phase offset between time `sample` and time `sample+1`. In the main loop, the result of this is added to each oscillator's instantaneous phase every sample.

Let's simplify some of these expressions.

264600, first of all, is just 44100 Hz (the sample rate) * 6 s (the length of the chirp)

Next, I'll factor a 1/7069.358259 out of the 3rd line, making it ```return ((1 - sample / 264600) * f1 + f2 * sample / 264600) / 7069.358259```. This is now clearly an inlined linear interpolation, so let's un-inline that to get ```return lerp(f1, f2, sample / 264600) / 7069.358259```.

Now we're getting somewhere, so the only remaining mystery is that 1/7069.358259 constant, which is just 44100/2pi. [It actually isn't 44100/2pi. It's equal to 44100/6.23819, and 2*pi equals 6.28319. I'm pretty sure this is just an error in the code and the correct value is 7018.72775.] Anyways, 44100/2pi is a conversion factor from Hz to radians/sample, because that's what the rest of the program expects from this function.

The conditional is just to stop the chirp at `f2` after six seconds. In the original program, it's a ternary instead.

``````k e(k a,k b) {
u b<s
? a*b* (.1238359+7e-8)
: (b<485100)
? a*(327*0144+67)
: (1-(b-485100)/88200)*a*32767;
}

double fade(double gain, double sample) {
if(sample < 264600) {
return gain * sample * 0.12383597;
} else {
if(sample < 485100) {
return gain * 32767;
} else {
return (1 - (sample - 485100) / 88200) * gain * 32767;
}
}
}

double fade(double gain, double sample) {
double vol;
if(sample < 44100 * 6) {
vol = sample / (44100 * 6);
} else if(sample < 44100 * 11) {
vol = 1.0;
} else {
vol = 1 - (sample - (44100 * 11)) / (44100 * 2);
}
return vol * gain * 32767;
}
``````

Even more magic numbers. As before, 264600 is the split point between the chirp and the sustain. 485100 is 11 seconds, which is the start of the fade out. We're using 16-bit signed samples here, so "full scale" is a gain of 32767.

This function essentially implements a piecewise-linear function that from sample 0 to 264600 goes from 0 to 32767, then holds at 32767 until sample 485100, and finally returns to 0 at sample 573300, and the whole function is multiplied by the `gain` parameter.

I've essentially factored the `gain * 32767` ouf of the conditional in the final deobfuscation. Hopefully this makes the structure of the function clear enough.

``````i,j,f,g,h;

int i, j, sample_stereo, amplitude, sample_total;
``````

You don't actually need to specify a type in C (unless you're compiling with `-Werror`). These are just a couple global variables that we'll use later. I'm pretty sure most of them don't need to be global.

``````p(int h) {
u b=m<s
? (1-m/s)*h[a][0]/v+h[a][1]*(m)/r
: h[a][1]/v,
g=(fabs((((b+a[h][3])/(6.28+.00319))-floor(0.5+(a[h][3]+b)/6.28319))*2)-0.5)
*e(a[h][2],m),
h[a][3]=fmod(h[a][3]+d(h[a][0],h[a][1],m),6.28319),
f=((l)(g*c(a[h][4],a[h][5],((k)m/t))/30))*65536
+(l)(g*c(a[h][4],a[h][5],1-((k)m/t))/30),
(*((int*)&a[++h][5]))++,
f;
}
``````

Yeah, this is the big one. I'm going to go through it in pieces rather than deobfuscate the whole thing in one go.

Go ahead and look back at what I redefined `a` as way at the start of the article. This is where some of that comes into play.

``````p(int h) {
u b = m<s
? (1 - m/s)*h[a][0]/v+h[a][1]*(m)/r
: h[a][1]/v,

if(oscs[30].o.time < 264600) {
phase_offset = (1 - oscs[30].o.time) * oscs[index].o.f1 / 7069.358259
+ oscs[index].o.f2 * oscs[30].o.time / 1870552195.42848;
} else {
phase_offset = oscs[index].o.f2 / 7069.358259;
}
``````

Oscillator number 30 isn't used, so the code reuses its `pan2` field for the current sample. I'm going to replace `oscs[30].o.pan2` with `current_sample` from here for clarity.

This is simply an inline of `chirp()` from earlier:

``````unsigned int advance(int index) {
phase_offset = chirp(oscs[index].o.f1, oscs[index].o.f2, current_sample);
``````

The original code for this section uses a ternary for the condition (which lets us avoid repeating the `b =`) and ends with a comma, which, if you aren't familiar with it (comma is a rather disused part of C syntax), essentially joins two or more expressions and evaluates to whatever the last one is. This function is made into one expression with the comma, which is why we aren't worrying about that `u` (`return`) at the start yet. It also interchanges `a[h]` and `h[a]`, which is a fairly common obfuscation technique.

``````g=(fabs((((b+a[h][3])/(6.28+.00319))-floor(0.5+(a[h][3]+b)/6.28319))*2)-0.5)
*e(a[h][2],m),
``````

For the deobfuscation, we're going to define:

``````double saw(double t) {
return fabs(t/6.28319 - floor(0.5 + t/6.28319)*2) - 0.5;
}
``````

which produces a sawtooth wave from -1 to 1, `t` in radians. With this definition, the original code becomes

``````amplitude = saw(phase_offset + oscs[index].o.phase)
``````

which calculates one sample of one oscillator. This line is where the actual audio synthesis happens.

``````h[a][3]=fmod(h[a][3]+d(h[a][0],h[a][1],m),6.28319),

oscs[index].o.phase = fmod(oscs[index].o.phase
+ chirp(oscs[index].o.f1, oscs[index].o.f2, current_sample),
6.28319);
``````

This line increments the phase of the oscillator and wraps it at 2*pi. The `chirp` call here can be replaced with the previously calculated `phase_offset`:

``````oscs[index].o.phase += phase_offset;
oscs[index].o.phase = fmod(oscs[index].o.phase, 6.28319);

f=((l)(g*c(a[h][4],a[h][5],((k)m/t)))/30))*65536
+(l)(g*c(a[h][4],a[h][5],1-((k)m/t))/30),

sample_stereo = ((unsigned short) (amplitude
* lerp(oscs[index].o.pan1, oscs[index].o.pan2, (double) current_sample / 573300) / 30)) * 65536
+ ((unsigned short) (amplitude
* lerp(oscs[index].o.pan1, oscs[index].o.pan2, 1 - ((double) current_sample / 573300)) / 30));
``````

The output is in stereo, and this is where the stereo happens. `advance()` returns an unsigned int consisting of the two packed single-channel 16-bit samples that make up one stereo sample. 573300 is the total length in samples of the audio, and each channel gets divided by 30 because there are 29 oscillators, so dividing by 30 makes clipping mathematically impossible. The multiplication by 65536 is really a left shift by 15. I'll simplify the panning calculation a bit as well.

``````double pan = lerp(oscs[index].o.pan1, oscs[index].o.pan2, (double) current_sample / 573300);
unsigned short sample_left = amplitude * pan / 30;
unsigned short sample_right = amplitude * (1 - pan) / 30;
sample_stereo = sample_left << 15 | sample_right;
``````

I've modified the panning code in this last version slightly and I'm pretty sure it's more correct. The original code pans incorrectly.

``````(*((int*)&a[++h][5]))++,

index++;
*((int *) &(oscs[index].o.pan2))++;
``````

The point of this code is to increment the `pan2` field of oscs[30]. Note that when `pan2` is being interpreted as a sample number it's being bitwise reinterpreted as an `int`, so we make sure to do that with the increment. This lets us get away with incrementing all 30 values here because twiddling the lower few bits of a `double` isn't going to make a noticeable impact on the panning.

``````      f;
}

return sample_stereo;
}
``````

`f`, being the last expression in the chain of commas, becomes the return value.

One last block to get through.

``````main() {
srand(time(0));
for(i=0; i < 054; )
fputc(((char *)a)[7-(i%8)+(i++&0370)],stdout);
a:
switch(i) {
case 44: i=0;
default: w[0] = (n 200)+200,
w[1] = pow(2, n 4+(n 2)*(7.0/12.0))*42,
w[1] += ((n 012)/11.0)*(w[1]/0xc8),
w[4]=w[5]=o,
w[2]=log(w[1]) / 3.5,
w[3]=0,
i++;
goto a;
case 30: w[5]=0; break;
}
i-=30;
b:
switch(i) {
default: h=0;
for(j=0; j<30; )
h+=p(j++)*2;
printf("%c%c%c%c",q(0),q(1),q(2),q(3),i++);
goto b;
case t:
break;
};
}
``````

Who needs control flow constructs when you have `switch` and `goto`?

``````int main(void) {
srand(time(0));
for(i = 0; i < 44; i++) {
putchar(((char *) oscs) [7 - (i%8) + (i & 248)]);
}
``````

This mess just writes the first 44 bytes of `oscs` to stdout; they make up the header for the WAV file. The strangeness in the index is just to reorder the bytes a bit.

``````  a:
switch(i) {
case 44: i=0;
default: w[0] = (n 200)+200,
w[1] = pow(2, n 4+(n 2)*(7.0/12.0))*42,
w[1] += ((n 012)/11.0)*(w[1]/0xc8),
w[4]=w[5]=o,
w[2]=log(w[1]) / 3.5,
w[3]=0,
i++;
goto a;
case 30: w[5]=0; break;
}
``````

This is a disguised `for` loop:

``````for(i = 0; i < 30; i++) {
oscs[i].o.f1 = rand() % 200 + 200;
oscs[i].o.f2 = pow(2, rand() % 4 + (rand() % 2) * (7.0 / 12.0)) * 42;
oscs[i].o.pan1 = oscs[i].o.pan2 = ((double) rand()) / RAND_MAX;
oscs[i].o.gain = log(oscs[i].o.f1) / 3.5;
oscs[i].o.phase = 0;
}
oscs[30].pan2 = 0;
``````

`rand() % 200 + 200` is the inital frequency; somewhere between 200 and 400 Hz. The whole bit with the `pow` and the 42 is the final frequency. This picks some harmonic of 42 Hz (the fundamental frequency of the final chord) and either shifts it up a 5th (7 / 12) or doesn't.

Both the initial and final pan are set to the same random value. I guess I didn't like how it sounded with the pan moving around. Maybe that has something to do with the panning being broken.

The gain is set so that higher-frequency oscillators are quieter, and each oscillator's phase is initialized at 0.

Finally, we set the current sample to 0. Since the bit patterns of (double) 0 and (unsigned long) 0 are the same, this behaves as expected.

``````  i -= 30;
b:
switch(i) {
default: h=0;
for(j=0; j<30; )
h+=p(j++)*2;
printf("%c%c%c%c",q(0),q(1),q(2),q(3),i++);
goto b;
case t:
break;
};
}
``````

Again, a disguised `for` loop:

``````    for(i = 0; i < 573300; i++) {
sample_total = 0;
for(j = 0; j < 30; j++) {
573300 is, again, the overall length of the file. We can simply add together the stereo samples (instead of separating it into left and right channels and adding those) because we made sure that clipping can't occurr, and we shift the sample left by one to put the split in the right place (`advance()` doesn't shift it far enough, just to make things a little more complicated). Then, it's a simple matter of outputting the sample values.