Lab 4: Frequency Response and Sampling

Lab 4: Frequency Response and Sampling

This lab will cover the frequency response of LSI systems, the frequency content of digital signals, and sampling basics. We have some interesting applications to get to, so let's get started!

Discrete Time Fourier Transform and Frequency Response
We will begin with a brief overview of the Discrete Time Fourier Transform (DTFT) and the frequency response of LSI systems.

The DTFT is the discrete-time version of our continuous-time Fourier transform (CTFT) from ECE 210. Like the CTFT, the DTFT is a complex-valued function that allows us to examine the frequency content of a signal or system. The DTFT is defined by

$$ X(\omega) = \sum_{n=-\infty}^{\infty}x[n]e^{-j\omega n} $$
We also must remember that the DTFT is completely represented by digital frequencies $-\pi$ to $\pi$ and is $2\pi$ periodic. Be careful labeling your frequency axis when taking the DTFT. In Python, we are able to take the DTFT of a signal using numpy's fft module. The two main functions we will use from this module are $\textrm{numpy.fft.fft()}$ and $\textrm{numpy.fft.rfft()}$. The "fft" in these functions is the Fast Fourier Transform, which is a computationally efficient way of computing the Discrete Fourier Transform of digital signals. You will learn more about the DFT and FFT in ECE 310 and Lab 5 of this course, but for this lab just think of it as a way of computing the DTFT of a digital signal.

Let's look at example usage for these two functions and what makes them different.

#import libraries first
import numpy as np
import matplotlib.pyplot as plt

from IPython.display import Audio
from scipy import signal
from scipy.io import wavfile
from skimage.io import imread
from sklearn.cluster import KMeans

%matplotlib inline
x = [0,1,0,2,0,2,0,1,0] #test signal

full_fft = np.fft.fft(x)
real_fft = np.fft.rfft(x)

omega_full = np.linspace(0,2*np.pi,len(full_fft)) #left limit, right limit, # pts
omega_real = np.linspace(0,np.pi,len(real_fft))

plt.figure(figsize=(15,6))
plt.subplot(121)
plt.title('Full FFT')
plt.xlabel('$\omega$')
plt.ylabel('Magnitude Response')
plt.plot(omega_full,np.absolute(full_fft))

plt.subplot(122)
plt.title('Real FFT')
plt.xlabel('$\omega$')
plt.ylabel('Magnitude Response')
plt.plot(omega_real,np.absolute(real_fft))
[<matplotlib.lines.Line2D at 0x7ff1696585c0>]

Observe the differences between the $\textrm{fft()}$ and $\textrm{rfft()}$ results. The way we have created the frequency axis points may have spoiled the answer, but we see that the $\textrm{fft()}$ function returns a DTFT with frequencies from 0 to $2\pi$ while $\textrm{rfft()}$ just gives us 0 to $\pi$: the real frequencies. It is important to acknowledge when the real frequencies are sufficient. If our signal is real-valued, we know that our spectrum will be Hermitian symmetric. In other words:

$$ x[n]\textrm{ real}\implies X(\omega) = X^*(-\omega), $$
where $X^*$ refers to the complex conjugate of the DTFT. Why is this important? Well, if our spectrum is Hermitian symmetric, then the spectrum's magnitude response is even symmetric and its phase response is odd symmetric:

$$ x[n]\textrm{ real}\implies |X(\omega)| = |X(-\omega)|\textrm{ and }\angle X(\omega) = -\angle X(-\omega). $$
Thus, if we want to look at the magnitude spectrum of a real-valued signal, it is sufficent to just look at the 0 to $\pi$ interval since it contains all unique information about the frequency content of our signal.

Now, let's make our plots look nicer too. We currently have a couple issues with them. First, they are very low resolution and coarse. Second, the full DTFT example is not zero-centered. Both problems can easily be fixed as follows:

x = [0,1,0,2,0,2,0,1,0]

