Slepian function

Slepian functions are a class of spatio-spectrally concentrated functions^[1][2] that form an orthogonal basis for bandlimited or spacelimited spaces.^[3][4][5][6] That is, they are concentrated in space or time while spectrally bandlimited, or concentrated in spectral band while space- or time-limited. They are widely used as basis functions for constructive approximation^[7] and in linear inverse problems,^[8][9] and as apodization tapers or window functions^[10] in quadratic problems^[11] of spectral density estimation.^[12][13]

Slepian function constructions exist in discrete (regular^[14] and irregular^[15]) and continuous^[16][17][18] varieties, in one, two, and three dimensions,^[19] in Cartesian and spherical geometry, on surfaces and in volumes,^[20] on graphs,^[21] and in scalar, vector,^[22] and tensor forms.^[23]

Without reference to any of these particularities,^[24] let f {\displaystyle f} be a square-integrable function of physical space, and let H {\displaystyle {\mathcal {H}}} represent Fourier transformation, such that F = H f {\displaystyle F={\mathcal {H}}f} and H − 1 F = f {\displaystyle {\mathcal {H}}^{-1}F=f} . Let the operators R {\displaystyle {\mathcal {R}}} and L {\displaystyle {\mathcal {L}}} project onto the space of spacelimited functions, S R {\displaystyle {\mathcal {S}}_{R}} , and the space of bandlimited functions, S L {\displaystyle {\mathcal {S}}_{L}} , respectively, whereby R {\displaystyle R} is an arbitrary nontrivial subregion of all of physical space, and L {\displaystyle L} an arbitrary nontrivial subregion of spectral space. Thus, the operator R {\displaystyle {\mathcal {R}}} acts to spacelimit, and the operator H − 1 L H {\displaystyle {\mathcal {H}}^{-1}{\mathcal {L}}{\mathcal {H}}} acts to bandlimit the function f {\displaystyle f} .

Slepian's quadratic spectral concentration problem aims to maximize the concentration of spectral power to a target region L {\displaystyle L} , for a function that is spatially limited to a target region R {\displaystyle R} . Conversely, Slepian's spatial concentration problem maximizes the spatial concentration to R {\displaystyle R} of a function bandlimited to L {\displaystyle L} . Using ⟨ ⋅ , ⋅ ⟩ {\displaystyle \langle \cdot ,\cdot \rangle } for the inner product both in the space and the spectral domains, both problems are stated equivalently using Rayleigh quotients in the form

λ = ⟨ R H − 1 L F , R H − 1 L F ⟩ ⟨ H − 1 L F , H − 1 L F ⟩ = ⟨ L H R f , L H R f ⟩ ⟨ H R f , H R f ⟩ = maximum . {\displaystyle \lambda ={\frac {\langle {\mathcal {R}}{\mathcal {H}}^{-1}{\mathcal {L}}\,F,{\mathcal {R}}{\mathcal {H}}^{-1}{\mathcal {L}}\,F\,\rangle }{\langle {\mathcal {H}}^{-1}{\mathcal {L}}\,F,{\mathcal {H}}^{-1}{\mathcal {L}}\,F\,\rangle }}={\frac {\langle {\mathcal {L}}{\mathcal {H}}{\mathcal {R}}f,{\mathcal {L}}{\mathcal {H}}{\mathcal {R}}f\rangle }{\langle {\mathcal {H}}{\mathcal {R}}f,{\mathcal {H}}{\mathcal {R}}f\rangle }}={\mbox{maximum}}.}

The equivalent spectral-domain and spatial-domain eigenvalue equations are

( L H R H − 1 L ) ( L F ) = λ ( L F ) {\displaystyle ({\mathcal {L}}{\mathcal {H}}{\mathcal {R}}{\mathcal {H}}^{-1\!}{\mathcal {L}})({\mathcal {L}}\,F\,)=\lambda ({\mathcal {L}}\,F\,)} and ( R H − 1 L H R ) ( R f ) = λ ( R f ) , {\displaystyle ({\mathcal {R}}{\mathcal {H}}^{-1\!}{\mathcal {L}}{\mathcal {H}}{\mathcal {R}})({\mathcal {R}}f)=\lambda ({\mathcal {R}}f),}

