Derivation of downsampling in the frequency domain
The expression 

$$s_d[m] = \frac{1}{D}\sum_{k=0}^{D-1}\exp(j2\pi \frac{km}{D})$$

is the Fourier series expansion of the periodic function $s_d[m]$ in the discrete case (i.e. only up to a certain amount of frequencies). 

$s_d[m]$ is a also called comb function, and I have written more about it (including the relation you ask for) in one of my articles.

Furthermore, here's some example Python code that verifies the correctness:

m = np.arange(64)
D = 16
sdm = 1/D*sum(np.exp(2j*np.pi*k*m/D) for k in range(D))
plt.plot(m, sdm.real)




The picture shows the Comb function, where every 16 samples the values becomes one, else zero.