Projection into Spherical Harmonic Basis

My understanding is that I should be able to just multiply each of my function's samples with the basis functions evaluated at the sample's XYZ direction, then accumulate the results to get a SH expansion for the function samples.

Practically, yes. In the continuous formalism, the expansion factors a for function f are obtained by

<img src="http://www.codecogs.com/eq.latex? a_{lm} = \int_0^{2\pi}\int_0^\pi Y^*_{lm}(\theta,\phi)f(\theta,\phi)sin(\theta)d\theta d\phi" />

Where the sin factor is the Jacobian. You can use which ever numerical integration you find most fit. For smooth functions, you could use the Simpson rule for both angles in a nested way. This is enabled by the Fubini theorem for smooth integrands. You might also use Gauss family of quadrature (e.g. Gauss-Legendre) to reduce the number of evaluation points.

My question is, how well will this work if I only have a few samples of my function, possibly with irregularly spaced samples? Will the SH basis smoothly interpolate just a few (e.g. 3 or 4 values) samples, or do I need to do a monte carlo integration scheme where I sample a separate function that uses the nearest sampled neighbor for the query direction.

I don't know the precise answer to this, but what you are most likely doing with SH basis is not interpolation, but approximation. When you evaluate your function in the truncated SH basis, there is no guarantee that it will interpolate your sample points. If you have only a few samples, you probably get aliasing issues as with the more traditional Fourier series.

I don't know the precise answer to this, but what you are most likely doing with SH basis is not interpolation, but approximation. When you evaluate your function in the truncated SH basis, there is no guarantee that it will interpolate your sample points. If you have only a few samples, you probably get aliasing issues as with the more traditional Fourier series.

OK, I see the distinction. In theory you could just increase the order of the expansion to better approximate the samples, right?

OK, I see the distinction. In theory you could just increase the order of the expansion to better approximate the samples, right?

Do you specifically need interpolation property, i.e., that the approximate representation will take on the sample values at their original locations? Or do you want to get a good overall representation of the function of which you have some limited set of samples?

Technically speaking, to get better approximation when increasing the order of the expansion, the original function needs to be in the L² function space over a sphere.

I'm using the spherical harmonics to represent the radiation pattern of sound sources. For these I may have measurements of frequency response at different directions relative to the source, such as for a loudspeaker. From these few measurements I need to be able to determine the radiation for an arbitrary direction that should roughly correspond to the measurements. The radiation is generally smoothly varying and should be easily represented by a low-order expansion. Interpolation is preferred, but approximation is ok. The most important part is fast evaluation of the expansion.

I think I will compute the expansion for increasing orders, until the error relative to the samples drops below some threshold, e.g. 1dB, or a max allowed order is reached.

It turns out that just using a few sample points to do the projection produces poor results that don't approximate or interpolate well. I think I will have to do some sort of numerical integration to do the projection so that it is more stable.

Edit: I am picking N random direction samples, interpolating between the closest 3 of the original measured samples, and then integrating over those instead. It seems to produce much better results and convergence.

I am picking N random direction samples, interpolating between the closest 3 of the original measured samples, and then integrating over those instead. It seems to produce much better results and convergence.

This is strange. If you just interpolate the original measurement points to create more data, you are creating data that is linearly dependent on the original data. No new information should emerge.

Did you take the original measurements at the correct quadrature points, i.e., the points where the integration weights are defined? Keep in mind that you don't "integrate over sample points". Surface integral over a set of isolated points vanishes. What you do is make and approximation to a surface integral in the form

int_S f( r ) dS = sum_n w_n*f(r_n)

where w_n are the weights and r_n the nodes of the quadrature rule.

For 2D surface integrals with smooth surface S, you can use the Fubini's theorem to write your spherical integration as

int_0^pi int_0^2pi f(theta, phi)*sin(theta) dphi dtheta = sum_n sum_m w_n*w_m*f(theta_n, phi_m)*sin(theta_n) + E

where E is some error that we are willing to neglect and the weights and nodes are now for some integration over finite 1D interval (trapezoidal, Simpson's, Gaussians, etc.)

The nodes theta_n and phi_m should coincide with your measurement points. If you can choose your measurement points, you should choose

them to match the nodes of the quadrature rule.