full_fft = np.fft.fft(x,512)
centered_fft = np.fft.fftshift(full_fft) #shifts central frequency to middle of array
real_fft = np.fft.rfft(x,512)

omega_full = np.linspace(-np.pi,np.pi,len(centered_fft)) #new frequency axis
omega_real = np.linspace(0,np.pi,len(real_fft))

plt.figure(figsize=(15,6))
plt.subplot(121)
plt.title('Full FFT')
plt.xlabel('$\omega$')
plt.ylabel('Real Part')
plt.plot(omega_full,np.absolute(centered_fft))

plt.subplot(122)
plt.title('Real FFT')
plt.xlabel('$\omega$')
plt.ylabel('Real Part')
plt.plot(omega_real,np.absolute(real_fft))
[<matplotlib.lines.Line2D at 0x7ff169874908>]

Those plots look so much better now! To fix our first problem with the resolution, we pass a second argument to the $\textrm{fft()}$ and $\textrm{rfft()}$ functions to specify the number of points. Do not worry about the underlying math for now, it will be covered later in ECE 310. For now, think of the number of points as dictating how many frequencies we would like to use in capturing the signal's DTFT. Try changing the number of points in the above code and observe the resolution of the DTFT. We fix the second problem by using $\textrm{np.fft.fftshift()}$ on our full DTFT result. This function zero-centers our DTFT for us along the frequency axis: how convenient!

For the rest of this lab and in the future, it is critical you remember these tips. Always make sure your FFT has a enough points to look clean, and always make sure to appropriately label your frequency axis. Also, for this lab, when we say to "take the DTFT of a signal", use the $\textrm{np.fft.rfft()}$ function unless noted otherwise since we will only work with real signals.

Lastly, let's see how we can look at the frequency response of an LSI system. The function we will use is $\textrm{signal.freqz()}$, which returns the normalized digital frequencies and frequency response given the numerator and denominator coefficients for an LSI system's transfer function. It is convention to plot a system's frequency response on a dB scale $\left(20\cdot\log_{10}(x)\right)$.

b = [np.sin((np.pi/2)*n)/(0.5*np.pi*n) if n != 0 else 1 for n in range(-100,101)] #numerator coefficients
a = [1,0] #denominator coefficients
w,h = signal.freqz(b,a) #w = omega/digital frequencies, h = frequency response
plt.figure(figsize=(10,6))
plt.title('Toy Frequency Response')
plt.plot(w,20*np.log10(np.absolute(h))) #plot magnitude of frequency response with db-scaling on y-axis
plt.xlabel('$\omega$')
plt.ylabel('Magnitude Response (dB)')
Text(0, 0.5, 'Magnitude Response (dB)')

Exercise 1: Implementing the DTFT
We will begin by implementing the DTFT according to the above definition for an arbitrary collection of frequencies. We will test using signals over a finite support, so we will modify our definition of the DTFT to simply say

$$ X(\omega) = \sum_{n=0}^{N-1}x[n]e^{-j\omega n}. $$
a. Implement the $\textrm{myDTFT()}$ function below, which returns the DTFT values for a given list of digital frequencies.

b. Use your DTFT function to compute the DTFT of $x[n] = \cos(\frac{\pi}{2}n)$, $0\leq n < 50$ for 50 evenly spaced frequencies from $-\pi$ to $\pi$ (non-inclusive). Also, compute the DTFT of $x[n]$ using $\textrm{np.fft.fft()}$. Plot the magnitude and phase of the two implementations in separate figures to verify you achieve the same result. Use $\textrm{np.absolute()}$ and $\textrm{np.angle()}$ for the magnitude and phase responses, respectively. Don't forget to zero-center the frequency axis of the $\textrm{np.fft.fft()}$ result using $\textrm{np.fft.fftshift()}$.

