Reverse Engineering an Opera

We'll start out by looking at the input side—here's our simple 5 Hz wave again:

Figure: 5 Hz wave.

As before, we'll look at the line as the fundamental signal, whether it's a microphone picking up sound waves, an antenna picking up radio waves—whatever. The points then are where we've actually sampled, in this case at 100 Hz, and since we've sampled for one second, we have 100 data points in this 1-second snapshot of the 5 Hz wave.

To the eye, it's pretty clear that the wave is very well sampled, and this leads to the question of just how many samples we need in a given time period to properly capture a waveform. This is answered by Nyquist Theory, which gives us that, for a given sampling rate , the maximum frequency of wave that can be accurately sampled, or Nyquist frequency, is .

However, waves with frequency greater than the Nyquist, known as undersampled waves, do not just disappear. The undersampled waves will show an effect called aliasing, in which they mirror down, appearing as potentially spurious signals in a frequency spectrum.

Figure: A line spectrum sampled at 500 Hz, showing ten input frequencies spaced 30 Hz apart. Note that at the top end, where the last two inputs—270 Hz and 300 Hz—are greater than the Nyquist frequency, the last two peaks have aliased, mirroring around and showing up at 230 Hz and 200 Hz.

Here, note that the input frequencies are evenly spaced by 30 Hz, and so when we reach the Nyquist frequency, they mirror around, showing as peaks at false frequencies starting near the Nyquist. To try to grasp why this happens, let's look at some waveforms.

Figure: Why aliasing occurs. A 7 Hz wave (orange line), sampled at 10 Hz (green points) looks exactly the same as a 3 Hz wave (blue line).

If we have the orange 7 Hz signal, but only sample it at 10 Hz (i.e. on the green points), then those samples end up looking exactly the same as what we'd see from the blue 3 Hz signal. As we saw from the spectrum above, the aliased wave frequency appears at the sampling frequency minus the real wave's frequency.

When a data acquisition system is designed for it, undersampling and aliasing can be used in various ways, such as bandpass downconversion. However, for a straight sampling system, aliased signals are a huge problem, and so hardware filters will commonly be used on the input side to remove frequencies above the Nyquist.

It's worth elaborating that the Nyquist frequency is an absolute limit, and so it applies most precisely to ideal, clean, phase-aligned signals. For real-world observations even in very good conditions, it's advisable to have a Nyquist well above—4x or more—the maximum frequency of interest, and strong filters that cut off frequencies well below the Nyquist. In cases where higher harmonics, noise, and other unknowns may be pertinent to the investigation, a Nyquist dozens of times the fundamental frequency involved would not be out of the question.

In general the Nyquist frequency is also a good measure of how well a given sampling system can reproduce a signal. Modern, consumer-grade ('DVD-quality') audio is sampled at 48 kHz, and the average upper limit of human hearing is generally considered to be about 20 kHz, with the most detailed perception below 10 kHz. Thus we see that consumer-grade systems can do a reasonable job of reproducing audio; however, production-grade often goes further, sampling at 96 kHz or even 192 kHz in order to better capture and reproduce the audio signal as well as allow for more mixing and filtering techniques.

Okay, so that's an answer to what the minimum is for how many samples we need in a given time period to capture waves at a given range of frequencies. Going in the other direction, do we get anything out of oversampling? This question leads into one of the main tradeoffs that occurs on the input side of sampling for FFTs.

No matter what our sampling rate, we can technically feed as many samples as we want to the FFT—we'll just need more space in memory to store the snapshot's input and output arrays. While we'll cover the relation more in the output section, for every two input samples we add to our snapshot, we'll get one more frequency bin in the output, and therefore more frequency resolution. Thus, larger snapshots can allow us to distinguish separate waves that are very close in frequency.

However, at a given sampling rate, we'll need to sample for a longer time period in order to get those additional samples. Thus, if we're trying to make a spectrogram, the tradeoff is between resolution on the frequency axis, and resolution on the time axis: longer snapshots will give us better frequency resolution, but our view of time will reciprocally get more coarse.

This is particularly important when you consider time-limited signals with potentially changing frequencies, and even more so when these signals do not themselves have any inherent perodicity over the course of their lifetime. This includes many manmade signals; e.g, instruments in a musical piece need not necessarily repeat any specific pattern. Other examples of such signals can be found among natural radio emissions, such as whistler waves.

