In [1]:

```
import numpy as np
signal = np.array([ 2.0000000e-03, 1.9069180e+00, -1.2827459e+00, -1.0278843e+00, 1.9599031e+00, -3.0928282e-01, -1.7105667e+00, 1.4564149e+00, 6.8071304e-01, -1.8799097e+00, 6.1372566e-01, 1.4111859e+00, -1.5603687e+00, -3.0450413e-01, 1.7222501e+00, -8.7117377e-01, -1.0823464e+00, 1.5749403e+00, -1.3548866e-02, -1.5065125e+00, 1.0101789e+00, 7.6873299e-01, -1.4585208e+00, 2.2623466e-01, 1.2166616e+00, -9.9519434e-01, -4.7264825e-01, 1.2053870e+00, -3.4519534e-01, -8.5756175e-01, 8.5622763e-01, 1.9940002e-01, -8.7211467e-01, 3.9123918e-01, 4.9281670e-01, -6.4777990e-01, 4.4045083e-03, 5.3522271e-01, -3.4587481e-01, -2.0867277e-01, 3.9176796e-01, -7.5815041e-02, -2.3235672e-01, 1.7076784e-01, 4.6201779e-02, -7.8306174e-02, -5.1144379e-03, -3.9879744e-02, 1.2998155e-01, 1.3860220e-02, -2.8030243e-01, 1.9759333e-01, 2.7404382e-01, -4.8857803e-01, 1.8222823e-03, 6.1879349e-01, -4.4649009e-01, -4.2028167e-01, 8.2464854e-01, -1.0688063e-01, -8.5761113e-01, 7.3211587e-01, 4.2190305e-01, -1.1056042e+00 ])
```

**Spectral analysis** is the task of analyzing a short snippet of signal, in order to find out what sinusoids it contains.

Spectral analysis is often used, for example, to study the chemical composition of some material sample (or of a distant supernova, or of a section of atmosphere). The way we do it is to heat the material, and then measure the electromagnetic signals it emits while it's heating. Any pure sinusoid, at any particular frequency, denotes the presence of some atom that has an oscillation at that frequency.

The big problem is that some of the sinusoids will be 20dB or even 60dB higher in amplitude than others. That makes it hard to determine exactly which sinusoids are present.

For example, consider the following signal. This signal is made of exactly 5 sinusoids:

$$x[n]=\sum_{q=1}^5 a_q \cos(\omega_q n+\theta_q)$$The task for today will be to try to identify those five frequencies, $\omega_1$ through $\omega_5$, from the 64-sample snippet of the signal that we have available.

In [2]:

```
import matplotlib.pyplot as plt
fig, axs = plt.subplots(figsize=(14,4))
axs.stem(signal)
len(signal)
```

Out[2]:

64

This is a 64-sample signal, with an obvious **beat frequency** of about 64 samples. So we can see that it probably is something like

Actually, this signal is the sum of **five sinusoids**, at five different frequencies, with five different amplitudes:

Some of the sinusoids have much lower amplitude, though, so we'll need to use clever application of windows in order to find all of their frequencies.

So far, we have learned the **discrete-time Fourier transform (DTFT)**. It looks like this:

Now we will learn the **discrete Fourier transform (DFT)**. It looks like this:

Notice that there are only two differences between the DFT and the DTFT:

**Finite length in time:**The DFT assumes that the input signal is only $N$ samples long.**Sampled in frequency:**Because the input only has $N$ independent numbers ($x[0]$ through $x[n_1]$), it's only really meaningful to calculate $N$ different samples in the frequency domain ($X[0]$ through $X[k]$). So, basically, the DFT is a sampled version of the DTFT:

The DFT is the only Fourier transform that is discrete in **both time and frequency**, so it is the only Fourier transform that we can compute on a computer.

We can simulate a DTFT on a computer, however, by **zero-padding** the input signal. This leaves a finite-length signal, but it makes $N$ much larger, so that we get a lot more samples in the frequency domain.

Here is a comparison between the 64-sample DFT of **signal**, and its 320-sample DFT. The 320-sample DFT can be considered similar to its DTFT. By plotting these two things together, we can see that the **DFT = samples of the DTFT in frequency**.

We will plot the magnitude DFT, and the log-magnitude DFT. The log-magnitude DFT allows us to learn more about low-amplitude samples.

In [3]:

```
N = len(signal)
print('N is ',N)
n_axis = np.arange(-2*N,3*N)
zero_padded_signal = np.concatenate((np.zeros(N),np.zeros(N),signal,np.zeros(N),np.zeros(N)))
fig,axs = plt.subplots(figsize=(14,4))
axs.plot(n_axis,zero_padded_signal)
axs.set_title('Signal, zero-padded to $5N=320$ samples in length')
```