c. Theoretically, the DTFT of $\cos(\omega_0n)$ should give us Kronecker deltas at $\pm\omega_0$. However, we see our implementation and the numpy DTFT function result in some non-ideal representation, like in the ramping behavior around the frequencies of the cosine. Why does this happen? Hint: consider how our practical definition of the DTFT in this exercise differse from the theoretical definition of the DTFT.

#Function to implement for part 1.a:
"""
Inputs:
x - input signal (list or np.array)
w - frequencies we want to compute the DTFT for (list or np.array)
Output:
dtft - value of the DTFT for signal x at each frequency specified in w (list or np.array)
"""
def myDTFT(x,w):
    dtft = [] #create empty list to append resulting computation
    #Iterate over each frequency in w
    for i in range (len(w)):
        DTFT_sum = []
        for j in range (len(x)):
            DTFT_sum.append(x[j]*np.exp(-1j*w[i]*j))
        
    #Compute summation for current frequency according to our DTFT definition, "1j" gives you the imaginary number
    
        dtft.append(sum(DTFT_sum))
        
    #Append result for current frequency
    
    return dtft

#Code for part 1.b:
x = np.array([np.cos(0.5*np.pi*n) for n in range(50)])
DTFT_xn = myDTFT(x,w)

#Code for magnitude and frequency response using the myDTFT function

full_fft = np.fft.fft(x,50)
centered_fft1 = np.fft.fftshift(DTFT_xn)
omega_full = np.linspace(-np.pi,np.pi,len(centered_fft1), endpoint = False)
plt.figure(figsize=(15,6))
plt.subplot(221)
plt.title('Magnitude Response using myDTFT function ')
plt.xlabel('$\omega$')
plt.ylabel('Real Part Magnitude')
plt.plot(omega_full,np.absolute(centered_fft1))
plt.subplot(222)
plt.title('Phase Response using myDTFt function ')
plt.xlabel('$\omega$')
plt.ylabel('Phase')
plt.plot(omega_full,np.angle(centered_fft1))

#Code for magnitude and frequency response using the numpy function for fast fourier transform


full_fft2 = np.fft.fft(x,50)
centered_fft2 = np.fft.fftshift(full_fft2)
omega_full = np.linspace(-np.pi,np.pi,len(centered_fft2))
plt.subplot(223)
plt.title('Magnitude Response using numpy function ')
plt.xlabel('$\omega$')
plt.ylabel('Real Part')
plt.plot(omega_full,np.absolute(centered_fft2))
plt.subplot(224)
plt.title('Phase Response using numpy function ')
plt.xlabel('$\omega$')
plt.ylabel('Phase')
plt.plot(omega_full,np.angle(centered_fft2))

#endpoint argument makes sure our definition aligns with the np.fft.fft result
w = np.linspace(-np.pi,np.pi,50,endpoint=False)

Answer for part 1.c:

One of the reasons why the representation is not ideal is because the sample length is finite whereas the theoretical definition is an infinite sum. This therefore does not give us the complete picture of the DTFT of the cosine function since our values for the eigenvector (exp^-jwn) will only sum through a finite amount of these eigenvectors, and not infinitely many n. Another reason why the sample is not depicted with the expected Kronecker deltas is because since we are trying to emulate a theoretical function (delta function) onto a real valued plot, it will not show us the actual function and therefore will cause the phase response plot to have ramping behavior near the (+-) omega values.

Exercise 2: Toy LSI Signals and Systems
For this exercise, we will get some practice with inspecting the DTFT of digital signals and frequency response of LSI systems. For parts 2.a and 2.b, plot the resulting magnitude and phase of the DTFT of each signal. Please place the two plots for each signal side-by-side using $\textrm{plt.subplot()}$. Don't forget to specify a larger number of points in your DTFT like 512. For parts 1.c and 1.d, plot the magnitude response of each system as shown above using $\textrm{signal.freqz()}$ use a dB scaling on the y-axis.

$\begin{align} a. x_1[n] = \frac{1}{4}\delta[n]+\frac{1}{2}\delta[n-1] + \delta[n-2] + \frac{1}{2}\delta[n-3] + \frac{1}{4}\delta[n-4] \end{align}$

