In this lab, you'll create sinusoids at a variety of sampling rates, both above and below Nyquist, and compare the resulting stem plots. Then you'll downsample a sinusoid, and interpolate it back to the original sampling rate, using a variety of different interpolation kernels.

In order to make sure everything works, you might want to go to the command line, and run

```
pip install -r requirements.txt
```

This will install the modules that are used on the autograder, including numpy, h5py, and the gradescope utilities.

In [48]:

```
import submitted
import importlib
import numpy as np
import matplotlib.pyplot as plt
import h5py
```

In this first part of the MP, you will create a sinusoid, and sample it at different sampling rates. The four parameters will be

`frequency`

$f=$ frequency of the sinusoid (Hertz)`phasor`

$z=re^j\phi=$ phasor encoding its amplitude and phase`duration`

$d=$ duration (seconds)`samplerate`

$F_s=$ sampling rate (samples/second)

The returned signal should be $$x[n] = \Re\{ z e^{j2\pi nf/F_s}\}$$

In [49]:

```
import importlib
importlib.reload(submitted)
help(submitted.sinusoid)
```

First, let's generate a sinusoid with plenty of samples, so we make sure it plots correctly. We'll use `stem`

to plot it, so we can see each sample.

In [50]:

```
importlib.reload(submitted)
timeaxis, signal = submitted.sinusoid(4, 1, 1, 32)
plt.figure(figsize=(14, 5))
plt.stem(timeaxis, signal)
plt.xlabel('Time (seconds)')
plt.title('4Hz cosine, sampled at 32 samples/second')
```

Out[50]:

Text(0.5, 1.0, '4Hz cosine, sampled at 32 samples/second')

Now let's see what happens when we gradually lower the sampling rate.

In [51]:

```
fig, ax = plt.subplots(4,1, figsize=(14,10))
samplerates = [16, 10, 8, 6]
for k in range(4):
timeaxis, signal = submitted.sinusoid(4, 1, 1, samplerates[k])
ax[k].stem(timeaxis, signal)
ax[k].set_title('4Hz sinusoid with a sampling rate of %d Hz'%(samplerates[k]))
```