given that H {\displaystyle {\mathcal {H}}} and H − 1 {\displaystyle {\mathcal {H}}^{-1}} are each others' adjoints, and that R {\displaystyle {\mathcal {R}}} and L {\displaystyle {\mathcal {L}}} are self-adjoint and idempotent.

The Slepian functions are solutions to either of these types of equations with positive-definite kernels, that is, they are bandlimited functions G = L F {\displaystyle G={\mathcal {L}}F} , concentrated to the spatial domain within R {\displaystyle R} , or spacelimited functions of the form h = R f {\displaystyle h={\mathcal {R}}f} , concentrated to the spectral domain within L {\displaystyle L} .

Let g ( t ) {\displaystyle g(t)} and its Fourier transform G ( ω ) {\displaystyle G(\omega )} be strictly bandlimited in angular frequency between [ − W , W ] {\displaystyle [-W,W]} . Attempting to concentrate g ( t ) {\displaystyle g(t)} in the time domain, to be contained within the time interval [ − T , T ] {\displaystyle [-T,T]} , amounts to maximizing

λ = ∫ − T T g 2 ( t ) d t ∫ − ∞ ∞ g 2 ( t ) d t = maximum , {\displaystyle \lambda ={\frac {\int _{-T}^{T}g^{2}(t)\,dt}{\int _{-\infty }^{\infty }g^{2}(t)\,dt}}={\text{maximum}},}

which is equivalent to solving either, in the frequency domain, the convolutional integral eigenvalue (Fredholm) equation

∫ − W W D T ( ω , ω ′ ) G ( ω ′ ) d ω ′ = λ G ( ω ) , D T ( ω , ω ′ ) = sin ⁡ T ( ω − ω ′ ) π ( ω − ω ′ ) , | ω | ≤ W , {\displaystyle \int _{-W}^{W}D_{T}(\omega ,\omega ')G(\omega ')\,d\omega '=\lambda G(\omega ),\qquad D_{T}(\omega ,\omega ')={\frac {\sin T(\omega -\omega ')}{\pi (\omega -\omega ')}},\qquad |\omega |\leq W,}

or the time- or space-domain version