$\begin{align} b. x_2[n] = -\delta[n] + 2\delta[n-2] - \delta[n-4] \end{align}$

$\begin{align} c. H_3(z) = \frac{z^2 -2z + 1}{z^2 -\frac{1}{2}z + \frac{1}{4}} \end{align}$

$\begin{align} d. H_4(z) = \frac{z^4 + 2z^3 + z^2}{z^4-\frac{1}{2}z^3+\frac{1}{4}z^2-\frac{1}{8}z + \frac{1}{16}} \end{align}$

#Code for 2.a:
#Remember to plot magnitude and phase of signals side-by-side with plt.subplot(nrows,ncols,plot number)
x1 = [1./4,1./2,1,1./2,1./4]

real_fft = np.fft.rfft(x1,512)
omega_real = np.linspace(0,np.pi,len(real_fft))

plt.figure(figsize=(20,20))

plt.subplot(321)
plt.title('Magnitude Response of X1')
plt.xlabel('$\omega$')
plt.ylabel('Real Part')
plt.plot(omega_real,np.absolute(real_fft))
plt.subplot(322)
plt.title('Phase Response of X1 ')
plt.xlabel('$\omega$')
plt.ylabel('Phase Response')
plt.plot(omega_real,np.angle(real_fft))
#Code for 2.b:
x2 = [-1,0,2,0,-1]

real_fft = np.fft.rfft(x2,512)
omega_real = np.linspace(0,np.pi,len(real_fft))

plt.subplot(323)
plt.title('Magnitude Response of X2')
plt.xlabel('$\omega$')
plt.ylabel('Real Part')
plt.plot(omega_real,np.absolute(real_fft))
plt.subplot(324)
plt.title('Phase Response of X2 ')
plt.xlabel('$\omega$')
plt.ylabel('Phase Response')
plt.plot(omega_real,np.angle(real_fft))

#Code for 2.c:
#Remember to plot magnitude response with dB-scaling on y-axis.
b = [1,-2,1]
a = [1,-1./2,1./4]

w,h = signal.freqz(b,a) #w = omega/digital frequencies, h = frequency response
plt.subplot(325)
plt.title('Magnitude Response of H3')
plt.plot(w,20*np.log10(np.absolute(h))) #plot magnitude of frequency response with db-scaling on y-axis
plt.xlabel('$\omega$')
plt.ylabel('Magnitude Response (dB)')
#Code for 2.d:
b = [1,2,1,0,0]
a = [1,-1./2,1./4,-1./8,1./16]

w,h = signal.freqz(b,a) #w = omega/digital frequencies, h = frequency response
plt.subplot(326)
plt.title('Magnitude Response of H4')
plt.plot(w,20*np.log10(np.absolute(h))) #plot magnitude of frequency response with db-scaling on y-axis
plt.xlabel('$\omega$')
plt.ylabel('Magnitude Response (dB)')
/Users/ryanyoseph/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:45: RuntimeWarning: divide by zero encountered in log10
Text(0, 0.5, 'Magnitude Response (dB)')

Sinusoidal Response of LSI Systems
We recall from ECE 210 that the response of LTI systems to sinusoidal inputs has a nice closed form using the frequency response of our system. The same is true for our discretized LSI systems! Let $H(\omega)$ be our frequency response and our input be some sinusoidal input with arbitrary amplitude $A$, frequency $\omega_0$, and phase $\theta$:

$$ A\sin(\omega_0n + \theta) \to \begin{bmatrix}H(\omega)\end{bmatrix} \to A|H(\omega_0)|sin(\omega_0n + \theta + \angle H(\omega_0)) $$
We see that the output is simply the input signal scaled by the magnitude response and we add phase according the value of the phase response at the sinusoid's frequency. Furthermore, we can extend this notion to sums of sinusoids by the linearity of our LSI systems. Now, let's verify this with an example system!

