Upsampling implementation
After upsampling with zero-insert, the higher frequency images need to be filtered out, which will serve the purpose of growing the zero insert values to the correct interpolated value in between samples. Once this is done, subsequent decimation operations can be properly done since every sample will then properly represent the signal.

For a signal that occupies BW $\in [0, \omega_1]$ with original sampling frequency at $\omega = 2\pi$ and is interpolated by $I$, after zero insert, the new sampling frequency will be at $I2\pi$ and $I-1$ images $2\omega_1$ wide will be centered at every multiple of the original sampling rate ($N2\pi$) for $N \in [1,2,\ldots I-1]$ within the new interpolated sampling rate $\omega \in [0, I2\pi]$.

In terms of the normalized angular output frequency with sampling rate $\omega = 2\pi$, the bandwidth of the original signal will be BW $\in [0, \omega_1/I]$ and the $I-1$ images will each be $2\omega_1/I$ wide and centered at multiples of $N2\pi/I$ for $N \in [1,2,\ldots I-1/I]$.

scipy.signal.firls() is a good choice for creating an interpolation filter and supports multiband filters to maximize rejection at the image frequencies. Decimation will follow a similar image folding where the $D-1$ locations that would be folded in should also be considered for filtering prior to decimation. 

For further details see related post:
Choosing right cut-off frequency for a LP filter in upsampler