MP7: Spectral Analysis¶

In [1]:
import numpy as np
with open('signal.dat') as f:
    signal = np.array([ float(x) for x in f.read().split() ])

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
No description has been provided for this image

Contents¶

  1. Zero-padding
  2. Peak picking
  3. Analysis by synthesis
  4. Windows
  5. Windowing
  6. Checking your results

1. Zero-padding¶

The basis of all spectral estimation is the discrete Fourier transform (DFT). It looks like this:

$$X[k] = \sum_{n=0}^{N-1} x[n]e^{-j\frac{2\pi kn}{N}},~~~0\le k\le N-1$$

The DFT is a computable approximation of the discrete-time Fourier transform (DTFT), which looks like

$$X(\omega) = \sum_{n=-\infty}^\infty x[n]e^{-j\omega n},~~~0\le\omega < 2\pi$$

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

  1. Finite length in time: The DFT assumes that the input signal is only $N$ samples long.
  2. 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:

$$X[k] = X(\omega_k),~~\omega_k = \frac{2\pi k}{N}$$

The most obvious method of spectral analysis is to just compute and plot the DFT of our signal. Let's try it.

In [3]:
X = np.fft.fft(signal)
fig, axs = plt.subplots(2,figsize=(14,4),layout='tight')
axs[0].plot(np.abs(X))
axs[0].set_title('DFT versus bin number')
axs[1].plot(2*np.pi*np.arange(len(X))/len(X),20*np.log10(np.abs(X)))
axs[1].set_title('log DFT versus frequency in radians/sample')
Out[3]:
Text(0.5, 1.0, 'log DFT versus frequency in radians/sample')
No description has been provided for this image

As you can see, the DFT by itself is not terribly informative. We can see that $X[19]$ and $X[20]$ are about the same amplitude. That's consistent with at least two possible explanations:

  1. There might be a single sinusoid at a frequency of $k=19.5$, or
  2. ...there might be two separate sinusoids, one at frequency $k=19$, another at frequency $k=20$.

We can find out which answer is correct by zero-padding the signal. If the original signal length is $L=64$, we will zero-pad it by a factor of five to a length of $N=5\times 64=320$. Notice that this doesn't change the values of $X(\omega)$ at all, because $x[n]=0$ for $64\le n\le 319$:

$$X(\omega)=\sum_{n=0}^{63} x[n]e^{-j\omega n}=\sum_{n=0}^{319} x[n]e^{-j\omega n}$$

What it does accomplish is to compute a DFT that has five times as many frequency samples as the original DFT, so that we can see the same spectrum with $5\times$ better frequency resolution:

$$\omega_k = \frac{2\pi k}{N} = \frac{2\pi k}{5L}$$

In [4]:
import submitted, importlib
importlib.reload(submitted)
padded_signal, padded_dft = submitted.todo_zeropadded(signal, 5)
fig,axs = plt.subplots(figsize=(14,4))
axs.plot(padded_signal)
axs.set_title('Signal, zero-padded to $5N=320$ samples in length')
Out[4]:
Text(0.5, 1.0, 'Signal, zero-padded to $5N=320$ samples in length')
No description has been provided for this image
In [5]:
omega = np.linspace(0,2*np.pi,320,endpoint=False)
omega_k = np.linspace(0,2*np.pi, 64, endpoint=False)
fig,axs = plt.subplots(2,1,figsize=(14,8),layout='tight')
axs[0].plot(omega,np.abs(padded_dft))
axs[0].stem(omega_k,np.abs(padded_dft[::5]))
axs[0].set_title('Magnitude DTFT, showing the DFT as samples in frequency')
axs[1].plot(omega,20*np.log10(np.maximum(0.001,np.abs(padded_dft))))
axs[1].stem(omega_k, 20*np.log10(np.maximum(0.001,np.abs(padded_dft[::5]))))
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()
No description has been provided for this image

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

$$\omega_1 = \frac{2\pi 19}{64},~~~\omega_2=\frac{2\pi 20}{64}$$

  • Can we define some criterion that automatically picks out those two peaks from the spectrum?
  • None of the other sinusoids are visible at all, because zero-padding is like windowing some ideal infinite-length signal with a length-64 rectangular window. Are there any windows that are better than the rectangular window?

2. Peak Picking¶