Notice that, as long as you have at least two samples per period ($F_s\ge 2f$), you can reconstruct the original cosine (e.g., there are four sign-reversals in one second, so you can see that it's a 4Hz sinusoid).

When $F_s < 2f$, you can no longer see the real underlying cosine. It has been "aliased" to a new "alias frequency," $f_a=F_s-f$. That's because, for any integer $n$,

$$\cos\left(\frac{2\pi f n}{F_s}\right) = \cos\left(\frac{2\pi (F_s-f) n}{F_s}\right)$$The highest frequency at which the cosine can be is called the "Nyquist rate", and it is $$F_N = \frac{1}{2}F_s$$

Let's adjust the sampling frequency so that 4Hz is just barely above or just barely below Nyquist.

In [52]:

```
fig, ax = plt.subplots(5,1, figsize=(14,10))
samplerates = [9, 8, 7, 6, 5]
for k in range(5):
timeaxis, signal = submitted.sinusoid(4, 1, 1, samplerates[k])
ax[k].stem(timeaxis, signal)
alias = min(4, samplerates[k]-4)
ax[k].set_title('4Hz, at $F_s$=%d Hz, looks like %dHz'%(samplerates[k], alias))
ax[k].set_ylim(-1,1)
fig.tight_layout()
```

Notice that, when $f > F_s/2$, it looks as though the sinusoid has a frequency of $F_s-f$. This is called its "aliased frequency." When the sampling rate gets even lower, aliasing can be even worse. The aliased frequency in general is

$$f_a = \min\left(f \mbox{mod} F_s, (F_s-f)\mbox{mod} F_s\right)$$where $\mbox{mod}$ means modulo (written in python as `f % Fs`

). Let's see some examples.

In [53]:

```
fig, ax = plt.subplots(5,1, figsize=(14,10))
samplerates = [4.5, 4, 3.5, 3, 2.5]
for k in range(5):
timeaxis, signal = submitted.sinusoid(4, 1, 2, samplerates[k])
ax[k].stem(timeaxis, signal)
alias = min(4 % samplerates[k], (samplerates[k]-4) % samplerates[k])
ax[k].set_title('4Hz, at $F_s$=%g Hz, looks like %gHz'%(samplerates[k], alias))
ax[k].set_ylim(-1,1)
fig.tight_layout()
```

Exactly the same thing happens to sine waves as to cosines, except that:

- When $f=F_s/2$ exactly, the sine wave disappears.
- When $F_s/2 < f < F_s$, the sign of the sine reverses. That's because $$\sin\left(\frac{2\pi f n}{F_s}\right) = -\sin\left(\frac{2\pi (F_s-f) n}{F_s}\right)$$

In order to plot sine waves, we'll use a phasor of $-j=e^{-j\pi/2}$, to give it a phase of $-\pi/2$.

In [54]:

```
fig, ax = plt.subplots(5,1, figsize=(14,10))
samplerates = [10,9,8,7,6]
for k in range(5):
timeaxis, signal = submitted.sinusoid(4, -1j, 1, samplerates[k])
ax[k].stem(timeaxis, signal)
alias = min(4 % samplerates[k], (samplerates[k]-4) % samplerates[k])
ax[k].set_title('4Hz, at $F_s$=%g Hz, looks like %gHz'%(samplerates[k], alias))
ax[k].set_ylim(-1,1)
fig.tight_layout()
```

These two equations: $$\cos\left(\frac{2\pi f n}{F_s}\right) = \cos\left(\frac{2\pi (F_s-f) n}{F_s}\right)$$ $$\sin\left(\frac{2\pi f n}{F_s}\right) = -\sin\left(\frac{2\pi (F_s-f) n}{F_s}\right)$$

...can be combined to give this equation: $$\Re\left\{z \exp\left(j\frac{2\pi f n}{F_s}\right)\right\}=\Re\left\{z^* \exp\left(j\frac{2\pi (F_s-f) n}{F_s}\right)\right\}$$

In other words,

- If $f\mbox{mod}F_s < (F_s-f)\mbox{mod}F_s$, then the aliased frequency is $$f_a=f\mbox{mod}F_s$$ and the phasor of the aliased sinusoid is the same as the phasor of the true sinusoid: $$z_a = z$$
- If $f\mbox{mod}F_s > (F_s-f)\mbox{mod}F_s$, then the aliased frequency is $$f_a=(F_s-f)\mbox{mod}F_s$$ and the phasor of the aliased sinusoid is the complex conjugate the phasor of the true sinusoid: $$z_a = z^*$$

For example, let's look at the signal $$x(t) = \cos\left(2\pi ft-\frac{\pi}{4}\right)$$

This is a cosine delayed by $\pi/4$ radians. It's positive for about $3/8$ of a period, then negative for $1/2$ period, and so on.

In [55]:

```
fig, ax = plt.subplots(5,1, figsize=(14,10))
samplerates = [30,11,10,6,5]
for k in range(5):
phasor = np.exp(-1j*np.pi/4)
timeaxis, signal = submitted.sinusoid(4, phasor, 1, samplerates[k])
ax[k].stem(timeaxis, signal)
alias = min(4 % samplerates[k], (samplerates[k]-4) % samplerates[k])
ax[k].set_title('4Hz, at $F_s$=%g Hz, looks like %gHz'%(samplerates[k], alias))
ax[k].set_ylim(-1,1)
fig.tight_layout()
```

Sometimes it's useful to explicitly calculate the aliased frequency and aliased phasor (HINT: you probably want to use the `np.mod`

method):

In [56]:

```
importlib.reload(submitted)
help(submitted.compute_aliasing)
```

In [57]:

```
importlib.reload(submitted)
samplerates = np.array([30,11,10,6,5])
freqs = np.repeat(4, 5)
phasors = np.repeat(np.exp(-1j*np.pi/4), 5)
aliased_freqs, aliased_phasors = submitted.compute_aliasing(freqs, phasors, samplerates)
print(aliased_freqs)
print(aliased_phasors)
```

If we plot a sinusoid at the aliased frequency and aliased phasor, it should look exactly the same as the original sinusoid.

In [58]:

```
importlib.reload(submitted)
time1, original = submitted.sinusoid(freqs[-1],phasors[-1],2,samplerates[-1])
time2, alias = submitted.sinusoid(aliased_freqs[-1],aliased_phasors[-1],2,samplerates[-1])
time3, oops = submitted.sinusoid(aliased_freqs[-1],aliased_phasors[0],2,samplerates[-1])
```

