SHT and CHT

[narahahn]

According to Eq. (3.34) in Rafaely, 2015, the forward and inverse discrete spherical harmonic transform are defined as

fnm = Y^-1 * f
f = Y * fnm

where fnm denotes the SHT coefficient vector, f the sound field vector, Y the spherical harmonics matrix, and Y^-1 its pseudo-inverse. In our toolbox, we are using the same matrix Y (see sht_matrix), but it is not used in a consistent way.

sfa-numpy/examples/modal_beamforming_rigid_array.py

Lines 20 to 26 in 43329a7

	# pressure on the surface of a rigid sphere for an incident plane wave
	bn = micarray.modal.radial.spherical_pw(N, k, r, setup='rigid')
	D = micarray.modal.radial.diagonal_mode_mat(bn)
	Y_p = micarray.modal.angular.sht_matrix(N, azi, elev)
	Y_pw = micarray.modal.angular.sht_matrix(N, azi_pw, np.pi/2)
	p = np.matmul(np.matmul(np.conj(Y_pw.T), D), Y_p)
	p = np.squeeze(p)

Here, Y_p is used for the synthesis (inverse transform) which follows the convention introduced in [Rafaely, 2015].

sfa-numpy/examples/modal_beamforming_rigid_array.py

Lines 28 to 40 in 43329a7

	# plane wave decomposition using modal beamforming
	Y_p = micarray.modal.angular.sht_matrix(N, azi, elev, weights)
	# get SHT matrix for a source ensemble of azimuthal plane waves
	azi_pwd = np.linspace(0, 2*np.pi, 91, endpoint=False)
	Y_q = micarray.modal.angular.sht_matrix(N, azi_pwd, np.pi/2)
	# get radial filters
	bn = micarray.modal.radial.spherical_pw(N, k, r, setup='rigid')
	dn, _ = micarray.modal.radial.regularize(1/bn, 100, 'softclip')
	D = micarray.modal.radial.diagonal_mode_mat(dn)
	# compute the PWD
	A_mb = np.matmul(np.matmul(np.conj(Y_q.T), D), Y_p)
	q_mb = np.squeeze(np.matmul(A_mb, np.expand_dims(p, 2)))
	q_mb_t = np.fft.fftshift(np.fft.irfft(q_mb, axis=0), axes=0)

In the following part (PWD), however, Y_p is used for the analysis (forward transform).

The user should also note that Y_p is different for the analysis (with weights)

Y_p = micarray.modal.angular.sht_matrix(N, azi, elev, weights)

and synthesis (without weights)

Y_p = micarray.modal.angular.sht_matrix(N, azi, elev)

The same applies to CHT.

In my opinion, we need to specify the usage of the SHT and CHT matrices in the document or in the docstrings.

[narahahn]

I just changed the title of this issue, since it's about the usage of the matrices. There is nothing wrong with the matrices themselves.

There seem to be (at least) two issues here:

1. The example is clearly inconsistent.
2. What convention shall we use for transform matrices?

Currently the output sht_matrix is meant to be interpreted as a forward transform (analysis) matrix.
This seems to be common for DFT matrices, e.g. scipy.linalg.dft or Matlab's dft_matrix.
Under this convention, the generation of p in the first example is wrong.

But it is maybe better to treat sht_matrix as the inverse transform (synthesis) matrix?
One could argue that the "basis functions" are "more basic" than the coefficient functionals.
Should we change to this convention?

Either way it really needs better documentation.

@narahahn and @spors agreed to use the convention of Rafaely, 2015:
The returns of sht_matrix and cht_matrix shall be interpreted as inverse transform (synthesis) matrices from now on.

Turns out the docstrings of sht_matrix and cht_matrix are inconsistent, too.
:(