博客的公式排版非常不好看，可以移步这里阅读pdf版。

如有错误敬请指正！

图形学中PRT使用球谐函数来预计算光照信息。

网上找到的球谐函数的相关数学原理比较凌乱，有些讲的也不是很清楚，所以自己以数学专业的角度整理了一下。

用到的数学知识：

• 微积分（数学分析）
• 线性代数
• 复变函数
• 泛函分析（较少）

注意：下文中的单个字母上标，如果不是标在括号上，一律不代表“幂次”。

内积空间

考虑空间L^2 (S^2)，即，在三维空间单位球表面上平方可积（模的平方在球面S^2上的积分存在且有限）的函数全体。在该空间上定义一个（复）内积（类比于向量的点乘、内积）
\langle f,\ g \rangle := \int_{S^2} f(w)\overline{g(w)}\,{\rm d} w.
这样的话，\forall f\in L^2(S^2)，
\Vert f\Vert^2:=\langle f,f\rangle = \int_{S^2} f(w)\overline{f(w)}\,{\rm d} w\ge 0.
即定义了一个范数\Vert\cdot \Vert（类比于向量的长度）。可以注意到，如上定义的内积总是有这样一些性质：
\begin{aligned}
\langle f,g\rangle &= \overline{\langle g,f\rangle}\\
k\langle f,g\rangle &= \langle kf,g\rangle = \langle f,\overline k g\rangle,\ \forall k\in\mathbb C\\
\langle f_1+f_2,g\rangle &= \langle f_1,g\rangle + \langle f_2,g\rangle
\end{aligned}
类比一般的线性空间，我们想找到一组正交系（正交基），这样就可以引入类似于坐标这样的概念。内积空间中的一组规范正交系，记为(e_i) _ { i \in \mathbb N}，需要满足：

\langle e_i, e_j\rangle = \delta _{ij},
其中
\begin{aligned}
\delta_{ij}=\begin{cases}
1,& i=j,\\
0,&i\neq j.
\end{cases}
\end{aligned}

连带勒让德多项式

连带（也称伴随）勒让德多项式是一组在L^2[-1,1]上的正交系。连带勒让德多项式记为P_l^m(x)，其中0\le m\le l，m,l\in\mathbb N。它没有简单的显式表达，但可以递推地求出。它满足
\int_{-1}^1 P_l^m(x)P_k^m(x) \,{\rm d}x= \frac2{2l+1}\frac{(l+m)!}{(l-m)!}\delta_{kl},
即，在m固定时，P_l^m与P_k^m正交，其中l\neq k。

对于m<0的情形，推广定义：
P_l^{-m}(x) = (-1)^m \frac{(l-m)!}{(l+m)!} P_l^m(x).

### 球谐函数（复数形式）

定义球谐函数如下：
Y_l^m (\theta,\varphi) = N_l^m P_l^{m}(\cos \theta){\rm e}^{{\rm i}m\varphi}, \ \ \ -l\le m\le l,
其中N_l^m是归一化常数，即，N_l^m是使得Y_l^m与自身的内积为1（范数为1）的常数（总是实数）：
(N_l^m)^{-2} = \langle P_l^m(\cos \theta){\rm e}^{{\rm i}m\varphi},\ P_l^m(\cos \theta){\rm e}^{{\rm i}m\varphi}\rangle.
考虑内积：
\langle Y_l^m ,Y_k^n\rangle ,
展开计算（注意到\overline{{\rm e}^{{\rm i}n\varphi}}={\rm e}^{-{\rm i}n\varphi}）：
\langle Y_l^m ,Y_k^n\rangle=N_l^m N_k^n\int_{S^2} P_l^m(\cos \theta)P_k^n(\cos\theta){\rm e}^{{\rm i}(m-n)\varphi}\,{\rm d} w.
拆开\theta和\varphi相关的项
\langle Y_l^m ,Y_k^n\rangle = N_l^m N_k^n \int_{0}^{2\pi} {\rm e}^{{\rm i}(m-n)\varphi}\,{\rm d}\varphi \int_0^\pi P_l^m(\cos\theta)P_k^n(\cos\theta)\sin\theta\,{\rm d}\theta.