In [59]:

```
fig, ax = plt.subplots(3,1, figsize=(14,10))
ax[0].stem(time1, original)
ax[0].set_title('Sinusoid at %dHz, %dHz samplerate'%(freqs[-1],samplerates[-1]))
ax[1].stem(time2, alias)
ax[1].set_title('Sinusoid at %dHz, %dHz samplerate'%(aliased_freqs[-1],samplerates[-1]))
ax[2].stem(time3, oops)
ax[2].set_title('Aliased freq, but original phasor')
fig.tight_layout()
```

Remember that, if you have a signal that's periodic with a period $T_0$, then you can find its Fourier series coefficients using the Fourier analysis formula: $$X_k = \frac{1}{T_0}\int_{0}^{T_0} x(t) e^{-j\frac{2\pi kt}{T_0}}dt$$

We can get a pretty good measurement of $X_k$ on a computer by setting the sampling rate high enough, replacing the integral by a sum, and multiplying the integrand by $F_s$ (to represent $dt=1/F_s$). After simplifying a little, we get: $$X_k = \frac{1}{N_0} \sum_{n=0}^{N_0-1} x[n]e^{-j\frac{2\pi kn}{N_0}}$$ where $$N_0 = T_0F_s$$ During the rest of this MP, it will be useful to have a good Fourier analysis function, so let's write it.

In [60]:

```
importlib.reload(submitted)
help(submitted.fourier_analysis)
```

Remember that, if a signal is periodic with period $N_0$, then it's also periodic with any period that's a multiple of $N_0$. We can use that fact to make the plot easier to read, e.g., by using four periods of a cosine as if they were just one period; that way the coefficient $X_1$ will be plotted as if it were $X_4$.

Let's choose a cosine with a really long period, say, 50 samples, and then compute the Fourier series coefficients over a space of 200 samples.

In [61]:

```
n = np.arange(200)
signal = np.cos(2*np.pi*n/50)
k = np.arange(41)/4
coefficients = submitted.fourier_analysis(signal, 41)
```

In [62]:

```
fig, ax = plt.subplots(2,1,figsize=(14,6))
ax[0].stem(n, signal)
ax[0].set_ylabel('$x[n]$')
ax[0].set_xlabel('$n$')
ax[1].stem(k, np.abs(coefficients))
ax[1].set_ylabel('$|X_k|$')
ax[1].set_xlabel('$k$')
```

Out[62]:

Text(0.5, 0, '$k$')

Let's define the triangle function as follows: $$h[n]=\begin{cases}1-\frac{|t|}{T}&-T+1\le n\le T-1\\0&\mbox{otherwise}\end{cases}$$

This function will be useful when we try to interpolate signals, in part 5. Go ahead and write a function to compute it. Notice that you only need to include the nonzero samples, not the zero-valued samples.

In [63]:

```
importlib.reload(submitted)
help(submitted.triangle)
```

In [64]:

```
importlib.reload(submitted)
t_timeaxis, triangle = submitted.triangle(25)
```

In [65]:

```
plt.figure(figsize=(14,4))
plt.plot(t_timeaxis, triangle)
```

Out[65]:

[<matplotlib.lines.Line2D at 0x7f852e7f9a60>]

We are now going to simulate a discrete-to-continuous time converter by downsampling our 50-sample cosine by a factor of 25, and then upsampling it by the same factor of 25. The upsampling operation will use the same four types of interpolation that are used in typical D/C converters: rectangle, triangle, spline, and sinc.

First, let's downsample:

In [66]:

```
lowrate_signal = signal[::25]
lowrate_n = n[::25]
plt.figure(figsize=(14,4))
plt.stem(lowrate_n, lowrate_signal)
plt.title('Cosine, sampled with only two samples per period')
```

Out[66]:

Text(0.5, 1.0, 'Cosine, sampled with only two samples per period')

Now let's pretend we're interpolating it into continuous time. The general formula for interpolation is

$$x(t) = \sum_{n=-\infty}^\infty x[n] h\left(t-nT\right)$$where $T$ is the sampling period.