These waves are generated in the upper atmosphere by lightning, and propagate down magnetic field lines. They have important effects in the Earth's plasma environment, including precipitating charged particles in the radiation belts. They were named because of their audible descending tone when the radio signal is directly converted to an audio signal.

So, what happens if our snapshots are so long that they end up encompassing an entire whistler wave? With zero time resolution along the period of the whistler, we'll just see a broad band of power along the entire range of the whistler's frequencies, and we won't grasp the true nature of the wave in time at all.

Of course, if make our snapshots too narrow in time, we won't be able to resolve the different frequencies that make up the whistler, either. When possible, it's often helpful to try different snapshot sizes, yielding multiple views of the frequency vs. time structure in the data.

That said, raising the sampling rate will also raise the overall resolution possible in either or both directions, but this requires additional data storage, and possibly better sampling equipment. While this isn't terribly burdensome in the audio frequency range, equipment capable of fully sampling well above the Nyquist for high-frequency radio waves (3 MHz up) can cost tens of thousands of dollars, and can easily require a gigabyte of storage for each minute of observation.

A complication of the FFT process lies in the fact that, while we only input a time-limited snapshot to the FFT algorithm, the transform is based on math on infinitely periodic signals. The practical effect this has is as if the waveform snapshot we feed in is repeated infinitely. So, if we take the simple 5 Hz waveform from earlier, what's actually being transformed is more like...

Figure: A depiction of the theoretical view of what's being transformed in an FFT: the input snapshot (orange line) extends infinitely in both directions (blue line, for particularly small values of infinity).

This presents a problem, because in general the endpoints of our input waveform are not going to be continuous when connected end-to-end. If we take our earlier example of a three-component waveform, with 1-second snapshots we see a significant discontinuity at the endpoints.

Figure: An example of boundary discontinuities when a snapshot (orange line) is repeated (blue line). Note the vertical jump where the orange and blue lines connect.

Reproducing discontinuities like this would add a bunch of extraneous frequencies to our spectrum, potentially ruining our analysis. This generation of spurious frequencies is known as spectral leakage.

So, what can we do about this? The canonical solution is to modify our input waveform samples in a way that preserves as much of the information we're trying to extract as possible, adds as little distortion as possible, and makes the endpoints match up when the snapshots are placed end-to-end. This is done by windowing the input so that the ends go to zero in a smooth manner, using a window function . With the ends at zero, there can be no discontinuity at the periodic boundaries.

For each sample in our -sample input, will produce a number by which we multiply to get our windowed input set, :

One of the more common and simple window functions is the Hann window, which makes use of a cosine function:

Figure: The shape and effect of a Hann window on an input snapshot. On the left a snapshot of a 5 Hz wave, in the center the function for the Hann window, and right the effect of applying the window, .

Adding a window is equivalent to impressing an additional signal onto the data (and thus the spectral components required to reproduce that signal), but as long as we remain aware of this, it can be controlled far more effectively than the boundary-discontinuity noise of spectral leakage.

There are a wide variety of window functions, with wider peaks, flat regions, or different roll-offs towards zero. Each one is developed to have particular characteristics affecting the signal, with different noise levels, bandwidth, and frequency resolution. There are even some specific cases where the 'rectangular' or 'boxcar' window (equivalent in most cases to applying no windowing function at all) has the most desirable characteristics.

Let's look at five different windows on two signals: first, a single 5 Hz sine wave, and second, a 3-wave composite signal, the most powerful signal at 23.7 Hz, a signal with 1/3rd of the amplitude at 8.92 Hz, and a weak signal at 1/15th the amplitude, at 5 Hz. We'll sample at 1 kHz for pretty good frequency resolution, and zoom in on the 0-50 Hz range.

Figure: A comparison between five window functions. In the columns, left to right, no window (aka rectangular or boxcar window), the Hann window, the Blackman-Harris window, a Flat-top window, and a Planck-taper window. On the top row, the shape of the window function itself; middle row line spectra of a windowed 5 Hz wave; and on the bottom line spectra of a composite of three waves of various parameters. Note for each spectrum the distance between peaks and the noise floor, and the sharpness or spreading of the peaks.