In the DTFT, a convex peak is a frequency, $\omega$, such that

$$\frac{d|X(\omega)|}{d\omega}=0~\text{and}~\frac{d^2|X(\omega)|}{d\omega^2}>0$$

When we have sampled the DTFT at samples $\omega_k=\frac{2\pi k}{N}$, if we have enough frequency samples (because of zero-padding), we can define convex peaks pretty reasonably as

$$X(\omega) > X(\omega-\frac{2\pi}{N})~\text{and}~X(\omega)>X(\omega+\frac{2\pi}{N})$$

Let's write some code to find all such frequencies in padded_dft:

In [6]:
importlib.reload(submitted)
help(submitted.todo_peakpick)

Help on function todo_peakpick in module submitted:

todo_peakpick(dft, omegaL, omegaH, L)
    Return a list of (freq,ampl,phase) pairs of convex peaks in abs(dft) between omegaL and omegaH.

    @param:
    dft (N,complex) - the DFT of a zero-padded signal
    omegaL (1,) - consider only peaks with omega > omegaL (radians/sample)
    omegaH (1,) - consider only peaks with omega < omegaH (radians/sample)
    L (1,) - the length of the nonzero part of the signal (window length)

    @return:
    peaklist (list of tuples) - each tuple is (freq,ampl,phase):
     - freq is in radians/sample, to the nearest 2pi/N
     - ampl is the Fourier series amplitude, i.e., 2/L times DFT amplitude
     - phase is the Fourier series cosine phase, i.e., phase of the positive-frequency DFT peak
In [7]:
peaklist = submitted.todo_peakpick(padded_dft, 0.1*np.pi, 0.9*np.pi, 64)
print('There are',len(peaklist),'peaks, the first one is',peaklist[0])

There are 47 peaks, the first one is (0.3337942194439155, np.float64(0.038565613258077204), np.float64(0.30453610524365804))

Well, it's not very useful to know that there is such a small peak at $\omega=0.33$. Fortunately, we can use sorted in python to find the two peaks with the highest amplitude:

In [8]:
s = sorted(peaklist, key=lambda t: -t[1]) # sort in descending order of amplitude
print('The peak with highest amplitude is',s[0])
print('... its DFT amplitude is its cosine DFS amplitude times L/2, which is',s[0][1]*64/2)
print('The peak with second highest amplitude is',s[1])
print('... its DFT amplitude is its cosine DFS amplitude times L/2, which is',s[1][1]*64/2)

The peak with highest amplitude is (1.8456856839840037, np.float64(0.8087487634732259), np.float64(-0.4101587978863801))
... its DFT amplitude is its cosine DFS amplitude times L/2, which is 25.87996043114323
The peak with second highest amplitude is (1.9438604544086846, np.float64(0.803300652182659), np.float64(-1.7991520847976563))
... its DFT amplitude is its cosine DFS amplitude times L/2, which is 25.705620869845088

OK, that looks right. The frequencies 1.85 and 1.94 match the spectrum plot pretty well, as do the DFT amplitudes 25.9 and 25.7. Phase we haven't measured.

Is there any way to verify these numbers? Yes, indeed: We can check them using analysis-by-synthesis.

3. Analysis by Synthesis¶

Analysis is the process of extracting key numerical parameters from a signal. Synthesis is the process of generating a signal from its parameters. One of the most accurate analysis methods, if we have enough computation available, is called analysis-by-synthesis. Analysis-by-synthesis is an iterative approach: (1) Start with an estimate of the parameters, (2) Use those parameters to synthesize a signal, (3) Measure the error between the synthetic signal and the original signal, (4) Use the error to improve your estimate of the parameters.

In this case, our parameters suggest that the signal should be the sum of two cosines, multiplied by a rectangular window:

