**1.2 Spectral decomposition**

- The most important topic in this book is
**spectral decomposition**, which is the idea that any signal can be expressed as the sum of sinusoids with different frequencies. - The most important mathematical idea in this book is the
**Discrete Fourier Transform**, or**DFT**, which takes a signal and produces its**spectrum**. The spectrum is the set of sinusoids that add up to produce the signal. - the most important algorithm in this book is the
**Fast Fourier Transform**, or**FFT**, which is an efficient way to compute the DFT.

**1.3 ~ 1.5 thinkdsp.py**

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 |
import thinkdsp import matplotlib.pyplot as plt #Signals cos_sig = thinkdsp.CosSignal(freq=440, amp=1.0, offset=0) sin_sig = thinkdsp.SinSignal(freq=880, amp=0.5, offset=0) mix = sin_sig + cos_sig #Use __add__ method wave = mix.make_wave(duration=0.5, start=0, framerate=11025) period = mix.period #get a period segment = wave.segment(start=0, duration=period*5) segment.plot() #draw plt.show() #Reading and writing Waves violin_wave = thinkdsp.read_wave('input.wav') wave.write(filename='output.wav') #Spectrum spectrum = wave.make_spectrum() spectrum.low_pass(cutoff=600, factor=0.01) spectrum.high_pass(cutoff=600, factor=0.01) spectrum.band_stop(low_cutoff=1000, high_cutoff=1500, factor=0.01) spectrum.plot() plt.show() #convert spectrum bato to wave wave = spectrum.make_wave() |

**1.6 Wave objects**

- Given a Signal, you can make a Wave.
- Given a Wave, you can make a Spectrum, and vice versa.
- A Wave object contains three attributes: ys is a NumPy array that contains the values in the signal.
- ts is an array of the times where the signal was evaluated or sampled.
- framerate is the number of samples per unit of time. The unit of time is usually seconds, but it doesn’t have to be. In one of my examples, it’s days.
- To modify a wave, you can access the ts and ys directly. For example: wave.ys *= 2, wave.ts += 1. The first line scales the wave by a factor of 2, making it louder. The second line shifts the wave in time, making it start 1 second later.

**Chapter 3. Non-periodic signals**

**3.4 Spectrogram**

- To recover the relationship between frequency and time, we can break the chirp into segments and plot the spectrum of each segment. The result is called a
**short-time Fourier Transform (STFT).** - The most common way to visualize a STFT is a
**spectrogram,**which shows time on the x-axis and frequency on the y-axis. Each column in the spectrogram shows the**spectrum**of a short segment,**using color or grayscale to represent amplitude.**

1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
import thinkdsp import matplotlib.pyplot as plt #Chirp : makes a sinusoid that sweeps linearly through a range of frequencies. signal = thinkdsp.Chirp(start=220, end=440) wave = signal.make_wave(duration=1, framerate=11025) ''' seg_length is the number of samples in each segment FFT is most efficient when the # of samples is a power of 2 ''' spectrogram = wave.make_spectrogram(seg_length=512) spectrogram.plot(high=700) plt.show() |

**3.6 Leakage**

- DFT treats waves as if they are periodic. that is, it assumes that the finite segment is a complete period from an infinite signal that repeats over all time. In practice, this assumption is often false, which creates problems.
- One common problem is discontinuity at the beginning and end of the segment. Because DFT assumes that the signal is periodic, it implicitly connects the end of the segment back to the beginning to make a loop. If the end does not connect smoothly to the beginning, the discontinuity creates additional frequency components in the segment that are not in the signal.

- Figure 3.4 (middle) shows the spectrum of this segment. Again, the peak is at 440 Hz, but now there are additional components spread out from 240 to 640 Hz. This spread is called
**spectral leakage**, because some of the energy that is actually at the fundamental frequency leaks into other frequencies. In this example, leakage happens because we are using DFT on a segment that becomes discontinuous when**we treat it as periodic.**

**3.7 Windowing**

- We can reduce leakage by smoothing out the discontinuity between the beginning and end of the segment, and one way to do that is
**windowing.** - A “window” is a function designed to transform a non-periodic segment into something that can pass for periodic. Figure 3.5 (top) shows a segment where the end does not connect smoothly to the beginning.
- Figure 3.5 (middle) shows a “Hamming window”, one of the more common window functions. No window function is perfect, but some can be shown to be optimal for different applications, and Hamming is a good, all-purpose window.
- Figure 3.5 (bottom) shows the result of multiplying the window by the original signal. Where the window is close to 1, the signal is unchanged. Where the window is close to 0, the signal is attenuated. Because the window tapers at both ends, the end of the segment connects smoothly to the beginning.
- Figure 3.4 (right) shows the spectrum of the windowed signal. Windowing has reduced leakage substantially, but not completely.
- Wave class provides window, and NumPy provides hamming like below.

1 2 |
wd = np.hamming(len(wave)) wave.window(wd) |

**Chapter 5. Autocorrelation**

**5.1 Correlation**

- In general, correlation between variables means that if you know the value of one, you have some information about the other. There are several ways to quantify correlation, but the most common is the Pearson productmoment correlation coefficient, usually denoted ρ. For two variables, x and y, that each contain N values:

