A function $f:\mathbb{R}\rightarrow\mathbb{R}$ is called band-limited if there exists $c>0$ and a function $\sigma\in L^2[-1, 1]$ s.t. $ f(x) = \int_{-1}^1e^{icxt}\sigma(t)dt. $ This is motivated by the inverse Fourier-transform. Let $ f(x) = \int_{-\infty}^\infty \hat{f}(k)e^{i kx}dk $ be the representation of a function $f$ via its Fourier transform $\hat{f}$. The function $f$ is band-limited if there exists $c>0$ s.t. $ f(x) = \int_{-c}^c e^{ikx}\hat{f}(k)dk $ A change of variable gives $ f(x) = \int_{-1}^1e^{icxt}\sigma(t)dt $ with $\sigma(t) = c\hat{f}(ct)$. Hence, the frequency range of a band-limited function spans a range of $2c$. Now define the operator $F_c: L^2[-1, 1]\rightarrow L^2[-1, 1]$ by $ F_c(\phi)(x) = \int_{-1}^1 e^{icxt}\phi(t)dt. $ The eigenfunctions $\psi_n$ of this operator are called *Prolate Spheroidal Wave Functions* (PSWFs). We denote the associated eigenvalues with $\lambda_n$ with the convention that $|\lambda_0|\geq |\lambda_1| \geq \dots$. The following facts hold: - For $c> 0$ the functions $\psi_j$ are purely real, orthonormal and complete in $L^2[-1, 1]$. - Even-numbered eigenfunctions are even. - Odd-numbered eigenfunctions are odd. - All eigenvalues are nonzero and simple. - Even-numbered eigenvalues are purely real. - Odd-numbered eigenvalues lie on the imaginary axis. ## Numerical evaluation The idea is to express PSWFs as a series of [[Legendre Polynomials]] $ \psi_n(x) = \sum_{j=0}^\infty \beta_j^{(n)} \overline{P}_j(x). $ Here, the $\overline{P}_j(x)$ are normalised Legendre polynomials with $\overline{P}_j = \sqrt{j + 1/2}P$ and $\|\overline{P}_j\|_{L^2[-1, 1]}$ =1. The coefficients $\beta_j^{(n)}$ are defined by $ \beta_j^{(n)} = \int_{-1}^1 \psi_n(x)\overline{P}_j(x)dx. $ It follows immediately from the even/odd-ness of the PSWFs and of the Legendre Polynomials that for all even-numbered PSWFs the odd coefficients are zero and that for all odd-numbered PSWFs the even-numbered coefficients are zero. To compute the coefficients we define the infinite matrix eigenvalue problem $ A\beta^{(n)} = \chi_n \beta^{(n)} $ with $A$ defined as follows: - $A[k, k] = k(k+1) + \frac{2k(k+1) - 1}{(2k+3)(2k-1)}c^2$ - $A[k, k + 2] = \frac{(k + 2)(k+1)}{(2k+3)\sqrt{(2k+1)(2k+5)}}c^2$ - $A[k+2, k] = A[k, k+2]$. This matrix is separable into two matrices $A_{even}$ and $A_{odd}$ by selecting only even/odd rows and columns, giving two tridiagonal symmetric matrix eigenvalue problems of the form $ \begin{align} A_{even}\begin{bmatrix}\beta_0^{(\tilde{n})}\\ \beta_2^{(\tilde{n})}\\ \vdots\end{bmatrix} &= \chi_\tilde{n}^{(even)} \begin{bmatrix}\beta_0^{(\tilde{n})}\\ \beta_2^{(\tilde{n})}\\ \vdots\end{bmatrix},\nonumber\\ A_{odd}\begin{bmatrix}\beta_1^{(\tilde{n})}\\ \beta_3^{(\tilde{n})}\\ \vdots\end{bmatrix} &= \chi_\tilde{n}^{(odd)} \begin{bmatrix}\beta_1^{(\tilde{n})}\\ \beta_3^{(\tilde{n})}\\ \vdots\end{bmatrix},\nonumber\\ \end{align} $ with $n = 2\tilde{n}$ for $n$ even and $n = 2\tilde{n}+1$ for $n odd. To compute the PSWFs in practice we simply restrict these infinite matrix problems to finite sections and numerically compute the coefficients. If $n$ is even we use $A_{even}$ and if $n$ is odd we use $A_{odd}$ to compute half the coefficients with the other half of the coefficients being zero. Once we have computed the coefficient vector $\beta^{(n)}$ it is useful to absorb the normalisation factor of the $\overline{P}_n$ into $\beta^{(n)}$ via $\alpha_j^{(n)} := \sqrt{j + 1/2}\beta_j^{(n)}$. We can then use the [[[Clenshaw Algorithm]]] to evaluate the associated Legendre series. The eigenvalues $\lambda_n$ can be computed from the equality $ \lambda_n\psi_n(0) = \int_{-1}^1 \psi_n(t)dt $ obtained from evaluating at $x=0$ the definition of PSWFs as eigenfunctions of the band-limited Fourier operator. A Python code to compute the Legendre coefficients $\alpha$ is given below. ```Python def pswf(pswf_order, m, c): """Compute the first m coefficients of the Legendre expansion of prolate spheroidal wave functions. - pswf_order: The order of the PSWF - m: The number of coefficients - c: The band-limit parameter """ ## If m is odd we increase by 1 as working with an even number of coefficients is easier. if m % 2 == 0: ncoeffs = m else: ncoeffs = m + 1 ## This is necessary to select the right rows/columns of the infinite matrix eigenvalue problem. ## For even order PSWFs we select the even rows/cols. For odd ones the odd rows/cols. if pswf_order % 2 == 0: increment = 0 else: increment = 1 ## The finite dimensional tridiagonal matrix whose eigenvectors we will compute. A = np.zeros([ncoeffs // 2 , ncoeffs //2], dtype='float64') for index in range(ncoeffs // 2 ): ## k denotes the right row in the infinite matrix problem. k = 2 * index + increment A[index, index] = k * (k + 1) + (2 * k * (k + 1) - 1) / ((2 * k + 3) * (2 * k - 1)) * c**2 if index + 1 < n // 2: A[index, index + 1] = (k + 2) * (k + 1) / ((2 * k + 3) * np.sqrt((2 * k + 1) * (2 * k + 5))) * c**2 A[index + 1, index] = A[index, index + 1] ## Now compute the eigenvalues and eigenvectors of A. w, v = scipy.linalg.eigh(A) beta = np.zeros(ncoeffs, dtype='float64') ## We have to assign the eigenvector to the right positions in the ## Legendre coefficients vector. ## We integer divide pswf_order by two because orders 0, 2, 4,... are mapped ## to the eigenvector numbers in the even problem and ## 1, 3, 5, ... are mapped to the eigenvector numbers in the odd problem. if pswf_order % 2 == 0: # Even case beta[::2] = v[:, pswf_order // 2] else: # Odd case beta[1::2] = v[:, pswf_order // 2] ## Fix the normalisation. We need the alpha coefficients of unnormalised ## Legendre polynomials. alpha = beta * np.sqrt(np.arange(ncoeffs) + .5) ## Finally, if we computed one coefficient too many we chop it off. if m % 2 == 0: return alpha else: return alpha[:-1] ``` Below is a plot of the fourth (n=3) PSWF for $c=10\pi$. ![[pswf.png]] ## Asymptotics For fixed $n$ and large $c$ [[#^c1443f| Xiao and Rokhlin]] derived asymptotic expansions of the form $ \psi_n^c(x) \sim e^{-cx^2/2}\sum_{j=0}^k(\alpha_j H_{n+4j}(\sqrt{c}x) + \beta_j H_{n-4j}(\sqrt{c}x)). $ Here, $H_m$ denotes a Hermite polynomial of order $m$. # References - [Wikipedia](https://en.wikipedia.org/wiki/Prolate_spheroidal_wave_function) - [Chebfun page on PSWFs](https://www.chebfun.org/examples/approx/Prolate.html) - H. Xiao, V. Rokhlin, and N. Yarvin, *Prolate spheroidal wave functions, quadrature and interpolation*, Inverse Problems 17, pp. 805--828 (2001), [doi: 10.1088/0266-5611/17/4/315](https://doi.org/10.1088/0266-5611/17/4/315) - A. Osipov, and V. Rokhlin, *On the evaluation of prolate spheroidal wave functions and associated quadrature rules*, Applied and Computational Harmonic Analysis 36, pp. 108--142 (2014), [doi: 10.1016/j.acha.2013.04.002](https://doi.org/10.1016/j.acha.2013.04.002) - H.Xiao, and V. Rokhlin, *High-Frequency Asymptotic Expansions for Certain Prolate Spheroidal Wave Functions*, The Journal of Fourier Analysis And Applications 9, pp. 575--596 (2003), [doi: 10.1007/s00041-003-0906-z](https://doi.org/10.1007/s00041-003-0906-z)] ^c1443f