∫ − T T D W ( t , t ′ ) g ( t ′ ) d t ′ = λ g ( t ) , D W ( t , t ′ ) = sin ⁡ W ( t − t ′ ) π ( t − t ′ ) = ( 2 π ) − 1 ∫ − W W e i ω ( t − t ′ ) d ω , t ∈ R . {\displaystyle \int _{-T}^{T}D_{W}(t,t')g(t')\,dt'=\lambda g(t),\qquad D_{W}(t,t')={\frac {\sin W(t-t')}{\pi (t-t')}}=(2\pi )^{-1}\int _{-W}^{W}e^{i\omega (t-t')}d\omega ,\qquad t\in \mathbb {R} .}

Either of these can be transformed and rescaled to the dimensionless

∫ − 1 1 D ( x , x ′ ) g ( x ′ ) d x ′ = λ g ( x ) , D ( x , x ′ ) = sin ⁡ T W ( x − x ′ ) π ( x − x ′ ) . {\displaystyle \int _{-1}^{1}D(x,x')g(x')\,dx'=\lambda g(x),\qquad D(x,x')={\frac {\sin TW(x-x')}{\pi (x-x')}}.}

The trace of the positive definite kernel is the sum of the infinite number of real and positive eigenvalues,

N = ∑ α = 1 ∞ λ α = ∫ − 1 1 D ( x , x ′ ) d x = 2 T W π , {\displaystyle N=\sum _{\alpha =1}^{\infty }\lambda _{\alpha }=\int _{-1}^{1}D(x,x')\,dx={\frac {2TW}{\pi }},}

that is, the area of the concentration domain in time-frequency space (a time-bandwidth product).

One-dimensional scalar Slepian functions or tapers^[25] are the workhorse of the Thomson multitaper method of spectral density estimation.

We use g ( x ) {\displaystyle g(\mathbf {x} )} and its Fourier transform G ( k ) {\displaystyle G(\mathbf {k} )} to denote a function that is strictly bandlimited to K {\displaystyle {\mathcal {K}}} , an arbitrary subregion of the spectral space of spatial wave vectors.^[26] Seeking to concentrate g ( x ) {\displaystyle g(\mathbf {x} )} into a finite spatial region R ∈ R 2 {\displaystyle R\in \mathbb {R} ^{2}} , of area A {\displaystyle A} , we must find the unknown functions for which

λ = ∫ R g 2 ( x ) d x ∫ − ∞ ∞ g 2 ( x ) d x = maximum . {\displaystyle \lambda ={\frac {\int _{R}g^{2}(\mathbf {x} )\,d\mathbf {x} }{\int _{-\infty }^{\infty }g^{2}(\mathbf {x} )\,d\mathbf {x} }}={\mbox{maximum}}.}

Maximizing this Rayleigh quotient requires solving the Fredholm integral equation

∫ K D R ( k , k ′ ) G ( k ′ ) d k ′ = λ G ( k ) , D R ( k , k ′ ) = ( 2 π ) − 2 ∫ R e i ( k ′ − k ) ⋅ x d x , k ∈ K . {\displaystyle \int _{\mathcal {K}}D_{R}(\mathbf {k} ,\mathbf {k} ')\,G(\mathbf {k} ')\,d\mathbf {k} '=\lambda G(\mathbf {k} ),\qquad D_{R}(\mathbf {k} ,\mathbf {k} ')=(2\pi )^{-2}\int _{R}e^{i(\mathbf {k} '-\mathbf {k} )\cdot \mathbf {x} }\,d\mathbf {x} ,\qquad \mathbf {k} \in {\mathcal {K}}.}

The corresponding problem in the spatial domain is

∫ R D K ( x , x ′ ) g ( x ′ ) d x ′ = λ g ( x ) , D K ( x , x ′ ) = ( 2 π ) − 2 ∫ K e i k ⋅ ( x − x ′ ) d k , x ∈ R 2 . {\displaystyle \int _{R}\!D_{\mathcal {K}}(\mathbf {x} ,\mathbf {x} ')\,g(\mathbf {x} ')\,d\mathbf {x} '=\lambda g(\mathbf {x} ),\qquad D_{\mathcal {K}}(\mathbf {x} ,\mathbf {x} ')=(2\pi )^{-2}\int _{\mathcal {K}}e^{i\mathbf {k} \cdot (\mathbf {x} -\mathbf {x} ')}\,d\mathbf {k} ,\qquad \mathbf {x} \in \mathbb {R} ^{2}.}

Concentration to the disk-shaped spectral band K = { k : ‖ k ‖ ≤ K } {\displaystyle {\mathcal {K}}=\{\mathbf {k} :\|\mathbf {k} \|\leq K\}} allows us to rewrite the spatial kernel as

D K ( x , x ′ ) = K J 1 ( K ‖ x − x ′ ‖ ) 2 π ‖ x − x ′ ‖ , {\displaystyle D_{\mathcal {K}}(\mathbf {x} ,\mathbf {x} ')={\frac {KJ_{1}(K\|\mathbf {x} -\mathbf {x} '\|)}{2\pi \|\mathbf {x} -\mathbf {x} '\|}},}

with J 1 {\displaystyle J_{1}} a Bessel function of the first kind, from which we may derive that

N = ∑ α = 1 ∞ λ α = ∫ R D K ( x , x ′ ) d x = K 2 A 4 π , {\displaystyle N=\sum _{\alpha =1}^{\infty }\lambda _{\alpha }=\int _{R}D_{\mathcal {K}}(\mathbf {x} ,\mathbf {x} ')\,d\mathbf {x} =K^{2}{\frac {A}{4\pi }},}

in other words, again the area of the concentration domain in space-frequency space (a space-bandwidth product).

We denote g ( r ^ ) {\displaystyle g(\mathbf {\hat {r}} )} a function on the unit sphere Ω {\displaystyle \Omega } and its spherical harmonic transform coefficient g l m {\displaystyle g_{lm}} at the degree l {\displaystyle l} and order m {\displaystyle m} , respectively,^[24] and we consider bandlimitation to spherical harmonic degree L {\displaystyle L} , that is, g ∈ S L {\displaystyle g\in {\mathcal {S}}_{L}} . Maximizing the quadratic energy ratio within the spatial subdomain R ⊂ Ω {\displaystyle R\subset \Omega } via

λ = ∫ R g 2 ( r ^ ) d Ω ∫ Ω g 2 ( r ^ ) d Ω = maximum {\displaystyle \lambda ={\frac {\int _{R}g^{2}(\mathbf {\hat {r}} )\,d\Omega }{\int _{\Omega }g^{2}(\mathbf {\hat {r}} )\,d\Omega }}={\mbox{maximum}}}

amounts in the spectral domain to solving the algebraic eigenvalue equation

∑ l ′ = 0 L ∑ m ′ = − l ′ l ′ D l m , l ′ m ′ g l ′ m ′ = λ g l m , D l m , l ′ m ′ = ∫ R Y l m ( r ^ ) Y l ′ m ′ ( r ^ ) d Ω {\displaystyle \sum _{l'=0}^{L}\sum _{m'=-l'}^{l'}D_{lm,l'm'}g_{l'm'}=\lambda g_{lm},\qquad D_{lm,l'm'}=\int _{R}Y_{lm}(\mathbf {\hat {r}} )Y_{l'm'}(\mathbf {\hat {r}} )\,d\Omega } ,

with Y l m {\displaystyle Y_{lm}} the spherical harmonic at degree l {\displaystyle l} and order m {\displaystyle m} . The equivalent spatial-domain equation, ∫ R D ( r ^ , r ^ ′ ) g ( r ^ ) d Ω = λ g ( r ^ ) , D ( r ^ , r ^ ′ ) = ∑ l = 0 L ∑ m = − l l Y l m ( r ^ ) Y l m ( r ^ ′ ) = ∑ l = 0 L ( 2 l + 1 4 π ) P l ( r ^ ⋅ r ^ ′ ) , {\displaystyle \int _{R}D(\mathbf {\hat {r}} ,\mathbf {\hat {r}} ')g(\mathbf {\hat {r}} )\,d\Omega =\lambda g(\mathbf {\hat {r}} ),\qquad D(\mathbf {\hat {r}} ,\mathbf {\hat {r}} ')=\sum _{l=0}^{L}\sum _{m=-l}^{l}Y_{lm}(\mathbf {\hat {r}} )Y_{lm}(\mathbf {\hat {r}} ')=\sum _{l=0}^{L}\left({\frac {2l+1}{4\pi }}\right)P_{l}(\mathbf {\hat {r}} \cdot \mathbf {\hat {r}} '),} is a homogeneous Fredholm integral equation of the second kind, with a finite-rank, symmetric, separable kernel.

The last equality is a consequence of the spherical harmonic addition theorem which involves P l {\displaystyle P_{l}} , the Legendre polynomial. The trace of this kernel is given by

N = ∑ α = 1 ( L + 1 ) 2 λ α = ∫ R D ( r ^ , r ^ ) d Ω = ∑ l = 0 L ∑ m = − l l D l m , l m = ( L + 1 ) 2 A 4 π , {\displaystyle N=\sum _{\alpha =1}^{(L+1)^{2}}\lambda _{\alpha }=\int _{R}D(\mathbf {\hat {r}} ,\mathbf {\hat {r}} )\,d\Omega =\sum _{l=0}^{L}\sum _{m=-l}^{l}D_{lm,lm}=(L+1)^{2}{\frac {A}{4\pi }},}

that is, once again a space-bandwidth product, of the dimension of S L {\displaystyle {\mathcal {S}}_{L}} and the fractional area of R {\displaystyle R} on the unit sphere Ω {\displaystyle \Omega } , namely A / ( 4 π ) {\displaystyle A/(4\pi )} .