- Pearson’s correlation is always between -1 and +1 (including both). If ρ is positive, we say that the correlation is positive, which means that when one variable is high, the other tends to be high. If ρ is negative, the correlation is negative, so when one variable is high, the other tends to be low.

**Chapter 6. Discrete Cosine Transform (DCT)**

**6.2 Synthesis with arrays**

- 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:
- M = cos(2πt ⊗ f) , y = Ma : where a is a vector of amplitudes, t is a vector of times, f is a vector of frequencies, and ⊗ is the symbol for the outer product of two vector

**6.3 Analysis**

- 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πt ⊗ f). Then we want to find a so that y = Ma; in other words, we want to solve a linear system. NumPy provides linalg.solve, which does exactly that.

1 2 3 4 5 |
def analyze1(ys, fs, ts): args = np.outer(ts, fs) M = np.cos(2*math.pi * args) amps = np.linalg.lstsq(M, ys) return amps |

- Like this, we can get amps from yt, fs, ts. but this algorithm that Solve a linear system of equations is so slow.

**6.4 Orthogonal matrices**

- y = Ma, so we can get a if we can compute inverse M (inverse matrix) fast.
- In particular, if M is orthogonal, the inverse of M is just the transpose of M, written . In Numpy transposing an array is a constant-time operation.
- Again, a matrix is orthogonal if its transpose is also its inverse.

**6.5 DCT-IV**

- 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.

1 2 3 4 5 6 |
def test2(): amps = np.array([0.6, 0.25, 0.1, 0.05]) N = 4.0 ts = (0.5 + np.arange(N)) / N fs = (0.5 + np.arange(N)) / 2 ys = synthesize2(amps, fs, ts) |

- Notice, First, I added 0.5 to ts and fs. Second, I canceled out time_units, which simplifies the expression for fs. (compare to previous version dct. not noted here)
- Then is almost , which means M is almost orthogonal. the inverse of M is just M/2. Now we can write a more efficient version of analyze.

1 2 3 4 5 |
def analyze2(ys, fs, ts): args = np.outer(ts, fs) M = np.cos(2*math.pi * args) amps = np.dot(M, ys) / 2 return amps |

- Instead of using np.linalg.solve, we just multiply by M/2.
- Combining test2 and analyze2, we can write an implementation of DCT-IV.

1 2 3 4 5 6 7 8 |
def dct_iv(ys): N = len(ys) ts = (0.5 + np.arange(N)) / N fs = (0.5 + np.arange(N)) / 2 args = np.outer(ts, fs) M = np.cos(2*math.pi * args) amps = np.dot(M, ys) / 2 return amps |

- Again, ys is the wave array. We don’t have to pass ts and fs as parameters; dct_iv can figure them out based on N, the length of ys. We can test it like this.

1 2 3 4 5 6 7 8 |
amps = np.array([0.6, 0.25, 0.1, 0.05]) N = 4.0 ts = (0.5 + np.arange(N)) / N fs = (0.5 + np.arange(N)) / 2 ys = synthesize2(amps, fs, ts) amps2 = dct_iv(ys) print(max(abs(amps-amps2))) |

- The result is very small (1e-16), which is what we expect due to floating-point errors.

**6.6 Inverse DCT**

- 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

1 2 3 4 5 6 7 8 9 10 11 |
def analyze2(ys, fs, ts): args = np.outer(ts, fs) M = np.cos(2*math.pi * args) amps = np.dot(M, ys) / 2 return amps def synthesize2(amps,fs, ts): args = np.outer(ts,fs) M = np.cos(2*math.pi * args) ys = np.dot(M,amps) return ys |

1 2 |
def inverse_dct_iv(amps): return dct_iv(amps) * 2 |

- inverse_dct_iv solves the synthesis problem: it takes the vector of amplitudes and returns the wave array, ys.

**Chapter 7. Discrete Fourier Transform**

**7.1 Complex exponentials**

e^φ = 1 + φ + (φ^ 2)/2! + (φ^ 3)/3! + ...

- This definition works with real numbers, imaginary numbers and, by a simple extension, with complex numbers.
- This definition works with real numbers, imaginary numbers and, by a simple extension, with complex numbers. Applying this definition to a pure imaginary number, iφ, we get

e^iφ= 1 + iφ − (φ ^2)/2! − i(φ ^3)/3! + ..

- By rearranging terms, we can show that this is equivalent to:

e^( iφ) = cos φ + i sin φ

- This formula implies that e iφ is a complex number with magnitude 1; if you think of it as a point in the complex plane, it is always on the unit circle. And if you think of it as a vector, the angle in radians between the vector and the positive x-axis is the argument, φ.

In the case where the exponent is a complex number, we have: e^( a+iφ )= (e^ a)*( e^(iφ)) = A*e^(iφ)

where '**A' is a real number that indicates amplitude and e^(iφ) is a unit complex number that indicates angle**.

- When the argument to np.exp is imaginary or complex, the result is a complex number.

