I have been trying to find an implementation of Spherical harmonics for higher dimensional data but I couldnt find anything in Sage, Mathematica, Matlab. Does anyone have any idea of a standard/fast implementation ?
MORE DETAILS: I am adding more information about the expression that I am interested in. This comes from equation 13 from this paper on the eigenfunctions of dot product kernels.
Moreover, the Legendre Polynomials may be expanded into an orthonormal basis of spherical harmonics $Y^d_{n,j}$ by the Funk-Hecke equation to obtain: $$P^n_d(x.y) = \frac{|S_{d-1}|}{N(d,n)} \sum_{j=1}^{N(d,n)} Y_{n,j}^d(x) Y_{n,j}^d(y)$$
Here $x,y$ are vectors in $R^d$ with unit norm and I am interested in computing the spherical harmonic $Y_{n,j}^d$ at the points $x$ and $y$.