Trying to understand downsampling and then upsampling
Downsampling can be defined as follows:

$s_{out}[n]= s_{in}[n\cdot M]$ (taking every M-th signal sample)

After applying Z-transform we obtain:

$\bar{s}_{out}[z]=\dfrac{1}{M}\sum_{m=0}^{M-1}\bar{s}_{in}\left[ z^{1/M} \cdot w_{M}^{m} \right]$, where $w_{M}=e^{-2\pi i /M}$

By substituting $z=e^{2\pi i f_n}$ ($f_n $ is the normalised frequency: $f_n = \dfrac{f}{f_s} \in [0, 0.5]$

This yields following relationship between input and output spectrum:

$\hat{s}_{out}[f_n]= \dfrac{1}{M}\sum_{m=0}^{M-1}\hat{s}_{in}\left[ (f_n-m)/M \right]$

You can understand that as a fact that your spectrum is being copied M times (even outside of your original sampling frequency) and stretched by factor of M.

Upsampling on the other hand is governed by following equations:

$s_{out}[n]=\begin{cases} s_{in}(n/M) &amp; \text{n} \; mod \;\text{M} =0 \\ 0 &amp; \text{otherwise}\end{cases}$

After applying Z-transform:

$ \sum_{n} s_{out}[n]\cdot z^{-n} = \sum_{n} s_{in}[n]\cdot z^{-nM}$, which gives: $\bar{s}_{out}[z]=\bar{s}_{in}[z^M] $

By doing same substitution of a unit circle instead of z variable:

$\hat{s}_{out}[f_n] = \hat{s}_{in}[M\cdot f_n] $

You can imagine that this equation is squeezing your spectrum by factor of M, bringing at the same time whatever is outside of your sampling frequency range.

Plots you've obtained are great example how it works! 

P.S. There is a dft tag.