1 2 3 4 5 6 7 |
phi = 1.5 z = np.exp(1j * phi) #Python uses j to represent the imaginary unit, rather than i. so 1j is just i. print(z) #(0.0707+0.997j) print(z.real)#0.0707 print(z.imag)#0.997 print(abs(z)) # 1.0, magnitude print(np.angle(z)) # 1.5 |

- This example confirms that e ^(iφ) is a complex number with magnitude 1 and angle φ radians

**7.3 The synthesis problem**

- we can create compound signals by adding up complex sinusoids with different frequencies. And that brings us to the complex version of the synthesis problem: given the frequency and amplitude of each complex component, how do we evaluate the signal?

**7.4 Synthesis with matrices**

- As we saw in Section 6.2, we can also express the synthesis problem in terms of matrix multiplication

1 2 3 4 5 6 7 |
PI2 = np.pi * 2 def synthesize2(amps, fs, ts): args = np.outer(ts, fs) M = np.exp(1j * PI2 * args) ys = np.dot(M, amps) return ys |

- Again, amps is a NumPy array that contains a sequence of amplitudes. fs is a sequence containing the frequencies of the components. ts contains the times where we will evaluate the signal.
- args contains the outer product of ts and fs, with the ts running down the rows and the fs running across the columns (you might want to refer back to Figure 6.1). Each column of matrix M contains a complex sinusoid with a particular frequency, evaluated at a sequence of ts.
- When we multiply M by the amplitudes, the result is a vector whose elements correspond to the ts; each element is the sum of several complex sinusoids, evaluated at a particular time.
- Remember that we can think of a complex number in two ways: either the sum of a real and imaginary part, x + iy, or the product of a real amplitude and a complex exponential, Ae^(iφ0) . Using the second interpretation, we can see what happens when we multiply a complex amplitude by a complex sinusoid. For each frequency, f , we have:

Ae^(iφ0) · e^( i2π f t) = Ae^(i2π f t+φ0)

- Multiplying by Ae^(iφ0) multiplies the amplitude by A and adds the phase offset φ0.
- Now that we have the more general solution to the synthesis problem – one that handles complex amplitudes – we are ready for the analysis problem.

**7.5 The analysis problem**

- The analysis problem is the inverse of the synthesis problem: given a sequence of samples, y, and knowing the frequencies that make up the signal, can we compute the complex amplitudes of the components, a?

**7.6 Efficient analysis**

- Unfortunately, solving a linear system of equations is slow. For the DCT, we were able to speed things up by choosing fs and ts so that M is orthogonal. That way, the inverse of M is the transpose of M, and we can compute both DCT and inverse DCT by matrix multiplication.
- We’ll do the same thing for the DFT, with one small change. Since M is complex, we need it to be
**unitary**, rather than orthogonal, which means that the inverse of M is the conjugate transpose of M, which we can compute by transposing the matrix and negating the imaginary part of each element. See http://en.wikipedia.org/wiki/Unitary_matrix. - The NumPy methods
*conj*and*transpose*do what we want. Here’s the code that computes M for N = 4 components:

1 2 3 4 5 |
N = 4 ts = np.arange(N) / N fs = np.arange(N) args = np.outer(ts, fs) M = np.exp(1j * PI2 * args) |

- If M is unitary, , where is the conjugate transpose of M, and I is the identity matrix. We can test wheter M is unitary like this:

1 |
MstarM = M.conj().transpose().dot(M) |

We can use this result to write a faster version of analyze1:

1 2 3 4 5 |
def analyze2(ys, fs, ts): args =np.outer(ts, fs) M = np.exp(1j * PI2 * args) amps = M.conj().transpose().dot(ys) / N return amps |

**7.7 DFT**

- As a function, analyze2 would be hard to use because it only works if fs and ts are chosen correctly. Instead, I will rewrite it to take just ys and compute freq and ts itself.
- First, I’ll make a function to compute the synthesis matrix, M:

1 2 3 4 5 6 7 8 9 10 11 12 13 |
def synthesis_matrix(N): ts = np.arange(N) / N fs = np.arange(N) args = np.outer(ts, fs) M = np.exp(1j * PI2 * args) return M def anylize3(ys): N = len(ys) M = synthesis_matrix(N) amps = M.conj().transpose().dot(ys) / N return amps |

- We are almost done; analyze3 computes something very close to the DFT, with one difference. The conventional definition of DFT does not divide by N

1 2 3 4 5 |
def dft(ys): N = len(ys) M = synthesis_matrix(N) amps = M.conj().tranpose().dot(ys) return amps |

- The inverse DFT is almost the same, except we don’t have to transpose and conjugate M, and now we have to divide through by N:

1 2 3 4 5 |
def idft(ys): N = len(ys) M = synthesis_matrix(N) amps = M.dot(ys) / N return amps |

- If I could go back in time, I might change the definition of DFT so it divides by N and the inverse DFT doesn’t. That would be more consistent with my presentation of the synthesis and analysis problems. Or I might change the definition so that both operations divide through by √ N. Then the DFT and inverse DFT would be more symmetric. But I can’t go back in time (yet!), so
**we’re stuck with a slightly weird convention.**For practical purposes it doesn’t really matter.