Exercise 3: Sinusoidal Response of an LSI System
For all parts of this problem, we will consider the following LSI system:

$$ h[n] = \frac{1}{2}\delta[n] + \delta[n-2] + \frac{1}{2}\delta[n-4] $$
a. Plot the magnitude and phase response of the above system. Do not use a dB-scale when plotting this system. Use 512 for the number of points in your DTFT. Verify your results by computing the frequency response by hand.

Consider the following two inputs:

$\begin{align} \bullet~x_b[n] = 1 + 2\sin\left(\frac{\pi}{4}n\right), \quad 0\leq n < 100 \end{align}$

$\begin{align} \bullet~x_c[n] = 2 + 10\sin\left(\frac{\pi}{2}n\right) + 4\cos\left(\pi n\right), \quad 0\leq n < 100 \end{align}$

b. Apply the filter $h[n]$ to input $x_b[n]$ using $\textrm{signal.lfilter()}$. Plot the magnitude of the DTFT for the input and filtered output, respectively, on separate subplots. Use 512 points again for your DTFT. Verify these results by hand!

c. Apply the filter $h[n]$ to input $x_c[n]$. Plot the magnitude of the DTFT for the input and filtered output, respectively, on separate subplots. Also verify these results by hand!

#Code for part 3.a:
h_n = [1./2,0,1,0,1./2]

plt.figure(figsize=(20,20))

real_fft = np.fft.rfft(h_n,512)
omega_real = np.linspace(0,np.pi,len(real_fft))

plt.subplot(321)
plt.title('Magnitude Response of H_n')
plt.xlabel('$\omega$')
plt.ylabel('Real Part')
plt.plot(omega_real,np.absolute(real_fft))
plt.subplot(322)
plt.title('Phase Response of H_n')
plt.xlabel('$\omega$')
plt.ylabel('Phase Response')
plt.plot(omega_real,np.angle(real_fft))

#Code for part 3.b:
x_b = np.array([1 + 2*np.sin(0.25*np.pi*n) for n in range(100)])

b = [1,0,2,0,1]
a = [2,0,0,0,0]

y1 = signal.lfilter(b,a,x_b)

real_fft_input = np.fft.rfft(x_b,512)
omega_real = np.linspace(0,np.pi,len(real_fft_input))
real_fft_filtered = np.fft.rfft(y1,512)


plt.subplot(323)
plt.title('Magnitude Response of X_b')
plt.xlabel('$\omega$')
plt.ylabel('Real Part')
plt.plot(omega_real,np.absolute(real_fft_input))
plt.subplot(324)
plt.title('Magnitude Response of Y1')
plt.xlabel('$\omega$')
plt.ylabel('Real Part')
plt.plot(omega_real,np.absolute(real_fft_filtered))

#Code for part 3.c:
x_c = np.array([2+ 10*np.sin(0.5*np.pi*n) + 4*np.cos(np.pi*n) for n in range(100)])

b = [1,0,2,0,1]
a = [2,0,0,0,0]

y2 = signal.lfilter(b,a,x_c)

real_fft_input = np.fft.rfft(x_c,512)
omega_real = np.linspace(0,np.pi,len(real_fft_input))
real_fft_filtered = np.fft.rfft(y2,512)


plt.subplot(325)
plt.title('Magnitude Response of X_c')
plt.xlabel('$\omega$')
plt.ylabel('Real Part')
plt.plot(omega_real,np.absolute(real_fft_input))
plt.subplot(326)
plt.title('Magnitude Response of Y2')
plt.xlabel('$\omega$')
plt.ylabel('Real Part')
plt.plot(omega_real,np.absolute(real_fft_filtered))
[<matplotlib.lines.Line2D at 0x7ff165b26518>]

