First of all very similar questions have been asked here and even I have answered some of them, but now I will provide an answer using yet another approach.
Our aim is to show the effect of down-sampling on the DFT.
Let $x[n]$ be a sequence of length $N$. And let's define $y[n]$ by downsampling $x[n]$ by an integer factor $M$ as shown:
$$ x[n] \longrightarrow \boxed{ \downarrow M } \longrightarrow y[n] = x[Mn] \tag{1}$$
$y[n]$ is of length $K = N/M$. At the moment assume $K$ is integer.
I will make use of the relation between DTFT $X(\omega)$ and DFT $X[k]$ for finite length sequence which states:
$$ X[k] = X(\omega)|_{w = \frac{2\pi}{N}k } = X( \frac{2\pi}{N}k ) ~~~,~~~ k = 0,1,2,...,N-1 \tag{2}$$
From the theoretical DTFT definition, it can be shown that the DTFT $Y(\omega)$ of $y[n]$ is given as:
$$ Y(\omega) = \frac{1}{M} \sum_{m=0}^{M-1} X( \frac{ \omega + 2\pi k}{M} ) ~~~,~~~ -\pi \leq \omega < \pi \tag{3}$$
And we define the $K$-point DFT of $Y[k]$ as:
$$ Y[k] = Y(\omega)|_{\omega = \frac{2\pi}{K} k } = Y(\frac{2\pi}{K} k) ~~~,~~~ k = 0,1,...,K-1 \tag{4}$$
Note the range of DFT index $k$ for $Y[k]$. Since $y[n]$ is a $K$-point sequence we have defined a $K$-point DFT of it.
Finally plug Eq(4) into Eq(3)
$$ Y[k] = Y(\frac{2\pi}{K} k) = \frac{1}{M} \sum_{m=0}^{M-1} X( \frac{ \frac{2\pi}{K} k + 2\pi m}{M} ) $$
$$ Y[k] = Y(\frac{2\pi}{K} k) = \frac{1}{M} \sum_{m=0}^{M-1} X( \frac{2\pi}{KM} k + \frac{2\pi}{M} m ) $$
Now $KM = N$ and we replace $\frac{2\pi}{M}$ with $\frac{2\pi}{N}(N/M) = \frac{2\pi}{N}K $ to get
$$ Y[k] = Y(\frac{2\pi}{K} k) = \frac{1}{M} \sum_{m=0}^{M-1} X( \frac{2\pi}{N} k + \frac{2\pi}{N} Km ) $$
$$ Y[k] = Y(\frac{2\pi}{K} k) = \frac{1}{M} \sum_{m=0}^{M-1} X( \frac{2\pi}{N} ( k + Km) ) $$
And finally we recognize that $X( \frac{2\pi}{N} ( k + Km) )$ is actually the $N$-point DFT $X[r]$ of $x[n]$, evaluated at $r = k + Km$ and insert it into the equation:
$$ \boxed{ Y[k] = \frac{1}{M} \sum_{m=0}^{M-1} X[ k + Km] ~~~,~~~k = 0,1,...,K-1 } \tag{5}$$
Note that DFT index k for $Y[k]$ is of modulo-K where as DFT index $r = k + Km$ for $X[r]$ is of modulo-N. Hence Eq(5) can alo be written like:
$$ \boxed{ Y[(k)_K] = \frac{1}{M} \sum_{m=0}^{M-1} X[ (k + Km)_N] } \tag{6}$$
The following is a MATLAB/ OCTAVE code to demonstrate the result:
N = 32; % Length of x[n]
M = 4; % Downsampling factor (integer)
K = N/M; % length of y[n] (assuming N divisble by M)
x = randn(1,N); % x[n]
y = x(1:M:end); % y[n] = x[Mn]
X = fft(x,N); % N-point DFT X[r] of x[n]
Y = zeros(1,length(y)); % K-point DFT Y[k] of y[n]
k = 0:length(Y)-1; % DFT index range k
for m = 0:M-1 % implements Eq.6 above
Y = Y + X( mod(k+K*m,N)+1 );
end
Y = Y/M;
Y2 = fft(y,K); % Comput Y[k] directly from y[n] for checking.
figure,stem(abs(Y)) % display that they are the same..
hold on
stem(abs(Y2),'r+');
