The topic of this chapter is the Discrete Cosine Transform (DCT), which is used in MP3 and related formats for compressing music; JPEG and similar formats for images; and the MPEG family of formats for video.
DCT is similar in many ways to the Discrete Fourier Transform (DFT), which we have been using for spectral analysis.
Once we learn how DCT works, it will be easier to explain DFT.
Here are the steps to get there:
We’ll start with the synthesis problem: given a set of frequency components and their amplitudes, how can we construct a wave?
Next we’ll rewrite the synthesis problem using NumPy arrays. This move is good for performance, and also provides insight for the next step.
We’ll look at the analysis problem: given a signal and a set of frequencies, how can we find the amplitude of each frequency component? We’ll start with a solution that is conceptually simple but slow.
Finally, we’ll use some principles from linear algebra to find a more efficient algorithm. If you already know linear algebra, that’s great, but I’ll explain what you need as we go.
Suppose I give you a list of amplitudes and a list of frequencies, and ask you to construct a signal that is the sum of these frequency components.
Using objects in the thinkdsp module, there is a simple way to perform this operation, which is called synthesis:
This example makes a signal that contains a fundamental frequency at 100 Hz and three harmonics (100 Hz is a sharp G2).
It renders the signal for one second at 11,025 frames per second and puts the results into a Wave object.
Conceptually, synthesis is pretty simple.
But in this form it doesn’t help much with analysis, which is the inverse problem: given the wave, how do we identify the frequency components and their amplitudes?
This function looks very different, but it does the same thing.
Let’s see how it works:
np.outer computes the outer product of ts and fs.
The result is an array with one row for each element of ts and one column for each element of fs.
Each element in the array is the product of a frequency and a time, \(f t\).
We multiply args by \(2 \pi\) and apply cos, so each element of the result is \(\cos (2 \pi f t)\).
Since the ts run down the columns, each column contains a cosine signal at a particular frequency, evaluated at a sequence of times.
np.dot multiplies each row of M by amps, element-wise, and then adds up the products.
In terms of linear algebra, we are multiplying a matrix, M, by a vector, amps.
In terms of signals, we are computing the weighted sum of frequency components.
The following figure shows the structure of this computation.
#TODO: Figure here
Each row of the matrix, M, corresponds to a time from 0.0 to 1.0 seconds; \(t_n\) is the time of the \(n\)th row.
Each column corresponds to a frequency from 100 to 400 Hz; \(f_k\) is the frequency of the \(k\)th column.
I labeled the \(n\)th row with the letters \(a\) through \(d\); as an example, the value of \(a\) is \(\cos [2 \pi (100) t_n]\).
The result of the dot product, ys, is a vector with one element for each row of M.
The \(n\)th element, labeled \(e\), is the sum of products:
\[e = 0.6 a + 0.25 b + 0.1 c + 0.05 d\]
And likewise with the other elements of ys.
So each element of ys is the sum of four frequency components, evaluated at a point in time, and multiplied by the corresponding amplitudes.
And that’s exactly what we wanted.
We can use the code from the previous section to check that the two versions of synthesize produce the same results.
The biggest difference between ys1 and ys2 is about 1e-13, which is what we expect due to floating-point errors.
Writing this computation in terms of linear algebra makes the code smaller and faster.
Linear algebra provides concise notation for operations on matrices and vectors.
For example, we could write synthesize like this:
\[\begin{split}\begin{aligned} M &=& \cos (2 \pi t \otimes f) \\ y &=& M a \end{aligned}\end{split}\]
where \(a\) is a vector of amplitudes, \(t\) is a vector of times, \(f\) is a vector of frequencies, and \(\otimes\) is the symbol for the outer product of two vectors.
Now we are ready to solve the analysis problem.
Suppose I give you a wave and tell you that it is the sum of cosines with a given set of frequencies.
How would you find the amplitude for each frequency component? In other words, given ys, ts and fs, can you recover amps?
In terms of linear algebra, the first step is the same as for synthesis: we compute \(M = \cos (2 \pi t \otimes f)\).
Then we want to find \(a\) so that \(y = M a\); in other words, we want to solve a linear system.
NumPy provides linalg.solve, which does exactly that.
The first two lines use ts and fs to build the matrix, M.
Then np.linalg.solve computes amps.
But there’s a hitch.
In general we can only solve a system of linear equations if the matrix is square; that is, the number of equations (rows) is the same as the number of unknowns (columns).
In this example, we have only 4 frequencies, but we evaluated the signal at 11,025 times.
So we have many more equations than unknowns.
In general if ys contains more than 4 elements, it is unlikely that we can analyze it using only 4 frequencies.
But in this case, we know that the ys were actually generated by adding only 4 frequency components, so we can use any 4 values from the wave array to recover amps.
For simplicity, I’ll use the first 4 samples from the signal.
Using the values of ys, fs and ts from the previous section, we can run analyze1 like this:
n=len(fs)amps2=analyze1(ys1[:n],fs,ts[:n])amps2
array([0.6 , 0.25, 0.1 , 0.05])
And sure enough, amps2 is the sequence of amplitudes we started with.
This algorithm works, but it is slow.
Solving a linear system of equations takes time proportional to \(n^3\), where \(n\) is the number of columns in \(M\).
We can do better.
One way to solve linear systems is by inverting matrices.
The inverse of a matrix \(M\) is written \(M^{-1}\), and it has the property that \(M^{-1}M = I\). \(I\) is the identity matrix, which has the value 1 on all diagonal elements and 0 everywhere else.
So, to solve the equation \(y = Ma\), we can multiply both sides by \(M^{-1}\), which yields:
\[M^{-1}y = M^{-1} M a\]
On the right side, we can replace \(M^{-1}M\) with \(I\):
\[M^{-1}y = I a\]
If we multiply \(I\) by any vector \(a\), the result is \(a\), so:
\[M^{-1}y = a\]
This implies that if we can compute \(M^{-1}\) efficiently, we can find \(a\) with a simple matrix multiplication (using np.dot).
That takes time proportional to \(n^2\), which is better than \(n^3\).
Inverting a matrix is slow, in general, but some special cases are faster.
In particular, if \(M\) is orthogonal, the inverse of \(M\) is just the transpose of \(M\), written \(M^T\).
In NumPy transposing an array is a constant-time operation.
It doesn’t actually move the elements of the array; instead, it creates a “view” that changes the way the elements are accessed.
Again, a matrix is orthogonal if its transpose is also its inverse; that is, \(M^T = M^{-1}\).
That implies that \(M^TM = I\), which means we can check whether a matrix is orthogonal by computing \(M^TM\).
So let’s see what the matrix looks like in synthesize2.
In the previous example, \(M\) has 11,025 rows, so it might be a good idea to work with a smaller example:
# suppress scientific notation for small numbersnp.set_printoptions(precision=3,suppress=True)amps=np.array([0.6,0.25,0.1,0.05])N=4.0time_unit=0.001ts=np.arange(N)/N*time_unitmax_freq=N/time_unit/2fs=np.arange(N)/N*max_freqys=synthesize2(amps,fs,ts)
amps is the same vector of amplitudes we saw before.
Since we have 4 frequency components, we’ll sample the signal at 4 points in time.
That way, \(M\) is square.
ts is a vector of equally spaced sample times in the range from 0 to 1 time unit.
I chose the time unit to be 1 millisecond, but it is an arbitrary choice, and we will see in a minute that it drops out of the computation anyway.
Since the frame rate is \(N\) samples per time unit, the Nyquist frequency is N/time_unit/2, which is 2000 Hz in this example.
So fs is a vector of equally spaced frequencies between 0 and 2000 Hz.
With these values of ts and fs, the matrix, \(M\), is:
You might recognize 0.707 as an approximation of \(\sqrt{2}/2\), which is \(\cos \pi/4\).
You also might notice that this matrix is symmetric, which means that the element at \((j, k)\) always equals the element at \((k, j)\).
This implies that \(M\) is its own transpose; that is, \(M^T = M\).
But sadly, \(M\) is not orthogonal.
If we compute \(M^TM\), we get:
But if we choose ts and fs carefully, we can make \(M\) orthogonal.
There are several ways to do it, which is why there are several versions of the discrete cosine transform (DCT).
One simple option is to shift ts and fs by a half unit.
This version is called DCT-IV, where “IV” is a roman numeral indicating that this is the fourth of eight versions of the DCT.
If you compare this to the previous version, you’ll notice two changes.
First, I added 0.5 to ts and fs.
Second, I canceled out time_units, which simplifies the expression for fs.
Some of the off-diagonal elements are displayed as -0, which means that the floating-point representation is a small negative number.
This matrix is very close to \(2I\), which means \(M\) is almost orthogonal; it’s just off by a factor of 2. And for our purposes, that’s good enough.
Because \(M\) is symmetric and (almost) orthogonal, the inverse of \(M\) is just \(M/2\).
Now we can write a more efficient version of analyze:
Starting with amps, we synthesize a wave array, then use dct_iv to compute amps2.
The biggest difference between amps and amps2 is about 1e-16, which is what we expect due to floating-point errors.
Finally, notice that analyze2 and synthesize2 are almost identical.
The only difference is that analyze2 divides the result by 2. We can use this insight to compute the inverse DCT:
definverse_dct_iv(amps):returndct_iv(amps)*2
inverse_dct_iv solves the synthesis problem: it takes the vector of amplitudes and returns the wave array, ys.
We can test it by starting with amps, applying inverse_dct_iv and dct_iv, and testing that we get back what we started with.
thinkdsp provides a Dct class that encapsulates the DCT in the same way the Spectrum class encapsulates the FFT.
To make a Dct object, you can invoke make_dct on a Wave.
The result is the DCT of a triangle wave at 400 Hz, as shown in the following figure.
The values of the DCT can be positive or negative; a negative value in the DCT corresponds to a negated cosine or, equivalently, to a cosine shifted by 180 degrees.
make_dct uses DCT-II, which is the most common type of DCT, provided by scipy.fftpack.
importscipy.fftpack# class Wave:defmake_dct(self):N=len(self.ys)hs=scipy.fftpack.dct(self.ys,type=2)fs=(0.5+np.arange(N))/2returnDct(hs,fs,self.framerate)
The results from dct are stored in hs.
The corresponding frequencies, computed as in Section [dctiv]{reference-type=“ref” reference=“dctiv”}, are stored in fs.
And then both are used to initialize the Dct object.
Dct provides make_wave, which performs the inverse DCT.
We can test it like this:
wave2=dct.make_wave()max(abs(wave.ys-wave2.ys))
8.881784197001252e-16
The biggest difference between ys1 and ys2 is about 1e-16, which is what we expect due to floating-point errors.
make_wave uses scipy.fftpack.idct:
# class Dctdefmake_wave(self):n=len(self.hs)ys=scipy.fftpack.idct(self.hs,type=2)/2/nreturnWave(ys,framerate=self.framerate)
By default, the inverse DCT doesn’t normalize the result, so we have to divide through by \(2N\).
In this chapter I claim that analyze1 takes time proportional to \(n^3\) and analyze2 takes time proportional to \(n^2\).
To see if that’s true, run them on a range of input sizes and time them.
In IPython, you can use the magic command %timeit.
If you plot run time versus input size on a log-log scale, you should get a straight line with slope 3 for analyze1 and slope 2 for analyze2.
You also might want to test dct_iv and scipy.fftpack.dct.
I’ll start with a noise signal and an array of power-of-two sizes
defanalyze1(ys,fs,ts):"""Analyze a mixture of cosines and return amplitudes. Works for the general case where M is not orthogonal. ys: wave array fs: frequencies in Hz ts: times where the signal was evaluated returns: vector of amplitudes """args=np.outer(ts,fs)M=np.cos(PI2*args)amps=np.linalg.solve(M,ys)returnamps
64
153 μs ± 0 ns per loop (mean ± std. dev. of 1 run, 10,000 loops each)
128
538 μs ± 0 ns per loop (mean ± std. dev. of 1 run, 1,000 loops each)
256
2.38 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 100 loops each)
512
14.4 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 100 loops each)
1024
73.5 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 10 loops each)
2048
416 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)
4096
2.93 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)
2.3866223708327747
The estimated slope is close to 2, not 3, as expected.
One possibility is that the performance of np.linalg.solve is nearly quadratic in this range of array sizes.
Here are the results for analyze2:
Show code cell content
Hide code cell content
defanalyze2(ys,fs,ts):"""Analyze a mixture of cosines and return amplitudes. Assumes that fs and ts are chosen so that M is orthogonal. ys: wave array fs: frequencies in Hz ts: times where the signal was evaluated returns: vector of amplitudes """args=np.outer(ts,fs)M=np.cos(PI2*args)amps=np.dot(M,ys)/2returnamps
64
6.58 μs ± 0 ns per loop (mean ± std. dev. of 1 run, 100,000 loops each)
128
7.46 μs ± 0 ns per loop (mean ± std. dev. of 1 run, 100,000 loops each)
256
8.46 μs ± 0 ns per loop (mean ± std. dev. of 1 run, 100,000 loops each)
512
10.6 μs ± 0 ns per loop (mean ± std. dev. of 1 run, 100,000 loops each)
1024
14.5 μs ± 0 ns per loop (mean ± std. dev. of 1 run, 100,000 loops each)
2048
24.1 μs ± 0 ns per loop (mean ± std. dev. of 1 run, 10,000 loops each)
4096
41.3 μs ± 0 ns per loop (mean ± std. dev. of 1 run, 10,000 loops each)
0.43247225316739907
This implementation of dct is even faster.
The line is curved, which means either we haven’t seen the asymptotic behavior yet, or the asymptotic behavior is not a simple exponent of \(n\).
In fact, as we’ll see soon, the run time is proportional to \(n \log n\).
The following figure shows all three curves on the same axes.
One of the major applications of the DCT is compression for both sound and images.
In its simplest form, DCT-based compression works like this:
Break a long signal into segments.
Compute the DCT of each segment.
Identify frequency components with amplitudes so low they are inaudible, and remove them.
Store only the frequencies and amplitudes that remain.
To play back the signal, load the frequencies and amplitudes for each segment and apply the inverse DCT.
Implement a version of this algorithm and apply it to a recording of music or speech.
How many components can you eliminate before the difference is perceptible?
thinkdsp provides a class, Dct that is similar to a Spectrum, but which uses DCT instead of FFT.
As an example, I’ll use a recording of a saxophone:
To compress a longer segment, we can make a DCT spectrogram.
The following function is similar to wave.make_spectrogram except that it uses the DCT.
Show code cell content
Hide code cell content
fromthinkdspimportSpectrogramdefmake_dct_spectrogram(wave,seg_length):"""Computes the DCT spectrogram of the wave. seg_length: number of samples in each segment returns: Spectrogram """window=np.hamming(seg_length)i,j=0,seg_lengthstep=seg_length//2# map from time to Spectrumspec_map={}whilej<len(wave.ys):segment=wave.slice(i,j)segment.window(window)# the nominal time for this segment is the midpointt=(segment.start+segment.end)/2spec_map[t]=segment.make_dct()i+=stepj+=stepreturnSpectrogram(spec_map,seg_length)
Now we can make a DCT spectrogram and apply compress to each segment:
