  Classical prolate spheroidal functions play an important role in the study of
time-band limiting, scaling limits of random matrices, and the distribution of
the zeros of the Riemann zeta function. We establish an intrinsic relationship
between discrete-continuous bispectral functions and the prolate spheroidal
phenomenon. The former functions form a vast class, parametrized by an infinite
dimensional manifold, and are constructed by Darboux transformations from
classical bispectral functions associated to orthogonal polynomials. Special
cases include spherical functions. We prove that all such Darboux
transformations which are self-adjoint in a certain sense give rise to integral
operators possessing commuting differential operators and to discrete integral
operators possessing commuting shift operators. One particularly striking
implication of this is the correspondence between discrete and continuous pairs
of commuting operators. Moreover, all results are proved in the setting of
matrix valued functions, which provides further advantages for applications.
Our methods rely on the use of noncommutative matrix valued Fourier algebras
associated to discrete-continuous bispectral functions. We produce the
commuting differential and shift operators in a constructive way with explicit
upper bounds on their orders and bandwidths, which is illustrated with many
concrete examples.
