This HTML version of is provided for convenience, but it is not the best format for the book. In particular, some of the symbols are not rendered correctly.
You might prefer to read the PDF version.
Chapter 2 Harmonics
2.1 Triangle waves
A sinusoid contains only one frequency component, so its spectrum has only one peak. More complicated waveforms, like the violin recording, yield DFTs with many peaks. In this section we investigate the relationship between waveforms and their spectrums.
I’ll start with a triangle waveform, which is like a straight-line version of a sinusoid. Figure 2.1 shows a triangle waveform with frequency 200 Hz.
To generate a triangle wave, you can start with a thinkdsp.TriangleSignal:
class TriangleSignal(Sinusoid): def evaluate(self, ts): cycles = self.freq * ts + self.offset / PI2 frac, _ = np.modf(cycles) ys = np.abs(frac - 0.5) ys = normalize(unbias(ys), self.amp) return ys
The only difference is evaluate. As we saw before, ts is the sequence of sample times where we want to evaluate the signal.
There are many ways to generate a triangle wave. The details are not important, but here’s how evaluate works:
Here’s the code that generates Figure 2.1:
signal = thinkdsp.TriangleSignal(200) signal.plot()
wave = signal.make_wave(duration=0.5, framerate=10000) spectrum = wave.make_spectrum() spectrum.plot()
Figure 2.2 shows two views of the result; the view on the right is scaled to show the harmonics more clearly. As expected, the highest peak is at the fundamental frequency, 200 Hz, and there are additional peaks at harmonic frequencies, which are integer multiples of 200.
But one surprise is that there are no peaks at the even multiples: 400, 800, etc. The harmonics of a triangle wave are all odd multiples of the fundamental frequency, in this example 600, 1000, 1400, etc.
Another feature of this spectrum is the relationship between the amplitude and frequency of the harmonics. Their amplitude drops off in proportion to frequency squared. For example the frequency ratio of the first two harmonics (200 and 600 Hz) is 3, and the amplitude ratio is approximately 9. The frequency ratio of the next two harmonics (600 and 1000 Hz) is 1.7, and the amplitude ratio is approximately 1.72 = 2.9. This relationship is called the harmonic structure.
2.2 Square waves
class SquareSignal(Sinusoid): def evaluate(self, ts): cycles = self.freq * ts + self.offset / PI2 frac, _ = np.modf(cycles) ys = self.amp * np.sign(unbias(frac)) return ys
Like TriangleSignal, SquareSignal inherits
And the evaluate method is similar. Again, cycles is the number of cycles since the start time, and frac is the fractional part, which ramps from 0 to 1 each period.
unbias shifts frac so it ramps from -0.5 to 0.5, then np.sign maps the negative values to -1 and the positive values to 1. Multiplying by amp yields a square wave that jumps between -amp and amp.
Like a triangle wave, the square wave contains only odd harmonics, which is why there are peaks at 300, 500, and 700 Hz, etc. But the amplitude of the harmonics drops off more slowly. Specifically, amplitude drops in proportion to frequency (not frequency squared).
I have a confession. I chose the examples in the previous section carefully to avoid showing you something confusing. But now it’s time to get confused.
Figure 2.5 shows the spectrum of a triangle wave at 1100 Hz, sampled at 10,000 frames per second. Again, the view on the right is scaled to show the harmonics.
The harmonics of this wave should be at 3300, 5500, 7700, and 9900 Hz. In the figure, there are peaks at 1100 and 3300 Hz, as expected, but the third peak is at 4500, not 5500 Hz. The fourth peak is at 2300, not 7700 Hz. And if you look closely, the peak that should be at 9900 is actually at 100 Hz. What’s going on?
The problem is that when you evaluate the signal at discrete points in time, you lose information about what happened between samples. For low frequency components, that’s not a problem, because you have lots of samples per period.
But if you sample a signal at 5000 Hz with 10,000 frames per second, you only have two samples per period. That turns out to be enough, just barely, but if the frequency is higher, it’s not.
To see why, let’s generate cosine signals at 4500 and 5500 Hz, and sample them at 10,000 frames per second:
framerate = 10000 signal = thinkdsp.CosSignal(4500) duration = signal.period*5 segment = signal.make_wave(duration, framerate=framerate) segment.plot() signal = thinkdsp.CosSignal(5500) segment = signal.make_wave(duration, framerate=framerate) segment.plot()
Figure 2.6 shows the result. I plotted the Signals with thin gray lines and the samples using vertical lines, to make it easier to compare the two Waves. The problem should be clear: even though the Signals are different, the Waves are identical!
When we sample a 5500 Hz signal at 10,000 frames per second, the result is indistinguishable from a 4500 Hz signal. For the same reason, a 7700 Hz signal is indistinguishable from 2300 Hz, and a 9900 Hz is indistinguishable from 100 Hz.
In this example, the highest frequency we can measure is 5000 Hz, which is half the sampling rate. Frequencies above 5000 Hz are folded back below 5000 Hz, which is why this threshold is sometimes called the “folding frequency”. It is sometimes also called the Nyquist frequency. See http://en.wikipedia.org/wiki/Nyquist_frequency.
The folding pattern continues if the aliased frequency goes below zero. For example, the 5th harmonic of the 1100 Hz triangle wave is at 12,100 Hz. Folded at 5000 Hz, it would appear at -2100 Hz, but it gets folded again at 0 Hz, back to 2100 Hz. In fact, you can see a small peak at 2100 Hz in Figure 2.4, and the next one at 4300 Hz.
2.4 Computing the spectrum
from np.fft import rfft, rfftfreq # class Wave: def make_spectrum(self): n = len(self.ys) d = 1 / self.framerate hs = rfft(self.ys) fs = rfftfreq(n, d) return Spectrum(hs, fs, self.framerate)
The parameter self is a Wave object. n is the number of samples in the wave, and d is the inverse of the frame rate, which is the time between samples.
The result of rfftfreq, which I call fs, is an array that contains frequencies corresponding to the hs.
The Spectrum class provides two read-only properties, amps and angles, which return NumPy arrays representing the magnitudes and angles of the hs. When we plot a Spectrum object, we usually plot amps versus fs. Sometimes it is also useful to plot angles versus fs.
Although it might be tempting to look at the real and imaginary parts of hs, you will almost never need to. I encourage you to think of the DFT as a vector of amplitudes and phase offsets that happen to be encoded in the form of complex numbers.
To modify a Spectrum, you can access the hs directly. For example:
spectrum.hs *= 2 spectrum.hs[spectrum.fs > cutoff] = 0
The first line multiplies the elements of hs by 2, which doubles the amplitudes of all components. The second line sets to 0 only the elements of hs where the corresponding frequency exceeds some cutoff frequency.
But Spectrum also provides methods to perform these operations:
You can read the documentation of these methods and others at http://greenteapress.com/thinkdsp.html.
At this point you should have a better idea of how the Signal, Wave, and Spectrum classes work, but I have not explained how the Fast Fourier Transform works. That will take a few more chapters.
Solutions to these exercises are in chap02soln.ipynb.
Exercise 1 If you use Jupyter, load chap02.ipynb and try out the examples. You can also view the notebook at http://tinyurl.com/thinkdsp02.
Exercise 2 A sawtooth signal has a waveform that ramps up linearly from -1 to 1, then drops to -1 and repeats. See http://en.wikipedia.org/wiki/Sawtooth_wave
Write a class called SawtoothSignal that extends Signal and provides evaluate to evaluate a sawtooth signal.
Exercise 3 Make a square signal at 1100 Hz and make a wave that samples it at 10000 frames per second. If you plot the spectrum, you can see that most of the harmonics are aliased. When you listen to the wave, can you hear the aliased harmonics?
Exercise 4 If you have a spectrum object, spectrum, and print the first few values of spectrum.fs, you’ll see that they start at zero. So spectrum.hs is the magnitude of the component with frequency 0. But what does that mean?
Try this experiment:
Exercise 5 Write a function that takes a Spectrum as a parameter and modifies it by dividing each element of hs by the corresponding frequency from fs. Hint: since division by zero is undefined, you might want to set spectrum.hs = 0.
Test your function using a square, triangle, or sawtooth wave.
Exercise 6 Triangle and square waves have odd harmonics only; the sawtooth wave has both even and odd harmonics. The harmonics of the square and sawtooth waves drop off in proportion to 1/f; the harmonics of the triangle wave drop off like 1/f2. Can you find a waveform that has even and odd harmonics that drop off like 1/f2?
Are you using one of our books in a class?We'd like to know about it. Please consider filling out this short survey.