Exercise 4: Yanny or Laurel? Why not both?
For this next exercise, we will visualize the effects of applying LSI systems as filters by looking in the frequency domain. But let's work with something more interesting than toy systems or signals: the infamous "Yanny or Laurel?" audio clip we used in Lab 1. Some people hear Yanny while others hear Laurel. In this activity, we will show using signal processing that both are audible and hopefully shed some light on this auditory illusion.

We have provided two filters in the files filter-one.npy and filter-two.npy in addition to the audio clip. Do not worry about the design of these filters. We will cover filter design in more detail in Lab 6.

a. Use $\textrm{np.load()}$ to load each filter. Note that these coefficients are simply the numerator coefficients of the transfer functions. We will assume a denominator of 1. Plot the magnitude of the DTFT for the audio signal using $\textrm{np.fft.rfft()}$ and the frequency response for each filter using $\textrm{signal.freqz()}$.

b. Apply each filter to the audio clip using $\textrm{signal.lfilter()}$. Listen to the results from each filter. As always, be careful with your volume before listening.

c. What sounds different in each filtered audio clip? Does this explain the auditory illusion? If so, how? If you are having trouble hearing a difference, try changing the playback sampling frequency for the filtered results a little ($\pm$10-20%).

fs,audio = wavfile.read('audiofile.wav') #load the data
print(audio.shape) #one channel
(34752,)
Audio(data = audio,rate = fs) #give it a listen for reference
#Code for 4.a
f1 = np.load('filter-one.npy')
b1 = f1
a1 = 1

real_fft_audio = np.fft.rfft(audio)
omega_real = np.linspace(0,np.pi,len(real_fft_audio))


plt.figure(figsize = (50,50))
plt.subplot(311)
plt.title('Magnitude Response of Audio Signal')
plt.xlabel('$\omega$')
plt.ylabel('Real Part')
plt.plot(omega_real,np.absolute(real_fft_audio))

w,h = signal.freqz(b1,a1) #w = omega/digital frequencies, h = frequency response
plt.subplot(312)
plt.plot(w,20*np.log10(np.absolute(h))) #plot magnitude of frequency response with db-scaling on y-axis
plt.title('Frequency Response of Filter One')
plt.xlabel('$\omega$')
plt.ylabel('Magnitude Response (dB)')


f2 = np.load('filter-two.npy')
b2 = f2
a2 = 1

w,h = signal.freqz(b2,a2) #w = omega/digital frequencies, h = frequency response
plt.subplot(313)
plt.plot(w,20*np.log10(np.absolute(h))) #plot magnitude of frequency response with db-scaling on y-axis
plt.title('Frequency Response of Filter Two')
plt.xlabel('$\omega$')
plt.ylabel('Magnitude Response (dB)')
#Code for 4.b
audio_result1 = signal.lfilter(b1,a1,audio)
audio_result2 = signal.lfilter(b2,a2,audio)

#Make sure to typecast audio before listening to it as follows:
audio_result1 = audio_result1.astype(np.int16)
audio_result2 = audio_result2.astype(np.int16)
#Listen to result 1 here!
Audio(data = audio_result1,rate = fs)
#List to result 2 here!
Audio(data = audio_result2,rate = fs)
Comments for 4.c here:

Question: What sounds different in each filtered audio clip? Does this explain the auditory illusion? If so, how? If you are having trouble hearing a difference, try changing the playback sampling frequency for the filtered results a little ( ± 10-20%).

Solution: The most noticable difference is that the two audio clips have very distinct volumes to them. The next noticable difference is how the audio clips sound at the start and the end of them respectively.

When we use filter one, this filter keeps the lower frequency noises while attenuating the higher frequency sound, and since the magnitude of the input signal has its greatest values when omega is < .35, this filtered signal will of course sound higher in volume. This is also why when I listen to the filtered one signal it sounds a lot like Laurel. Laurel is a sound that contains a lot more lower frequencies (based on just how it sounds when I say it) and therefore would be more easily heard with a low pass filter.