In order to simulate D-to-C conversion, we will use $t$ to represent time in the highrate signal (e.g., $0\le t<200$) instead of time in the continuous-time signal. We'll use $n$ to represent time in the lowrate signal (e.g., in the example, $0\le n< 8$), while $T$ is the ratio between the two sampling rates (e.g., $T=25$).

In [67]:

```
importlib.reload(submitted)
help(submitted.interpolate)
```

Notice those last few lines from the docstring. The triangle function has a length of $2T-1$. The lowrate signal has a length of $N$. If you put in a triangle function in place of every lowrate sample, you'd get a total length of $NT+(T-1)$, but that's not what we want -- what we want is a signal of length $NT$. In order to make that happen, we'll use modulo indexing, i.e., the interpolation formula should be implemented as

`for n in range(N):`

You might find the function `np.mod`

to be useful.

In [68]:

```
importlib.reload(submitted)
highrate_signal = submitted.interpolate(lowrate_signal, 25, t_timeaxis, triangle)
```

In [69]:

```
fig,ax = plt.subplots(3,1,figsize=(14,10))
ax[0].stem(n,signal)
ax[0].set_title('Original Signal')
ax[1].stem(lowrate_n, lowrate_signal)
ax[1].set_title('Lowrate "Discrete Time" Signal')
ax[2].stem(n, highrate_signal)
ax[2].set_title('Highrate "Continuous Time" Signal using Triangle Interpolator')
```

Out[69]:

Text(0.5, 1.0, 'Highrate "Continuous Time" Signal using Triangle Interpolator')

Well, it doesn't look very much like a cosine! In fact, it's just a linear interpolation between the lowrate samples.

In order to better understand what's going on, let's look at its spectrum. In the following plot you should see that most of the spectral energy is at $X_1$, but there is significant nonzero energy also at $X_3$, $X_5$, and so on:

In [70]:

```
highrate_coefficients = submitted.fourier_analysis(highrate_signal, 41)
```

In [71]:

```
fig, ax = plt.subplots(2,1,figsize=(14,6))
ax[0].stem(k, np.abs(coefficients))
ax[0].set_title('Fourier Series Coefficients of the Original Sinusoid')
ax[1].stem(k, np.abs(highrate_coefficients))
ax[1].set_title('Fourier Series Coefficients of the Triangle-Interpolated Sinusoid')
```

Out[71]:

Text(0.5, 1.0, 'Fourier Series Coefficients of the Triangle-Interpolated Sinusoid')

Rectangle interpolation is much, much worse than triangle interpolation, but it's much easier to implement in hardware, so it's actually much more common in the real world. Let's try it, and see how it turns out:

In [72]:

```
importlib.reload(submitted)
help(submitted.rectangle)
```

In [73]:

```
importlib.reload(submitted)
r_timeaxis, rectangle = submitted.rectangle(25)
```

In [74]:

```
plt.figure(figsize=(14,4))
plt.plot(r_timeaxis, rectangle)
```

Out[74]:

[<matplotlib.lines.Line2D at 0x7f8531937640>]

In [75]:

```
importlib.reload(submitted)
highrate_signal = submitted.interpolate(lowrate_signal, 25, r_timeaxis, rectangle)
```

In [76]:

```
fig,ax = plt.subplots(3,1,figsize=(14,10))
ax[0].stem(n,signal)
ax[0].set_title('Original Signal')
ax[1].stem(lowrate_n, lowrate_signal)
ax[1].set_title('Lowrate "Discrete Time" Signal')
ax[2].stem(n, highrate_signal)
ax[2].set_title('Highrate "Continuous Time" Signal using Rectangle Interpolator')
```

Out[76]:

Text(0.5, 1.0, 'Highrate "Continuous Time" Signal using Rectangle Interpolator')

In [77]:

```
highrate_coefficients = submitted.fourier_analysis(highrate_signal, 41)
```

In [78]:

```
fig, ax = plt.subplots(2,1,figsize=(14,6))
ax[0].stem(k, np.abs(coefficients))
ax[0].set_title('Fourier Series Coefficients of the Original Sinusoid')
ax[1].stem(k, np.abs(highrate_coefficients))
ax[1].set_title('Fourier Series Coefficients of the Rectangle-Interpolated Sinusoid')
```

Out[78]:

Text(0.5, 1.0, 'Fourier Series Coefficients of the Rectangle-Interpolated Sinusoid')