When using Wannier functions to study the electronic structure of multi-parameter Hamiltonians $H^{(\boldsymbol k,\bf λ)}$ carrying a dependence on crystal momentum $\boldsymbol k$ and an additional periodic parameter $\bfλ$, one usually constructs several sets of Wannier functions for a set of values of $\bfλ$. We present the concept of higher dimensional Wannier functions (HDWFs), which provide a minimal and accurate description of the electronic structure of multi-parameter Hamiltonians based on a single set of HDWFs. The obstacle of non-orthogonality of Bloch functions at different $\bfλ$ is overcome by introducing an auxiliary real space, which is reciprocal to the parameter $\bfλ$. We derive a generalized interpolation scheme and emphasize the essential conceptual and computational simplifications in using the formalism, for instance, in the evaluation of linear response coefficients. We further implement the necessary machinery to construct HDWFs from ab initio within the full-potential linearized augmented plane-wave method (FLAPW). We apply our implementation to accurately interpolate the Hamiltonian of a one-dimensional magnetic chain of Mn atoms in two important cases of $\bfλ$: (i) the spin-spiral vector $\boldsymbol q$, and (ii) the direction of the ferromagnetic magnetization $\hat{\boldsymbol m}$. Using the generalized interpolation of the energy, we extract the corresponding values of magneto-crystalline anisotropy energy, Heisenberg exchange constants, and spin stiffness, which compare very well with the values obtained from direct first principles calculations. For toy models we demonstrate that the method of HDWFs can also be used in applications such as the virtual crystal approximation, ferroelectric polarization and spin torques. ; Comment: 23 pages, 11 figures