N is 64

Out[3]:

Text(0.5, 1.0, 'Signal, zero-padded to $5N=320$ samples in length')

In [4]:

```
omega = np.linspace(0,2*np.pi,5*N,endpoint=False)
omega_k = np.linspace(0,2*np.pi, N, endpoint=False)
absDTFT = np.abs(np.fft.fft(zero_padded_signal))
absDFT = np.abs(np.fft.fft(signal))
fig,axs = plt.subplots(2,1,figsize=(14,8))
axs[0].plot(omega,absDTFT)
axs[0].stem(omega_k,absDFT)
axs[0].set_title('Magnitude DTFT, showing the DFT as samples in frequency')
axs[1].plot(omega,20*np.log10(np.maximum(0.001,absDTFT)))
axs[1].stem(omega_k, 20*np.log10(np.maximum(0.001,absDFT)))
axs[1].set_title('Log Magnitude DTFT, showing DFT as samples in frequency')
axs[1].set_ylabel('decibels')
axs[1].set_xlabel('Frequency (radians/sample)')
fig.tight_layout()
```

From this view, we can see clearly that the first two sinusoids are at exactly $k=19$ and $k=20$, i.e.,

$$\omega_1 = \frac{2\pi 19}{64},~~~\omega_2=\frac{2\pi 20}{64}$$The levels of these two tones are both about +30dB, so their coefficients are something like

$$a_1 = a_2 = 10^{30/20}$$But none of the other sinusoids are visible at all. All of the other sinusoids are masked by the sidelobes of the windowing function. In order to find the other sinusoids, we need to learn more about the sidelobes of the windowing function.

If we take the DTFT of a short signal, that's the same thing as taking the DTFT of a long signal multiplied by a rectangular window:

$$w_R[n] = \begin{cases}1&0\le n\le N-1\\0&\mbox{otherwise}\end{cases}$$In [5]:

```
N = len(signal)
w_R = np.ones(N)
w_R_zp = np.concatenate((np.zeros(N),np.zeros(N),np.ones(N),np.zeros(N),np.zeros(N)))
n_axis = np.arange(-2*N, 3*N)
fig, axs = plt.subplots(figsize=(14,4))
axs.stem(n_axis, w_R_zp)
axs.set_title('Rectangular Window of length $N=64$')
```

Out[5]:

Text(0.5, 1.0, 'Rectangular Window of length $N=64$')

We've already calculated the DTFT of a rectangular window a few times. It is

$$W_R(\omega) = \sum_{n=-\infty}^\infty w_R[n]e^{-j\omega n}=\sum_{n=0}^{N-1}e^{-j\omega n} = e^{-j\omega\frac{N-1}{2}} \frac{\sin(\omega N/2)}{\sin(\omega/2)}$$The real-valued part of that is called the "Dirichlet form." It's very much like a sinc function, except that it's periodic, with a period of $2\pi$.

$$D_N(\omega) = \frac{\sin(\omega N/2)}{\sin(\omega/2)}\approx \frac{\sin(\omega N/2)}{\omega/2}$$Among other things, we can see that

$$D_N(\omega) \approx \begin{cases} N&\omega=0\\\frac{N}{2}&\omega=\frac{\pi}{N}\\0&\omega=\frac{2\pi}{N}\\-\frac{2N}{3\pi}&\omega=\frac{3\pi}{N}\\0&\omega=\frac{4\pi}{N}\\\frac{2N}{5\pi}&\omega=\frac{5\pi}{N}\\\vdots & \vdots\end{cases}$$The sidelobes are not exactly $-2N/3\pi$, $2N/5\pi$, $-2N/7\pi$ and so on, but they're pretty close. Notice that $2N/3\pi$ is a really large number, for something that we would really prefer to be zero.

In [6]:

```
W_R_DTFT = np.abs(np.fft.fft(w_R_zp))
W_R_DFT = np.abs(np.fft.fft(w_R))
fig, axs = plt.subplots(2,1,figsize=(14,8))
axs[0].plot(omega, W_R_DTFT)
axs[0].stem(omega_k, W_R_DFT)
axs[0].set_title('Magnitude DFT of the rectangular window (64-point DFT shown as samples)')
axs[1].plot(omega, 20*np.log10(np.maximum(0.001,W_R_DTFT)))
axs[1].stem(omega_k, 20*np.log10(np.maximum(0.001, W_R_DFT)))
axs[1].set_title('Log Magnitude DTFT (decibels) of rectangular window (64-point DFT shown as samples)')
axs[1].set_ylabel('decibels')
axs[1].set_xlabel('Frequency (radians/sample)')
fig.tight_layout()
```