Spherical equivalent of Nyquist frequency

Question

Let $\phi$ be a scalar function defined on the surface of a sphere. I have samples of $\phi$ at various locations on the sphere. I want to apply a spherical harmonic transform. I know that $\phi$ is 'band-limited' in the sense that it can be accurately represented by a truncated spherical harmonic expansion, which only includes terms up to a certain angular order, $l_{max}$. (Using this notation: http://mathworld.wolfram.com/SphericalHarmonic.html.)

The sampling is much denser than I require to accurately perform the spherical harmonic transform. I would like to discard some of the samples before doing the transform, so that the computation will be faster. How many samples do I need (in terms of $l_{max}$) to ensure 'good' results?

I realise that the answer will depend on the definition of 'good' and the distribution of the sample points (in my case, the samples are randomly distributed). However I am looking for a result like Nyquist's theorem from which I can make a safe estimate of how many samples to keep.

Answer 1

Looks like a potential application for blue noise also known as Poisson disk sampling, which is random placement of samples but with a guaranteed minimum distance between sample locations. I think that would give more accurate transform results than independently located random samples.

Various spherical harmonic transform algorithms with specific sampling schemes take as input $4l_\text{max}^2,$ $2l_\text{max}^2$ or $l_\text{max}^2$ samples, see Khalid et al., 2014, "An Optimal-Dimensionality Sampling Scheme on the Sphere with Fast Spherical Harmonic Transforms", arXiv:1403.4661. There's no going lower than $l_\text{max}^2$ because it exactly equals the degrees of freedom in the spherical harmonic domain.

Perhaps you can determine which factor you need to use in front of $l_\text{max}^2$ for your distribution of sample locations, to get reliable results. This can be done based on transform–inverse transform numerical accuracy testing

Answer 2

I'm a little confused about what you are asking and how it relates to the set of spherical basis functions. Those are the same one that electron orbitals are based on, right? I have never understood them to be constrained to the surface of a sphere. Then again, that is not my usual stomping grounds.

Nevertheless, I think I can generalize your question and provide an answer.

$\phi$ is a poor choice for function name in this context. Yes, I know it is the standard for "potential function" in many contexts, but within spherical coordinates it has a different meaning. I will use $lon$ and $lat$ for your point location, and retain $\phi$ for the value so there is no confusion.

You have a large collection of points that fit the model:

$$ \phi_n = \phi(lat_n,lon_n) $$

You would like to find $\phi$ as a linear combination of pre-defined basis functions.

$$ \phi(lat,lon) = \sum_{k=0}^{K-1} c_k B_k(lat,lon) $$

Presumably, $B_0(lat,lon) = 1$, but that is not necessary, but I wanted to throw in the "DC" term to help orient folks in relation to a DFT.

Anyway, for each $k$, you have a computable function definition of $B_k(lon,lat)$.

This is now a standard best fit problem. You want to find the values of the $c_k$s for the best fit of $\phi$ to your data.

Set up the following arrays:

$$ B = \begin{bmatrix} B_0(lat_0,lon_0) & B_1(lat_0,lon_0) & B_2(lat_0,lon_0) & ... & B_{K-1}(lat_0,lon_0) \\ B_0(lat_1,lon_1) & B_1(lat_1,lon_1) & B_2(lat_1,lon_1) & ... & B_{K-1}(lat_1,lon_1) \\ : & : & : & ... & : \\ B_0(lat_{N-1},lon_{N-1}) & B_1(lat_{N-1},lon_{N-1}) & B_2(lat_{N-1},lon_{N-1}) & ... & B_{K-1}(lat_{N-1},lon_{N-1}) \\ \end{bmatrix} $$

$$ C = \begin{bmatrix} c_0 \\ c_1 \\ c_2 \\ : \\ c_{K-1} \\ \end{bmatrix} $$

$$ P = \begin{bmatrix} \phi_0 \\ \phi_1 \\ : \\ \phi_{N-1} \\ \end{bmatrix} $$

$B$ (NxK) and $P$ (Nx1) are known from your data points and basis function evaluations. $C$ (Kx1) contains the unknowns you are solving for.

Ideally,

$$ BC = P $$

If $N=K$, simply multiply both sides by $B^{-1}$ and you have a solution.

$$ C = B^{-1}P $$

Generally, you want more points than that, but that is the minimum requirement unless you want to be underfitted.

If $N>K$, there is a standard approach to find the best fit $C$:

$$ B^* B C = B^* P $$

$$ C = (B^* B)^{-1} B^* P $$

This is basically an interpolation technique. Interestingly, on the domain of a sphere surface, interpolation vs extrapolation seems to be meaningless. You will want your points to be as spread out as possible as to not to leave any large areas of your surface unsampled. The size of $B^* B$ is KxK, and its inversion is probably your biggest computational load depending on how intensive your basis functions are to evaluate at the surface points.

In a conventional DFT, $N=K$ and $B^{-1}$ is $\frac{1}{N}B^*$ so the inverse is implicitly solved.

The concept of a "Nyquist frequency" doesn't have a direct correlation in this scenario in the general case. It may in terms of your choice of $B_k$ functions, but it would still be irrelevant because in your situation you have randomly distributed points instead of regularly spaced ones.

Your ultimate answer is then you are going to want to have at least as many sample points as you have basis functions, but are likely going to want significantly more. The number of points you select will not impact the size of the inverse that needs to be taken.

Hope this helps.

---

This is in response to Fat32's comment.

Looking at the reference given by the OP, it seems the $B_k$s are defined by equations (18) through (33), though it looks like some others are listed as well. Of course, the list can be extended. Since the OP's points are randomly distributed, I don't see any other feasible method than the one I provided. Note that the first one is indeed a "DC" function.

It does raise the question though if there is a set of definable and spread out locations for which matrix $B^{-1}$ is a multiple of $B^*$ when $N=K$, or $ B^* B $ is a multiple of the identity matrix when $N>K$ . (Calling Olli) Just because continuous functions are orthogonal does not necessarily mean that discrete sampled versions of them are. With the trig functions (the DFT) they are, with Legendre polynomials, for instance, they aren't. [Correction: Being diagonal is sufficient for a trivial inverse.]

Regardless, for a fixed set of sample locations $(B^* B)^{-1} B^*$ can always be pre-calculated. Then finding the coefficients is merely a simple matrix multiplication, just like with the DFT. This is not what the OP has. [Correction: It is possible to select a set of points with $N \ge K$ so the matrix is not invertible. With a random selection this is not likely, but the columns of $B$ may end up not being very orhogonal which would lead to poor numerical behavior. Thanks for the paper reference, Olli. That looks spot on and the OP should definitely read it.]

The concept of a "Nyquist frequency" is strictly dependent on the spacing of evenly spaced sampling locations on a 1D scale, and completely independent of the number of points in the frame. We tend to think of it as "the last valid basis function before we get into aliases" (assuming the upper half is considered negative frequencies and are used first). The OP's situation is similar to unevenly sampled points of a signal, so I really don't see the Nyquist concept as being meaningfully translatable.

I interpreted his question as being more of a the necessary number of points for a good answer and just speculation that it was analogous to the Nyquist.

I do think that it is instructive to understand that the DFT is indeed an example of a best fit solution just like this.

Answer 3

The paper Computing Fourier Transforms and Convolutions on the 2-Sphere (Driscoll & Healy, 1994) develops, in section 4, a sampling theorem analogous to Nyquist's theorem, showing that for a band-limited function on the sphere (with non-zero coefficients up to $\ell_{max}$), the function can be exactly transformed into spherical harmonics (and then recovered at any point) if samples are provided on a grid of points whose longitude and latitude are equally spaced ($\theta_{j} = \pi j/2\ell_{max}$, $\phi_{k} = \pi k/2\ell_{max}$).