There's a lot going on above. One of the most important things to note is that the spectra are all in a logarithmic scale with units of decibels (dB), and shifted so that the highest peak is at 0 dB. The top row of spectra (for the 5 Hz wave) are all set to the same 800 dB y-axis range, while the lower row are individually set to show the most important data. From just this relatively simple comparison, we can see a clear tradeoff exists.

While all four of the cosine-based windows have very low noise floors, approximately 600 dB below peak for the single, pure sine wave, the rectangular window is clearly the best, with a narrow, sharp peak. The more complex the construction of a cosine-based window, the more the peak spreads out, and the Planck-taper window is just shy of nonsense.

However, when we add in just two other waves, the rectangular window becomes comparatively bad due to the spectral leakage: the spreading from the peaks is huge, and the noise floor rises to more like 50 dB off-peak, almost entirely hiding the weak 5 Hz peak, which only shows up as a tiny blip above the spread. For our signal, the best window appears to be the Hann window, which shows all three peaks well, with significantly less spread than the rectangle.

Aside: Is the noise and spreading of peaks entirely due to windowing? No, the discretization/sampling process itself introduces noise and artifacts. One way to look at it is that each frequency bin covers a range of frequencies dependent on the snapshot size and sampling rate, and it's fairly unlikely that a given fundamental wave in our signal will fall directly at the center of a bin. Effectively, each frequency in a bin (or sample in time) is a little fuzzy in frequency/time, and this introduces noise.

From this small test comparison we can see at least a hint towards the need for case-by-case consideration of various window functions, depending on factors such as the signal source, probable peak separation, dynamic range, and noise levels. Particularly careful consideration is necessary in cases where the transformed data is what is stored, such as systems with limited storage (remote collection) or bandwidth (space probes) that can only store and/or transmit a selection of frequencies. In such systems, you only have one chance to get it right.

A final complication of windowing is that the window can affect the absolute powers in the result. This may be fairly obvious simply by noting that the the sum over the amplitudes in a Hann window— —comes to 0.5, i.e. half of the amplitude in a waveform is going to be damped by the Hann window. A simple and oft-used fix for this is to simply globally correct the resultant powers by dividing out . This is a useful rough approximation, but is only exact for a pure-noise signal. In general, a window function will have different effects over the range of the spectrum, and the distortion will depend upon the input signal.

Briefly, both the input and output of the Fourier transform ideally contain all of the information about the wave—amplitude, frequency, and phase—in the form of a complex number . But wait, we didn't feed it a complex number—we only had real samples, right? This is true, and the simplest way to think about it is that this means we need two samples on the input side to yield a single complex number on the output side. There are ways of sampling a signal that yield complex data, but for real samples such as ours the output of a standard DFT/FFT will be 'mirrored', i.e. one half of the output will be complex numbers arrayed by frequency, and the other half will be the same numbers, in reverse order.

The specific FFTW function we're using here, rfft(), knows that we're inputing real numbers, and so it returns the correct-order half of the output (plus a 0-frequency offset bin), and we don't have to deal with discarding half of the outputs. This also saves some processing time on the transform.

The full complex-number output is useful in various ways. Taking will yield the phase angle of each frequency component. With multiple perpendicular sampling antennas, the complex components can be mixed in various ways to find things like wave polarization and direction of arrival. For our purposes, let's look at the details of turning our complex numbers into powers.

To get the power that we're interested in, we'll have to go through a number of steps to process a given output datum .

Putting it all together yields

The last pieces we need to get some nicely formatted output are to associate our output bins to frequencies, and make sure our times are aligned correctly.

We already know that, for an FFT snapshot of size , our output will have bins. The only other thing we need is the range of frequencies we're dealing with. This range has its max at the Nyquist frequency, which at our sampling frequency will just be . Thus, we're spreading a range of to Hz over bins, meaning each bin will have an span, and so our frequencies are the range 0:F/N:F/2.

The times assigned to each of a spectrogram's snapshots are easier, but contain a choice. For each sample in our total data set, we have a from the set of all times . With samples per snapshot, to get our times we just have to go through with a step size of . The choice left to the user is where within each step you take the time for a given snapshot: from the beginning, middle, or end. For the following we'll use the beginning time, and so snapshot_times = [ T[i] for i in 1:N:length(T) ]