As for the second filtered, when I upped the rate to 1.2*fs, I could hear the sound of Yanni very slightly. The reason why this is true is because we applied a high pass filter to the audio signal as well as increasing the sampling rate will directly increase pitch as well. Yanni is a signal that contains higher frequency components (the Y in Yanni, the "nn" that is lengthed too sounds high pitch) and therefore the high pass filter will make the initial audio signal sound a lot like Yanni instead of Laurel. The volume of this signal also sounds lower than filter one since the respective omega values that this filter keeps are on average a lower magnitude on the audio signal graph.

Sampling and Analog-to-Digital Conversion
For the second half of this lab, we will focus on the process of sampling and storing digital signals. We will begin with some review.

When sampling a continuous time signal, we must be careful to sample at an appropriate frequency. For any bandlimited signal with a maximum frequency of $f_{max}$ or bandwidth $B$, we can guarantee no aliasing if we sample above twice $f_{max}$. This is known as the Nyquist criterion:

$$ f_s > 2B = f_{Nyquist}. $$
How can we relate the analog and digital frequencies before and after sampling? There is a simple equation for that!

$$ \omega_d = \Omega_aT, $$
where $T$ is the sampling period, $\omega_d$ is our digital frequency, and $\Omega_a$ is the analog frequency. Recall that the DTFT of a digital signal is $2\pi$ periodic and our digital frequencies are bounded between $-\pi$ and $\pi$. For example, suppose we have a signal $x(t)$ with $f_{max}$ = 35kHz and we sample at $f_s = 30kHz$. Where will this maximum frequency lie in the digital spectrum?

$$ \omega_d = 2\pi\cdot35000\cdot\frac{1}{30000} $$$$ \omega_d = \frac{7\pi}{3} $$
Clearly, we have sampled below the Nyquist rate, and thus the signal has aliased. The max frequency directly maps to $\frac{7\pi}{3}$; however, by the $2\pi$ periodicity of the DTFT, we will also have a frequency component at $\frac{7\pi}{3}-2\pi = \frac{\pi}{3}$ in the central copy of the DTFT. Excercise 5 will give you some practice with different sampling rates and how to explain aliasing.

In ECE 310, we mainly focus on the sampling part of the analog-to-digital conversion process. Don't forget that in practice we must consider quantization effects. We cannot store every possible analog value of a signal, so we must select a finite number of levels to represent our data. Most simply in ECE 110, we learn about uniform quantizers where the levels are evenly spaced throughout the dynamic range (range of possible values). Each captured sample is "rounded" or quantized to its nearest level. But we should also consider other quantization schemes. For example, what if most of our analog samples densely range between 0V-1V, while a few noisy samples spike up to 5V. A uniform quantizer would lose resolution at lower voltages and accomodate the noisy samples too much. Perhaps there is a better way? Excercise 6 will show you an example of a non-uniform quantizer and let you compare it with a uniform quantizer on a couple test images.

Exercise 5: It's a bird! It's a plane! No, it's just aliasing!
We will now get some hands-on experience with aliasing and sampling effects. Python has a helpful function in the scipy.signal package called $\textrm{signal.chirp()}$ that generates a swept cosine signal. This means we can create a sweeping tone between a start and end frequency. Unfortunately, the documentation for this function is a bit confusing, so let's briefly demonstrate its usage:

Fs = 44100 #sampling rate for audio clip in Hz
t1 = 5 #make clips 5 seconds
t = np.linspace(0,t1,t1*Fs) #make sure to specify the number of points to match desired sampling frequency!!!
f0 = 0 #start frequency (Hz)
f1 = 22050 #end frequency (Hz)

"""
instantanous frequency, f(t) = f0 + (f1-f0)*(t/t1)
"""

chirp_original = signal.chirp(t,f0 = f0, t1 = t1, f1 = f1)
#BE CAREFUL WITH YOUR VOLUME! CHIRP SEQUENCES CAN BE LOUD!

Audio(data = chirp_original,rate = Fs) #give it a listen