/* Try to match pitch of input sound (adult male voice) with a bank of PLLs.
Input and output are 8ksps 8-bit unsigned.
Compile in C99 mode, e.g.
gcc -std=c99 -g -Wall pll.c -o pll
Two example ways to run it on Linux with ALSA:
sox voices.wav -t raw -r 8000 -u -1 -c 1 - | ./pll | aplay
arecord | ./pll | aplay # preferably with headphones
The best source I’ve found for learning about PLLs is chapter 12 of
Hans Camenzind’s book “Designing Analog Chips”, which can be
downloaded from
This works pretty well for the octave I’ve been trying it on, and it’s
enormously simpler and less computationally intensive than a Fourier
Transform. Each iteration of the inner loop has been unrolled and
inlined by GCC to 28 instructions, all straight-line except for two
jumps to conditionally negate, and only very simple instructions:
movl, addl, subl, leal, sarl, testb, jne. On my netbook, it runs at
about 500× faster than real time, about 4 megasamples per second,
which suggests that further trickery will be needed to run it for a
wider range of pitches on an Arduino.
Although it more or less works (at least for a limited range of my
voice) this code has some problems:
- It doesn’t work reliably, even for things right smack in the middle
of its theoretical range. Sometimes it misses a strong signal. I
don't know why. See comment at end of file for how to reproduce the
problem.
- The output “square wave” isn’t square. I don’t know why.
- The code is relatively fragile to modifications because of the
interrelationships between the bit shift in current_om, the number
of PLLs, and the low-pass filter rolloff, which together determine
the sensitivity of the system to noise, its responsiveness to
frequency changes, its ability to capture new signals, and the lock
range once a signal is captured.
- If you have more than an octave’s worth in the PLL bank (which I do,
barely), you sometimes end up with it locking onto a second or third
harmonic instead of a fundamental.
- The lock ranges of the PLLs are perhaps too wide; I've seen the
lowest one lock onto the DC component of the signal!
A “crad” is a made-up unit of angular measure designed to avoid
multiplication and division in the main loop of the algorithm.
*/
#include
enum {
pi = 0x400, /* Crads per pi radians. must be a power of 2 for & */
rate = 8000, /* Samples per second. */
iir_shift = 9, /* >> for low-pass filter rolloff at 512 samples (16Hz). */
};
typedef struct pll {
int om0; /* Const base angular frequency, omega, in crads/sample. */
int pha; /* Current phase in crads. */
int err; /* Accumulated low-pass-filtered phase error. */
int amp; /* Accumulated low-pass-filtered amplitude. */
} pll;
/* Compute the actual current angular frequency of a PLL. The bit
shift “10”, which determines how much frequency shift we get, was
determined by trial and error because I am ignorant of theory. */
static int current_om(pll f) { return f.om0 + (f.err >> 10); }
static void iir_low_pass(int *out, int in) { *out += in - (*out >> iir_shift); }
/* Update PLL state for a new input sample. */
static void run_pll(pll *p, int cc)
{
p->pha += current_om(*p);
iir_low_pass(&p->err, (p->pha ) & pi ? cc : -cc);
iir_low_pass(&p->amp, (p->pha - pi/2) & pi ? cc : -cc);
}
/* Best integer approximation of a/b. */
static int over(int a, int b) { return (a + b/2) / b; }
/* Convert Hz to an omega value in crads per sample. */
static int om_from_hz(int hz) { return over(pi * 2 * hz, rate); }
static int hz_from_om(int om) { return over(rate * om, pi * 2); }
static int min(int a, int b) { return a < b ? a : b; }
static int max(int a, int b) { return a > b ? a : b; }
static void debug_dump(pll ps[], int len, int max_pll)
{
for (int ii = 0; ii < len; ii++) {
fprintf(stderr, "%3d", hz_from_om(current_om(ps[ii])));
fprintf(stderr, ii == max_pll ? "* " : " ");
}
fprintf(stderr, "%6d\n", ps[max_pll].amp);
}
int main()
{
pll ps[4] = {{0}};
int cc; /* Last input sample. */
int max_pll; /* Index of PLL with max amplitude. */
int ph = 0; /* Phase of output signal. */
int tt = 0; /* Time, for periodic debugging output. */
/* Set PLL bank base frequencies, covering an octave around the
range of my voice’s usual fundamental. GCC compiles this loop to
four movl instructions with -O5. */
#define LEN(x) (sizeof(x)/sizeof((x)[0]))
for (int ii = 0; ii < LEN(ps); ii++) {
ps[ii].om0 = om_from_hz(88 + ii * 88.0 / (LEN(ps)-1) + 0.5);
}
while ((cc = getchar()) != EOF) {
max_pll = 0;
for (int ii = 0; ii < LEN(ps); ii++) {
run_pll(&ps[ii], cc);
if (ps[ii].amp > ps[max_pll].amp) max_pll = ii;
}
/* Generate a square wave using omega from the strongest PLL for
our output, attenuated when the amplitude is too low. */
ph += current_om(ps[max_pll]);
putchar(ph & pi ? 0 : max(0, min(255, ps[max_pll].amp >> 5)));
if (0 == tt++ % 1024) debug_dump(ps, LEN(ps), max_pll);
}
return 0;
}
/* Here’s a strong 124Hz signal of my voice that it misses, even if
repeated indefinitely, in a base64-encoded .wav file:
UklGRiIDAABXQVZFZm10IBAAAAABAAEAQB8AAEAfAAABAAgAZGF0Yf4CAAB+fnVveXBuW1BLRmuC
nameoY6XjHdvQk5EZ4iPs5etn6GfeG9BQz1Tdn2ci5+fpKyQimdsaXqRkJ+RmpieqJqUenVsbXNp
ZVJMQ0pvh6CjmpqIloNzZD9OR3SImrGbsZullXdqPUc/YX2Ko5Opoq6tk4pnbm+Gl5mklqGfqqyj
mYWDen+BeXNfW1Roj6G7r6+lmqaFhmBVW2GVmL20sbijsYiAWkFLSnuDn6Gdr6O2oZd8bHR1lJWh
nJqgm6mdloN2cmpxamhaTkpGbYafq5yhipaKeXJJVEdth5Syma6ZnZZ4cEFHQVt7hZ+Ro56jppCK
a29rfZCToZafmqClmpR9em9tcGdoVVBITnWHpqWgnYiYfHpiRlNHfYeirZywk6SGdmI8SD5ndoyd
j6aYqp2Qg2Z0aIaNlJuNnZGinJWNdnhnbmtkYk1MP1R1iKaYmouEjHByTkZHTX5/ppydo4+deHJP
PURDb3WVlZWjlqmSjnVmcG2Ni5yWlJ2WqJuZh3h4a3VtbGBRUEprhZ2uoaaPko94dU5WUGqOlbii
r6Wdonx5SUhJWH+EppikpqOvlZNwcXN7mJWol56fn6yamYF7dXF8bW5cVVFUfI2oqJ6eiJiAd2dJ
WE9+i6Cwm62UoIt1ZT1JQGV5iJ6PpJqkn42EZnNuh5eYoZGemKGklpF5enB0dGpnUVBKV3yLp6Cg
nIqafHhdTFZRhIeoq5+uk6SBdl4+TENwepGbkaeVqJeLfGFuaIuSmp2Qm5Ogm5CFcXZqc25oX01L
Qk1tgZmalZWDjndqWT9MRXV9lKSSpIqVfWVWOkU8YnGCk4agjZ2RfnRXZ2KDipCThpaNnJSJf2xw
aHJsZ15LS0BPb4OcmJmUh5N2b1hFUlGChqWpnq2WpYR2Wz5NR3aAmJ6ZrJywm499ZnJylZympJqm
oK6jmYl7fXqEfnprW1hQWnaOoKShoJSbhnZpTFpXf4ufr52wmKOLc2A8TEVugZChlKmbqZyJfWBv
bo2YnaOWpZ6qopOJ
Run it as follows to see this program utterly fail to detect the
strong signal even when repeated infinitely:
while : ; do sox 124hz-8bit.wav -t raw - ; done | ./pll | aplay
*/