$$x[n]=a_1\cos(\omega_1 n+\theta_1)w[n]+a_2\cos(\omega_2 n+\theta_2)w[n]$$ $$w[n]=\left\{\begin{array}{ll} 1 & 0\le n\le 63\\0 & \mbox{otherwise}\end{array}\right.$$

There are two different ways to think about $x[n]$: (1) It is the sum of two windowed cosines, or (2) It is the sum of two modulated rectangles. Modulation is the process of multiplying any signal ($w[n]=$np.ones(64), in this case) by a complex exponential or a cosine. Here's the docstring:

In [9]:
importlib.reload(submitted)
help(submitted.todo_modulate)

Help on function todo_modulate in module submitted:

todo_modulate(signal, amplitude, frequency, phase)
    Modulate signal with a cosine at specified amplitude, frequency, and phase.

    @param:
    signal (L,) - the signal
    amplitude (1,) - the amplitude of the modulator cosine
    frequency (1,) - frequency, in radians/sample
    phase (1,) - phase of the modulator cosine

    @return:
    modulated (L,) - the modulated signal

OK, let's try it, and plot the result:

In [10]:
importlib.reload(submitted)
x1 = submitted.todo_modulate(np.ones(64), s[0][1], s[0][0], s[0][2])
x2 = submitted.todo_modulate(np.ones(64), s[1][1], s[1][0], s[1][2])
synthetic_padded, synthetic_dft = submitted.todo_zeropadded(x1+x2, 5)
fig,axs = plt.subplots(2,2,figsize=(14,8),layout='tight')
axs[0,0].stem(signal)
axs[0,0].set_title('Original Signal')
axs[0,1].plot(omega,20*np.log10(np.maximum(0.001,np.abs(padded_dft))))
axs[0,1].set_title('Spectrum of Original (dB)')
axs[1,0].stem(x1+x2)
axs[1,0].set_title('Synthetic Signal')
axs[1,1].plot(omega,20*np.log10(np.maximum(0.001,np.abs(synthetic_dft))))
axs[1,1].set_title('Spectrum of Synthetic (dB)')
Out[10]:
Text(0.5, 1.0, 'Spectrum of Synthetic (dB)')
No description has been provided for this image

That's pretty close, but not exactly the same! To get a better handle on the difference, let's try subtracting the two to find the error:

In [11]:
err = x1+x2 - signal
err_padded, err_dft = submitted.todo_zeropadded(err, 5)
fig,axs = plt.subplots(1,2,figsize=(14,4),layout='tight')
axs[0].stem(err)
axs[0].set_title('Error signal')
axs[1].plot(omega,20*np.log10(np.maximum(0.001,np.abs(err_dft))))
axs[1].set_title('Spectrum of Error (dB)')
Out[11]:
Text(0.5, 1.0, 'Spectrum of Error (dB)')
No description has been provided for this image

Oops. The error is complicated, and its spectrum is still completely dominated by the sidelobes of the rectangular window. It seems that, in order to find the components that are missing, we need to better understand the window spectrum, and possibly we should use a different window!

4. Windows¶

4.1 Rectangular Window¶

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 L-1\\0&\mbox{otherwise}\end{cases}$$

In [13]:
importlib.reload(submitted)
wR_padded, wR_dft = submitted.todo_zeropadded(np.ones(64), 5)
n_axis = np.arange(len(wR_padded))
fig, axs = plt.subplots(figsize=(14,4))
axs.stem(n_axis, wR_padded)
axs.set_title('Rectangular Window of length $N=64$')
Out[13]:
Text(0.5, 1.0, 'Rectangular Window of length $N=64$')
No description has been provided for this image

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}^{L-1}e^{-j\omega n} = e^{-j\omega\frac{L-1}{2}} \frac{\sin(\omega L/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_L(\omega) = \frac{\sin(\omega L/2)}{\sin(\omega/2)}\approx \frac{\sin(\omega L/2)}{\omega/2}$$

Among other things, we can see that

$$D_L(\omega) \approx \begin{cases} L&\omega=0\\\frac{L}{2}&\omega=\frac{\pi}{L}\\0&\omega=\frac{2\pi}{L}\\-\frac{2L}{3\pi}&\omega=\frac{3\pi}{L}\\0&\omega=\frac{4\pi}{L}\\\frac{2L}{5\pi}&\omega=\frac{5\pi}{L}\\\vdots & \vdots\end{cases}$$

The sidelobes are not exactly $-2L/3\pi$, $2L/5\pi$, $-2L/7\pi$ and so on, but they're pretty close. Notice that $2L/3\pi$ is a really large number, for something that we would really prefer to be zero. For example, if $L=64$, then the main lobe has an amplitude of $L=64$ ($20\log_10 64=36$dB), and the first sidelobe has an amplitude of $2L/3\pi=13.6$ ($20\log_10 13.6=23$dB), so the first sidelobe is lower than the mainlobe by only a factor of 0.21 (-13dB).

In [15]:
importlib.reload(submitted)
fig, axs = plt.subplots(2,1,figsize=(14,8), layout='tight')
axs[0].plot(omega, np.abs(wR_dft))
axs[0].stem(omega_k, np.abs(wR_dft[::5]))
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,np.abs(wR_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)')
Out[15]:
Text(0.5, 0, 'Frequency (radians/sample)')
No description has been provided for this image

To help us see the sidelobe amplitudes relative to the mainlobe amplitude, it's useful to divide the spectrum by $W(0)$:

In [16]:
importlib.reload(submitted)
fig, axs = plt.subplots(2,1,figsize=(14,8), layout='tight')
axs[0].plot(omega, np.abs(wR_dft/wR_dft[0]))
axs[0].stem(omega_k, np.abs(wR_dft[::5]/wR_dft[0]))
axs[0].set_title('Magnitude DFT of the rectangular window, normalized')
axs[1].plot(omega, 20*np.log10(np.maximum(0.001,np.abs(wR_dft/wR_dft[0]))))
axs[1].set_title('Log Magnitude DTFT (decibels) of rectangular window, normalized')
axs[1].set_ylabel('decibels')
axs[1].set_xlabel('Frequency (radians/sample)')
Out[16]:
Text(0.5, 0, 'Frequency (radians/sample)')
No description has been provided for this image

We've already seen how the sidelobes of a rectangular window make it really, really hard to identify low-amplitude sinusoidal components in your input signal. The first sidelobe is down by only 0.21 (-13dB), and even the distant sidelobes are only down by about 0.03 (-30dB).

4.2. Hamming window¶

The Hamming window was designed with the following idea:

$$w_H[n] = \begin{cases}a + b\cos\left(\frac{2\pi n}{N-1}\right)&0\le n\le N-1\\0&\mbox{otherwise}\end{cases}$$

$$a,b = \arg\min \left|W_H\left(\omega=\frac{3\pi}{N}\right)\right|$$

In words, the coefficients in the Hamming window are chosen in order to minimize the amplitude of the first sidelobe. It would be possible to choose $a$ to set the first sidelobe to exactly zero, but by convention, we round it off to only two significant figures ($a=0.54$), so the first sidelobe is only guaranteed to have an amplitude less than 0.01. The resulting coefficients are

$$w_H[n] = 0.54 - 0.46\cos\left(\frac{2\pi n}{N-1}\right)$$

The main lobe now has a half-width of $4\pi/N$, and the first sidelobe is at $5\pi/N$, and the resulting sidelobe amplitude is

$$\left|W_H\left(\frac{5\pi}{N}\right)\right| = 0.006,~~~20\log_{10}\left|W_H\left(\frac{5\pi}{N}\right)\right|=-44\mbox{dB}$$

In [17]:
importlib.reload(submitted)
wH_padded, wH_dft = submitted.todo_zeropadded(np.hamming(64), 5)
fig,axs = plt.subplots(figsize=(14,4))
axs.stem(n_axis, wH_padded)
axs.set_title('Hamming window, zero-padded to a length of $5N$')
Out[17]:
Text(0.5, 1.0, 'Hamming window, zero-padded to a length of $5N$')
No description has been provided for this image
In [18]:
fig, axs = plt.subplots(2,1,figsize=(14,8))
axs[0].plot(omega, np.abs(wH_dft/wH_dft[0]))
axs[0].stem(omega_k, np.abs(wH_dft[::5]/wH_dft[0]))
axs[0].set_title('Magnitude DFT of the Hamming window, normalized')
axs[1].plot(omega, 20*np.log10(np.maximum(0.001,np.abs(wH_dft/wH_dft[0]))))
axs[1].set_title('Log Magnitude DTFT (decibels) of Hamming window, normalized')
axs[1].set_ylabel('decibels')
axs[1].set_xlabel('Frequency (radians/sample)')
fig.tight_layout()
No description has been provided for this image

Notice that the first sidelobe is now down by 44dB! We pay two prices for this very low sidelobe. First, the main lobe is twice as wide. Second, by minimizing the first sidelobe, we didn't really minimize the rest of the sidelobes.

4.3 Hann Window¶

A Hann window is like a Hamming window, except that it's made to taper all the way down to zero at either end:

$$w_{N}[n]=\begin{cases}0.5-0.5\cos\left(\frac{2\pi n}{L+1}\right)&0\le n\le L-1\\0&\mbox{otherwise}\end{cases}$$

Remember that the Hamming coefficients $(a,b)=(0.54,-0.46)$ were optimized for the smallest possible first sidelobe. The Hann coefficients $(a,b)=(0.5,-0.5)$ will therefore have a larger first sidelobe, but some of the higher-frequency sidelobes will be much smaller.

In particular:

  • The Hamming window's first sidelobe is at -44dB, and all of its other sidelobes are at about the same level.

  • The Hann window's first sidelobe is at only -30dB, but its distant sidelobes decay to less than -80dB.

In [19]:
importlib.reload(submitted)
wN_padded, wN_dft = submitted.todo_zeropadded(np.hanning(64), 5)
fig,axs = plt.subplots(figsize=(14,4))
axs.stem(n_axis, wN_padded)
axs.set_title('Hann window, zero-padded to a length of $5N$.  Notice how it tapers all the way down to zero.')
Out[19]:
Text(0.5, 1.0, 'Hann window, zero-padded to a length of $5N$.  Notice how it tapers all the way down to zero.')
No description has been provided for this image
In [20]:
fig, axs = plt.subplots(2,1,figsize=(14,8))
axs[0].plot(omega, np.abs(wN_dft/wN_dft[0]))
axs[0].stem(omega_k, np.abs(wN_dft[::5]/wN_dft[0]))
axs[0].set_title('Magnitude DFT of the Hann window, normalized')
axs[1].plot(omega, 20*np.log10(np.maximum(0.00001,np.abs(wN_dft/wN_dft[0]))))
axs[1].set_title('Log Magnitude DTFT (decibels) of Hann window, normalized')
axs[1].set_ylabel('decibels')
axs[1].set_xlabel('Frequency (radians/sample)')
fig.tight_layout()
No description has been provided for this image

4.4 Bartlett (Triangular) Window¶

A Bartlett window, or triangle window, is just what its name says: a triangle. For example, if $L$ is an even number, the Bartlett window can be defined as

$$w_B[n]=\begin{cases}\min\left(\frac{2n+1}{L},\frac{2(L-n)-1}{L}\right)&0\le n\le L-1\\0&\mbox{otherwise}\end{cases}$$

Notice that a Bartlett window is the convolution of a rectangular window with itself:

$$w_{B,L}[n] = w_{R,L/2}[n] \ast w_{R,L/2}[n]$$

Therefore its DTFT is the product:

$$W_{B,L}(\omega) = W_{R,L/2}(\omega)^2$$

In particular, its sidelobes happen at frequencies like $\frac{6\pi}{L},\frac{10\pi}{L},\ldots$ instead of frequencies like $\frac{3\pi}{L},\frac{5\pi}{L},\ldots$. Because of this, windowing a signal using a Bartlett window can sometimes avoid masking components that would be masked by any other window.

In [21]:
importlib.reload(submitted)
wB_padded, wB_dft = submitted.todo_zeropadded(np.bartlett(64), 5)
fig,axs = plt.subplots(figsize=(14,4))
axs.stem(n_axis, wB_padded)
axs.set_title('Bartlett window, zero-padded to a length of $5N$.')
Out[21]:
Text(0.5, 1.0, 'Bartlett window, zero-padded to a length of $5N$.')
No description has been provided for this image
In [22]:
fig, axs = plt.subplots(2,1,figsize=(14,8))
axs[0].plot(omega, np.abs(wB_dft/wB_dft[0]))
axs[0].stem(omega_k, np.abs(wB_dft[::5]/wB_dft[0]))
axs[0].set_title('Magnitude DFT of the Bartlett window (64-point DFT shown as samples)')
axs[1].plot(omega, 20*np.log10(np.maximum(0.0001,np.abs(wB_dft/wB_dft[0]))))
axs[1].set_title('Log Magnitude DTFT (decibels) of Bartlett window.')
axs[1].set_ylabel('decibels')
axs[1].set_xlabel('Frequency (radians/sample)')
fig.tight_layout()
No description has been provided for this image

5. Windowing¶

When we multiply two signals in the time domain, that's the same as convolving them in the frequency domain:

$$x[n]=w[n]s[n] \Leftrightarrow X(\omega)=\frac{1}{2\pi} W(\omega)\ast S(\omega)$$

Convolution tends to blur signal components. In particular, if $s[n]$ is a sum of cosines, then $S(\omega)$ is a sum of sharp impulses --- but $X(\omega)$ is the sum of shifted window spectra.

$$s[n]=a_1\cos(\omega_1 n)+a_2\cos(\omega_2 n)\Leftrightarrow S(\omega)=a_1\pi\left(\delta(\omega-\omega_1)+\delta(\omega+\omega_1)\right)+a_2\pi\left(\delta(\omega-\omega_2)+\delta(\omega-\omega_2)\right)$$

$$x[n]=a_1\cos(\omega_1 n)w[n]+a_2\cos(\omega_2 n)w[n]\Leftrightarrow S(\omega)=\frac{a_1}{2}\left(W(\omega-\omega_1)+W(\omega+\omega_1)\right)+\frac{a_2}{2}\left(W(\omega-\omega_2)+W(\omega-\omega_2)\right)$$

If $a_2$ is much, much less than $a_1$, it might be completely masked by the sidelobes of $W(\omega-\omega_1)$. In order to avoid that, we can multiply by some window other than the rectangular window:

  • If $\omega_2$ is near the first sidelobe of $W(\omega-\omega_1)$, try using a Hamming window
  • If $\omega_2$ is far away from $\omega_1$ but $a_2$ is very, very small, try using a Hann or Bartlett window
In [23]:
importlib.reload(submitted)
help(submitted.todo_applywindow)

Help on function todo_applywindow in module submitted:

todo_applywindow(signal, window)
    Window the signal.

    @param:
    signal (L,) - the signal to be windowed
    window (L,) - the window to be applied

    @return:
    windowed (L,) - the windowed signal

5.1 Rectangular Window¶

We've already computed the spectrum with a rectangular window. Let's try something else.

5.2 Hamming window¶

In [24]:
importlib.reload(submitted)
signal_hamming = submitted.todo_applywindow(signal, np.hamming(64))
signal_hamming_zp, signal_hamming_dft = submitted.todo_zeropadded(signal_hamming, 5)
fig,axs = plt.subplots(2,1,figsize=(14,8))
axs[0].plot(n_axis,padded_signal)
axs[0].set_title('Rectangular-windowed signal, zero-padded to a length of $5N$')
axs[1].plot(n_axis, signal_hamming_zp)
axs[1].set_title('Hamming windowed signal, zero-padded to a length of $5N$')
axs[1].set_xlabel('Time (samples)')
fig.tight_layout()
No description has been provided for this image
In [62]:
fig,axs = plt.subplots(2,1,figsize=(14,8))
axs[0].plot(omega,np.abs(signal_hamming_dft))
axs[0].stem(omega_k,np.abs(signal_hamming_dft[::5]))
axs[0].set_title('Magnitude DTFT of Hamming-windowed signal')
axs[1].plot(omega,20*np.log10(np.maximum(0.001,np.abs(signal_hamming_dft))))
axs[1].stem(omega_k, 20*np.log10(np.maximum(0.001,np.abs(signal_hamming_dft[::5]))))
axs[1].set_title('Log Magnitude DTFT of Hamming-windowed signal')
axs[1].set_ylabel('decibels')
axs[1].set_xlabel('Frequency (radians/sample)')
fig.tight_layout()
No description has been provided for this image

In this plot, we can see one more frequency component that wasn't visible using the rectangular window! It is right in between $k=16$ and $k=17$, so we could say that it's at frequency $k=16.5$.

$$\omega_1=\frac{2\pi 19}{64},~~~\omega_2=\frac{2\pi 20}{64},~~~\omega_3=\frac{2\pi 16.5}{64}$$

with amplitudes of

$$a_1=a_2=10^{30/20},~~~a_3 = 10^{-5/20}$$

Notice that the Hamming window has a very low first sidelobe, but the disadvantage of the Hamming window is that its main lobe is twice as wide. In this figure, for example, we can't see the separation between the $\omega_1$ and $\omega_2$; they look like one large peak. We need the rectangular window to distinguish $\omega_1$ and $\omega_2$, but we need the Hamming window to detect $\omega_3$.

5.3. Hann window¶

In [63]:
importlib.reload(submitted)
signal_hann = submitted.todo_applywindow(signal, np.hanning(64))
signal_hann_zp, signal_hann_dft = submitted.todo_zeropadded(signal_hann, 5)
fig,axs = plt.subplots(2,1,figsize=(14,8))
axs[0].plot(n_axis,padded_signal)
axs[0].set_title('Rectangular-windowed signal, zero-padded to a length of $5N$')
axs[1].plot(n_axis, signal_hann_zp)
axs[1].set_title('Hann windowed signal, zero-padded to a length of $5N$')
axs[1].set_xlabel('Time (samples)')
fig.tight_layout()
No description has been provided for this image
In [64]:
fig,axs = plt.subplots(2,1,figsize=(14,8))
axs[0].plot(omega,np.abs(signal_hann_dft))
axs[0].stem(omega_k,np.abs(signal_hann_dft[::5]))
axs[0].set_title('Magnitude DTFT of Hann-windowed signal')
axs[1].plot(omega,20*np.log10(np.maximum(0.001,np.abs(signal_hann_dft))))
axs[1].stem(omega_k, 20*np.log10(np.maximum(0.001,np.abs(signal_hann_dft[::5]))))
axs[1].set_title('Log Magnitude DTFT of Hann-windowed signal')
axs[1].set_ylabel('decibels')
axs[1].set_xlabel('Frequency (radians/sample)')
fig.tight_layout()
No description has been provided for this image

There's two more tones! There is certainly a tone at $k=29.5$, with a level of about -30dB. There also seems to be another tone at about $k=6.5$, with a level of about -40dB.

$$\omega_1=\frac{2\pi 19}{64},~~\omega_2=\frac{2\pi 20}{64},~~\omega_3=\frac{2\pi 16.5}{64},~~\omega_4=\frac{2\pi 29.5}{64},~~\omega_5=\frac{2\pi 6.5}{64}$$

with amplitudes of

$$a_1=a_2=10^{30/20},~~a_3 = 10^{-5/20},~~a_4=10^{-25/40},~~a_5=10^{-40/20}$$

5.4. Bartlett Window¶

In [65]:
importlib.reload(submitted)
signal_bartlett = submitted.todo_applywindow(signal, np.bartlett(64))
signal_bartlett_zp, signal_bartlett_dft = submitted.todo_zeropadded(signal_bartlett, 5)
fig,axs = plt.subplots(2,1,figsize=(14,8))
axs[0].plot(n_axis,padded_signal)
axs[0].set_title('Rectangular-windowed signal, zero-padded to a length of $5N$')
axs[1].plot(n_axis, signal_bartlett_zp)
axs[1].set_title('Bartlett windowed signal, zero-padded to a length of $5N$')
axs[1].set_xlabel('Time (samples)')
fig.tight_layout()
No description has been provided for this image
In [66]:
fig,axs = plt.subplots(2,1,figsize=(14,8))
axs[0].plot(omega,np.abs(signal_bartlett_dft))
axs[0].stem(omega_k,np.abs(signal_bartlett_dft[::5]))
axs[0].set_title('Magnitude DTFT of Bartlett-windowed signal')
axs[1].plot(omega,20*np.log10(np.maximum(0.001,np.abs(signal_bartlett_dft))))
axs[1].stem(omega_k, 20*np.log10(np.maximum(0.001,np.abs(signal_bartlett_dft[::5]))))
axs[1].set_title('Log Magnitude DTFT of Bartlett-windowed signal.  Notice how each sidelobe is split by the double peaks at $2\pi 19/64$ and $2\pi 20/64$')
axs[1].set_ylabel('decibels')
axs[1].set_xlabel('Frequency (radians/sample)')
fig.tight_layout()
No description has been provided for this image

Well, no new signal components obvious in this window -- but we've already found five, so let's call it done.

6. Checking your results¶

Before you submit your code to the autograder, you can check your results using grade.py:

In [73]:
!python grade.py

....
----------------------------------------------------------------------
Ran 4 tests in 0.015s

OK
In [ ]: