We study unitary orthonormal bases in the sense of Pimsner and Popa for inclusions $(\mathcal{B}\subseteq \mathcal{A}, E),$ where $\mathcal{A}, \mathcal{B}$ are finite dimensional von Neumann algebras and $E$ is a conditional expectation map from $\mathcal{A}$ onto $\mathcal{B}$. It is shown that existence of such bases requires that the associated inclusion matrix satisfies a spectral condition forcing dimension vectors to be Perron-Frobenius eigenvectors and the conditional expectation map preserves the Markov trace. Subject to these conditions, explicit unitary orthonormal bases are constructed if either one of the algebras is abelian or simple. They generalize complex Hadamard matrices, Weyl unitary bases, and a recent work of Crann et al which correspond to the special cases of $\mathcal{A}$ being abelian, simple, and general multi-matrix algebras respectively with $\mathcal{B}$ being the algebra of complex numbers. For the first time $\mathcal{B}$ is more general. As an application of these results it is shown that if $(\mathcal{B}\subseteq \mathcal{A}, E),$ admits a unitary orthonormal basis then the Connes-Størmer relative entropy $H(\mathcal{A}_1|\mathcal{A})$ equals the logarithm of the square of the norm of the inclusion matrix, where $\mathcal{A}_1$ denotes the Jones basic construction of the inclusion. As a further application, we prove the existence of unitary orthonormal bases for a large class of depth 2 subfactors with abelian relative commutant.