## 3 Proposed 3-D Sound Image Reproduction Based on Spherical Harmonic Expansion
### Algorithm of Proposed Method  
We explain the procedures of the proposed method to generate an accurate 3-D sound image. 1. Generation of a target directivity pattern \(D^{\text{Tar}}(\theta,\phi)\). The proposed method requires a target directivity pattern \(D^{\text{Tar}}(\theta,\phi)\) including a peak at the target sound image position. This directivity pattern affects the clarity of the 3-D sound image at the target position. In this paper, a provisional directivity pattern is generated by multiple signal classification (MUSIC) [16]. MUSIC estimates the direction of arrival of a sound source and generates the sharp spatial spectrum to the direction of the sound source. In Step 1, we generate the spatial spectrum \(P(\theta,\phi)\) for the target sound image at the position of \((r^{\text{Tar}},\theta^{\text{Tar}},\phi^{\text{Tar}})\) by using MUSIC. Then, the generated spatial spectrum \(P(\theta,\phi)\) is normalized to the target directivity pattern \(D^{\text{Tar}}(\theta,\phi)\) as: \[D^{\text{Tar}}(\theta,\phi)=\frac{P(\theta,\phi)-P(\theta^{\text{Min}},\phi^{ \text{Min}})}{P(\theta^{\text{Tar}},\phi^{\text{Tar}})-P(\theta^{\text{Min}}, \phi^{\text{Min}})},\] (8) where \((\theta^{\text{Min}},\phi^{\text{Min}})\) is the position with the smallest spatial spectrum. An example of the target directivity pattern \(D^{\text{Tar}}(\theta,\phi)\) is shown in Figure 8.
2. Calculation of the target directivity pattern in the SH domain \(A^{\text{Tar}}_{nm}\). We calculate the target directivity pattern in the SH domain \(A^{\text{Tar}}_{nm}\) using Equations (5) and (7) as: \[A^{\text{Tar}}_{nm}=\int_{0}^{\pi}\int_{0}^{2\pi}D^{\text{Tar}}(\theta,\phi)Y ^{m}_{n}(\theta,\phi)^{*}\sin\theta d\theta d\phi.\] (9)
3. Calculation of the weighting factor \(w_{i}(k)\). Figure 6: Spherical coordinate system. Figure 7: Shapes of SH function at order \(n=2\) and degree \(m=2\). We calculate the weighting factor \(w_{i}(k)\) for \(i\)th loudspeaker at the position of \((r_{i}^{\text{LS}},\varrho_{i}^{\text{LS}},\phi_{i}^{\text{LS}})\) to reproduce the target directivity pattern \(A_{nm}^{\text{Tar}}\) in the real sound field. The weighting factor \(w_{i}(k)\) is obtained by:  
\[w_{i}(k)=\sum_{n=0}^{\infty}\sum_{m=-n}^{n}w_{nm,i}(k)Y_{n}^{m}(\varrho_{i}^{ \text{LS}},\phi_{i}^{\text{LS}}), \tag{10}\]  
where \(k\) is the wave number and \(i\in\{1,2,\cdots,22\}\) is the loudspeaker index. On the basis of the mode matching, the target directivity pattern \(A_{nm}^{\text{Tar}}\) has a relationship between the weighting factor \(w_{nm,i}(k)\) in the SH domain as:  
\[w_{nm,i}(k) =\frac{A_{nm}^{\text{Tar}}}{b_{n}^{\text{S}}(k,r_{i}^{\text{LS}})}, \tag{11}\] \[b_{n}^{\text{S}}(k,r_{i}^{\text{LS}}) \triangleq\frac{4\pi}{(r_{i}^{\text{LS}})^{2}k}\frac{j_{n}(kr^{ \text{Tar}})}{j_{n}^{{}^{\prime}}(kr_{i}^{\text{LS}})}, \tag{12}\]  
where \(j_{n}(\cdot)\) and \(j_{n}^{{}^{\prime}}(\cdot)\) are the spherical Bessel function and its derivative, respectively, and \(b_{n}^{\text{S}}(\cdot)\) is the mode strength [17]. The mode strength \(b_{n}^{\text{S}}(\cdot)\) theoretically represents the radial strength of the directivity. Substituting Equation (11) into Equation (10), the weighting factor \(w_{i}(k)\) can be obtained as follows:  
\[w_{i}(k)=\sum_{n=0}^{\infty}\sum_{m=-n}^{n}\frac{A_{nm}^{\text{Tar}}}{b_{n}^{ \text{S}}(k,r_{i}^{\text{LS}})}Y_{n}^{m}(\varrho_{i}^{\text{LS}},\phi_{i}^{ \text{LS}}). \tag{13}\]
4. Generation of the input signals for 22 loudspeakers. The weight factor \(w_{i}(k)\) in the spatial domain can be used as the frequency domain filter, that is, \[W_{i}(\omega) =w_{i}\Big{(}\frac{\omega}{c}\Big{)},\] (14) Figure 8.0 (15) where \(\omega\) is the angular frequency and \(c\) is the speed of sound. Then, the input signal for \(i\)th loudspeaker \(y_{i}(t)\) is obtained by: Figure 8.1 (16) Figure 8.2 (17) where \(x(t)\) is the sound source, \(\text{IDFT}[\cdot]\) and \(\text{DTFT}[\cdot]\) represent the inverse discrete time Fourier transform and discrete time Fourier transform operators, respectively.