• 如果m\neq n，那么
  \int_{0}^{2\pi} {\rm e}^{{\rm i}(m-n)\varphi}\,{\rm d}\varphi = \int_0^{2\pi} \cos(m-n)\,{\rm d}\varphi+{\rm i}\int_0^{2\pi} \sin(m-n)\,{\rm d}\varphi = 0.
  m=n时被积项退化为1，积分值为2\pi。于是
  \langle Y_l^m ,Y_k^n\rangle = 2\pi \delta_{mn} N_l^m N_k^n \int_0^\pi P_l^m(\cos\theta)P_k^n(\cos\theta)\sin\theta\,{\rm d}\theta.

• 现在设m=n，由于连带勒让德函数的性质以及t=\cos\theta换元可以立即得出
  \langle Y_l^m ,Y_k^m\rangle = 2\pi N_l^m N_k^m \delta_{lk}.

通过以上讨论可以知道，当且仅当l=k且m=n时上述内积非0。又由N_l^m定义，我们得到
\langle Y_l^m ,Y_k^n\rangle = \delta_{lk}\delta_{mn},
从而全体复球谐函数成为一组规范正交系。

球谐函数直角坐标表示

引理. Y_l^m(x,y,z)是关于x,y,z的l次多项式，可以写成关于z的l-|m|次多项式与关于x,y的|m|次齐次多项式(x\pm{\rm i}y)^{|m|}的乘积。

Proof. 先考虑m\ge 0的情形。
Y_l^m(\theta,\varphi) =N_l^m P_l^{m}(\cos \theta){\rm e}^{{\rm i}m\varphi},
根据递推公式
P_l^m(t) = (-1)^m (1-t^2)^{m/2} \frac{{\rm d}^m}{{\rm d}t^m}P_l(t),
其中P_l(t)为一般的勒让德多项式。于是
P_l^m(\cos \theta) = (-1)^m \sin^m \theta F(\cos \theta),
其中
F(t) :=\frac{{\rm d}^m}{{\rm d}t^m}P_l(t).
因为勒让德多项式
P_l(t) =\frac{1}{2^l l!}\frac{{\rm d}^l}{{\rm d}t^l}((t^2-1)^l)
是一个关于t的l次多项式，从而F(t)是一个关于t的l-m次多项式。因为在S^2上z=\cos\theta，所以F(\cos\theta)是关于z的l-m次多项式。于是
\begin{aligned}
Y_l^m(\theta,\varphi) &=N_l^m (-1)^m F(z) \sin^m \theta {\rm e}^{{\rm i}m\varphi}\\
&=N_l^m (-1)^m F(z) (\sin \theta {\rm e}^{{\rm i}\varphi})^m\\
&=N_l^m (-1)^m F(z) (\sin \theta \cos\varphi+{\rm i}\sin\theta \sin \varphi)^m\\
& = N_l^m (-1)^m F(z) (x+{\rm i}y)^m,
\end{aligned}
这就说明了对于m\ge 0的情况成立。对于m<0的情况，与-m同理可得，结果可写为kG(z)(x-{\rm i}y)^{|m|}。证毕。

无穷维坐标与帕塞瓦尔公式

