Nyquist's Theorem in Action with Audio::Wav

[jeffa]

So, how long have CD's been around now? Since the early eighties, right? Right, in 1982 the first digital audio 5-inch CD discs were marketed - Philips/Sony finalized their standard in 1980, they had begun work in 1969.

So, how long has digital audio been around? Shortly before that, right? Nope, digital audio was born in the early 20th century when Alec H. Reeves invented pulse-code modulation PCM at the International Telephone and Telegraph Co. in France. PCM wouldn't be used until 1962 after the first practical integrated circuits were developed. The first analog sound wave was digitized at 8000 hertz.

But, the notion of digital audio and the notion of 'digitizing' and analog wave has been around since the early 1830's, when Samual Morse was expermenting with telegraph relays.

So, who is Nyquist and what does he have to do with all of this ob-trivia?

Henry Nyquist was a chap that worked at Bell Telephone Laboratories in the early 20th century. In 1928, he discovered that to acurately represent an analog wave digitally, that analog wave needs to be sampled by at least TWICE the amount of the highest perceived frequency that is desired for playback. I mentioned above that the first analog wave was digitized at 8KHz, this means that highest perceived frequency would only be half of that, 4KHz - but this perfect for the human voice, which is contained in the range of roughly 90 to 1200 Hz.

So, what happens when you break Nyquist's Theorem? Frequency aliasing - new frequecies appear and distort the original signal. The same phenoma appears when you watch the blades of a helicopter rotate - you seem to see new 'blades' which rotate the other direction at a much slower frequency.

According to Nyquist's Theorem, if S is the sampling rate, F is a frequency greater than 2S, and N is some integer, then a new sampled frequency is also created at:

F_a = ± NS ± F

Where Perl Comes In

So, let's demonstrate Nyquist's Theorem with the CPAN module Audio::Wav. Our sampling rate will be 44.1KHz at 16 bits - the CD standard. According to the formula above, if we try to sample a 36KHz wave, we should hear an 8KHz aliased wave:

S - F = F_a
44 - 36 = 8

And here is the code to demonstrate. Two.wav files will be created (each one is two seconds - one at 36KHz and the other at 8KHz) - play them back to back in your favorite audio player for a comparison. Have fun, and for an extra bonus, call a friend on your telephone and play the results for them - they won't be able to hear it, because your conversion is being filtered to prevent aliasing.

---

use strict;
use Audio::Wav;
use Math::Complex;

use constant TWO_PI => pi() * 2;

my $wav = Audio::Wav->new();
my $sample_rate = 44100;
my $bits_sample = 16;

write_out('nyquist_good.wav',8000,2);
write_out('nyquist_bad.wav',36000,2);

sub write_out {
   my $filename = shift;
   my $write = $wav->write($filename, {
      bits_sample => $bits_sample,
      sample_rate => $sample_rate,
      channels    => 1,
   });
   add_sine($write,@_);
   $write->finish();
}

sub add_sine {
   my ($write,$hz,$length) = @_;
   my $max_no = (2 ** $bits_sample) / 2;
   $length   *= $sample_rate;

   for my $pos (0..$length) {
       my $time = ($pos / $sample_rate) * $hz;
       my $val  = sin(TWO_PI * $time);
       my $samp = $val * $max_no;
       $write->write($samp);
   }
}

I updated the code with jryan's suggestion, he also told me to remind you folks to turn down your volume before you listen, as an 8KHz sine wave is not a rather pleasing tone. :)

---

References:

• Audio::Wav
• http://history.acusd.edu/gen/recording/digital.html
• http://history.acusd.edu/gen/recording/notes.html
• http://www.tnt-audio.com/topics/frequency_e.html
• Principles of Digital Audio
• Code

jeffa

L-LL-L--L-LL-L--L-LL-L--
-R--R-RR-R--R-RR-R--R-RR
F--F--F--F--F--F--F--F--
(the triplet paradiddle)

[Zaxo]

Good writeup, thanks.

It's interesting to note that 44.1kHz == 44100 == 2**2 * 3**2 * 5**2 * 7**2. All those small factors are good for sample rate conversions, reducing interpolation for natural harmonics, and fast Fourier transforms.

After Compline,
Zaxo

[jryan]

Instead of use constant PI => (22 / 7) * 2;, give the core use Math::Complex; a try. Just use the module, and a sub which defines PI more accurately is defined (via 4*arctan(1)). The sub is even called pi, so you wouldn't even have to change most of your code.

Other than that, an EXCELLENT node, nice work!

[tilly]

Some notes about Nyquist's Theorem which I have mentioned in the past.

The result was first suggested by Nyquist, but the first proof was due to Claude Shannon, for which reason you will encounter it with one, the other or both names attached.

It also gives rise to a piece of popular terminology. The width of frequencies you are using determines the rate at which you can send data points, which is how that rate became known as your bandwidth. Technically, though, this is a misnomer because your frequency band limits how fast you can send analog datapoints, not digital ones. But given measurement accuracy, it does put a practical limit on how fast you can send bits.

Now there is an obvious question. Where does the factor of 2 come from?

Well a naive understanding can be gained from saying that you need both sin and cos coefficients.

Someone with more math, comfortable with thinking about complex exponentials (ie exp(i*$x) = cos($x) + i*sin($x)) can remember that if you are going up to frequency f then you are actually dealing with complex exponentials in the range -f to f. So the actual width of the band is twice the frequency range. (Of course out of a positive and a negative exponential term you produce a sin and a cos term, so this is really the same as what we had before.)

Someone with a lot more math can actually produce a very good heuristic as follows. We want to understand the space of L2 functions (L2 means that the square is integrable) which are limited in complex exponentials to the frequency range from $a to $b. Well the first step is to find ourselves a nice orthonormal basis. Well an orthonormal basis is going to remain orthonormal under a Fourier transform, so take the Fourier transform of our functions and in the frequency domain we are talking about L2 functions which are 0 except inside of the interval from $a to $b.

Well L2 functions with support in the range from $a to $b are well understood, and an obvious orthonormal basis is the usual complex exponential basis (which is the usual Fourier series). So apply the inverse Fourier transform to that and (after much calculation) you find that they come back with a basis in terms of functions:

$C*exp(i*($b-$a)/2) * sinc($k + $x*($b-$a))

(I forget the constant $C, and I don't have time to work it out, sorry) where sinc($x) is defined as sin(PI * $x)/(PI * $x). Which very conveniently happens to be 1 at 0, and 0 at every other integer.

So we have our orthonormal basis for this band limited signal, and if we sample at time intervals of 1/($b-$a) we are measuring our basis on the nose. (Because one basis term is 1, and the rest are all 0 there.) How convenient!

Now back to our problem. In the case where we have a maximum frequency f, the range is from -f to f, and voila! The Nyquist rate falls out! (This would likely convince a physicist, but a full formal proof is a lot harder.)

Now if you understood that explanation in full, you have far more math and/or physics than you need. If you kind of got a sense, that is what I was hoping for. If it left you feeling dazed and confused, sorry. This one requires a lot of background that isn't very relevant to Perl, but I find this kind of stuff fun. (And now most of you probably think of me as a sick puppy, but that's OK. I was a math geek long before I ever heard of Perl...)

UPDATE
I miswrote a + as a * in the statement of the basis. My apologies for the typo. I also added a word of clarification about how sinc relates to samples.