We know $\left\{\left(\exp\left(-\frac{2\pi imk}{M}\right)\right)_{k=0}^{M-1}\right\}_{m=0}^{M-1}$ is an orthogonal basis of the space of complex sequences of length $M$. The orthogonality can be proved by the well-known equation
$$ \frac{1}{M}\sum_{k=0}^{\textrm{M}-1}\exp\left(-2i\pi \frac{mk}{M}\right)=\begin{cases}1\quad & m|M=0 \\ 0\quad &m|M\neq0 \end{cases} $$ This equation is simple to prove because $m$ and $k$ are all integers and the sum can be calculated using geometric series sum formula. Now I want to prove a more general result where $k$ and $m$ are vectors of length $d$ and $M$ is $d\times d$ integer nonsingular matrix.
$$ \frac{1}{\textrm{det}(M)}\sum_{k\in \mathcal{N}(M^T)}\exp\left(-2i\pi m^TM^{-1}k\right)=\begin{cases}1\quad & m\in\mathcal{L}(M^T) \\ 0\quad &m\notin\mathcal{L}(M^T),m\in\mathbb{Z}^d \end{cases} $$ where $\mathcal{L}(M^T)=\{M^Tn;n\in\mathbb{Z}^d\}$ and $\mathcal{N}(M^T)=\{ M^Tt; t\in [0,1)^d\}\cap\mathbb{Z}^d$. $\mathcal{L}(M^T)$ is said to be a lattice generated by $M^T$.
We can see that this equation is equivalent to the first one if $d=1$. Take $M=4$ for example. Then $\mathcal{L}(M^T)=\{0,1,2,3\}$ and $\mathcal{N}(M^T)=\{0,\frac{1}{4},\frac{2}{4},\frac{3}{4}\}$.
I find a proof that uses fourier transform to convert the problem into one in frequency domain, sketched in the problem $12.29$, page 657, multirate systems and filter banks by Vaidyanathan. I am wondering if there exists any elementary proof similar to the $1$-dimensional case where a geometric series sum formula is used.