假设我们有一组希尔伯特（范数完备的内积空间）空间X上的规范正交系(e_i) _ {i\ge 0}。对于任意一个点f \in X，我们可以将f写为如下形式（定义与为什么可以这样做的证明稍繁琐，可以参考任意泛函分析教材，不详细展开了）：
f = \sum_{j=0}^\infty a_j e_j,
写成无穷维坐标形式即(a_0, a_1, a_2, \dots). 利用规范正交系的性质我们可以得到
\langle f, e_i\rangle = \sum_{j=0}^\infty a_j\langle e_j, e_i\rangle = a_i,
以及内积空间上的勾股定理的推广帕塞瓦尔等式
\Vert f\Vert^2 = \langle f, f\rangle = \sum_{j=0}^\infty \langle a_je_j, a_je_j\rangle = \sum_{j=0}^\infty |a_j|^2.

对于球谐函数系，即X = L^2(S^2)，我们可以写成
f = \sum_{l=0}^\infty \sum_{m=-l}^l a_l^m Y_l^m,
以及帕塞瓦尔等式
\Vert f\Vert^2 = \sum_{l=0}^\infty \sum_{m=-l}^l |a_l^m|^2.

球谐系数
a_l^m = \int f(w) \overline{Y_l^m (w)}\,{\rm d}w.

球谐函数（实数形式）

复数形式的球谐函数一般而言并不好用，就比如上面计算系数的积分，如果球谐函数是实值的，能省下很多麻烦。

我们用球谐函数的复数形式的Y_{1}^1, Y_1^0,Y_1^{-1}举例如何构造实值球谐函数。前面我们已经证明了它们满足规范正交系的性质，我们任取Y_l^m，l\neq 1。这时候一定有
\begin{aligned}
\langle Y_l^m, Y_1^1\rangle &= 0,\\
\langle Y_l^m, Y_1^0\rangle &= 0,\\
\langle Y_l^m, Y_1^{-1}\rangle &= 0.
\end{aligned}
也就是说，对于函数f，如果满足
f = c_{-1}Y_1^{-1}+c_0Y_1^0 +c_1 Y_1^1,
那么
\langle Y_l^m ,f\rangle = 0.
如果我们从全体可以写为c_{-1}Y_1^{-1}+c_0Y_1^0 +c_1 Y_1^1的函数（即Y_1^m张成的线性子空间）中，重新选出三个相互正交的、实值且范数为1的函数来替代掉原本的三个Y_1^m，这样构造的“球谐函数族”仍然是规范正交基，虽然它可能不再是物理或数学上普遍使用的球谐函数。

实际上这样做的时候，我们总是把Y_l^m和Y_l^{-m}这一对函数按照上面的方法重新选择。注意到Y_l^m与Y_l^{-m}形式非常相似，我们总是可以将它们两个作为一对派生出新的另一对，比如相加之后除以\sqrt2再取虚（实）部，相减之后除以\sqrt 2再取实（虚）部。另外，Y_l^0总是实值的，不需要重新选取。以及，-Y_l^m总是可以替换掉Y_l^m。

上面的讨论说明了实值球谐函数的选取方式非常多，主要看实际应用中你想如何定义它。实值球谐函数的例子，可以参考en wiki的实值球谐函数表。

投影

对于任意一个定义在球面S^2上的函数f，我们有分解
f = \sum_{l=0}^\infty \sum_{m=-l}^l a_l^m Y_l^m,
如果我们只取其中的前几层，比如l\le 2，就得到了f对子空间的投影
f_p = \sum_{l=0}^2 \sum_{m=-l}^l a_l^m Y_l^m.
f_p是在l\le 2的所有球谐函数张成的空间里，最接近 f的那个。因为根据帕塞瓦尔等式
\Vert f-f_p\Vert^2 = \sum_{l=3}^\infty \sum_{m=-l}^l |a_l^m|^2,
对于l\le 2球谐函数张成的子空间，这一差值的范数至少是\Vert f-f_p\Vert^2.可以注意到，如果层数l取值越大，f_p对f的逼近就越好。

在图形学中我们一般用投影f_p来代替原本的f，这样做一是因为我们不可能算出无穷个球谐系数a_l^m，二是取前几层得到的投影在大多数情况下，已经非常近似于原本的f。

球谐函数的旋转

根据施密特正交化，我们总是能轻松构造出一组规范正交基，除了球谐函数，当然也可以选择别的规范正交基来做分解。使用球谐函数是因为它们的另一个优秀性质：快速旋转。

考虑一个旋转R。f与旋转R的复合，仍然是一个定义在S^2的函数。

对于Y_l^m与R的复合，必然也有
Y_l^m\circ R(\theta,\varphi) = \sum_{k=0}^\infty \sum_{n=-k}^k b_{l,m,k}^n(R) Y_k^n(\theta, \varphi).
但是，对于球谐函数而言，实际上可以写成
Y_l^m\circ R(\theta,\varphi) = \sum_{n=-l}^l b_{mn} Y_l^n(\theta, \varphi),
其中b_{mn}和R有关。也就是说，第l层的球谐函数与任意旋转的复合，仍然在第l层球谐函数张成的子空间内。这一性质与三维旋转群SO(3)的不可约表示密不可分，严格的证明可以参考物理学中的群论（第2版）第122页。

至于如何解出上述的a_l^n，实际上我们可以把它看成解一个2l+1元线性方程组，那么我们自然需要2l+1个方程。取2l+1个不同的球面上的采样点(\theta_i, \varphi_i),\ i=-l,\dots, l，得到2l+1个等式：（有的说法是将采样点投影为球谐系数，这一说法是非常不严谨的。球谐函数分解的是L^2(S^2)上的函数，而不是S^2上的点）
Y_l^m\circ R(\theta_i, \varphi_i) = \sum_{n=-l}^l b_{mn} Y_l^n(\theta_i, \varphi_i),\ i=-l,\dots, l.
记矩阵A = \left(Y_l^n(\theta_i, \varphi_i)\right) _ {i,n}，向量b_{m} = (b_{m,-l},\ \dots,\ b_{ml})^{\rm T}，向量c_m = (Y_l^m\circ R(\theta_{-l}, \varphi_{-l}), \dots, Y_l^m\circ R(\theta_{l}, \varphi_{l}))^{\rm T}。则得到方程组
Ab_m = c_m.
于是（若A可逆）
b_m = A^{-1} c_m.
这个等式只解出了Y_l^m的旋转。现在假设我们待旋转的f在第l层上的球谐系数是a_l^m，这一层的投影f_p可以分解为
f_p= \sum_{m=-l}^l a_l^m Y_l^m,
则
f_p\circ R = \sum_{m=-l}^l a_l^m \sum_{n=-l}^l b_{mn} Y_l^n = \sum_{n=-l}^l Y_l^n \sum_{m=-l}^l a_l^m b_{mn}.
把所有列向量b_m合并，写成矩阵形式B = (b_m) _ {m=-l}^l，同理C=(c_m) _ {m=-l}^l，得到
B = A^{-1}C
Y_l^n项新的系数是\sum_{m=-l}^l a_l^m b_{mn}，写成一列，即为
Ba_l = A^{-1}C a_l.
实际应用中，矩阵A的值是可以预先算好的，因为A与输入的旋转R没有关联，也就是说，只要选取的2l+1个采样点，能够让算出来的A是可逆矩阵，那么这些采样点就是良好的。我们可以直接计算出A^{-1}，需要进行旋转时，只需要用R计算出C，即可得到旋转后的球谐系数。

如果你在别处见到的结果是写成SA^{-1}的，而不是A^{-1}C，那是因为，在上面的讨论中采用的A，相对于那种写法是转置过的，这样的目的在于，最终计算时用的是矩阵左乘列向量（数学上比较通用的写法），而不是矩阵右乘行向量。

其他分解方式

球面上当然也可以选取其他的方式分解，比如Zonal Harmonics以及Wavelets. Zonal Harmonics是一些特殊的球谐函数旋转后的结果。Wavelets不适合快速旋转，但它应用于高频信息的保留时效果比球谐基好。