Hi everyone,

In my work I do a lot with Slepian functions, and I find scipy's current capability to make those functions very limited, so I propose to expand it.

I'll start with a little motivation for why I think Slepians are great and why it's worth to expand scipy's Slepian capabilities, and then I propose the changes I would make and the new code I would add.

1. Motivation

Slepians (or, using the official term: "discrete, prolate, spheroidal sequences (DPSS)") are the provably optimal solution to the so-called spectral concentration problem [1], so they have a lot of application in frequency analysis [2] and recently my group and I started working on them in the context of ambient noise analysis for quantum computers [3, 4].

Slepians are the eigensolutions of the so-called Toelpitz matrix and there are therefor many of them. The different eigenvectors are called "orders" of the Slepians. In frequency analysis and in my work, it is very common to use multiple orders because (given the right bandwidth parameter) a set of Slepian orders forms a maximally spectrally concentrated "rectangle" in frequency space. Based on those Slepian orders, Thompson [2] developed the "multitaper" spectral reconstruction technique, in which multiple orders are applied as window functions and the optimal guess for the frequency spectrum is obtained by combining individual order estimates.

2. What's the current implementation

The current implementation of Slepians in scipy's signal module only offers the zeroth order, so it is not possible to use for e.g. the multitaper spectral analysis. In addition to that, because of the way the eigensystem is solved, it breaks with a MemoryError when using high numbers of points (i.e. a few thousand). The request for high numbers of points is not unreasonable: people do take long time-traces and the FFT makes it possible to compute Fourier transforms on those.

3. How could it be expanded?

I suggest to re-implement the existing function in order to:

a) calculate a variable number of Slepian orders

b) allow for higher point numbers by using a Newtonian solver for the eigenvalues

I have working code for this, which is based on the implementation suggested in [5] but is in pure Python. It consists of three functions: the main Slepian function plus two little helper functions which would not be exposed to the user. Most of the code is linted and pep8'd but I haven't got unit tests yet.

Once we have the full set of Slepians available, from there we (or I) could also implement the multitaper spectral analysis technique, but let's do one step at a time.

On a side node, not sure if this is a good argument for expanding the library, but both a Slepian function and a multitaper function exist in Matlab as well.

4. What's next?

So, what do you guys think of my proposal? If you like it, please advise me on how to proceed as this is the first time I'm planning to contribute. From what I took of the web page, I'll start by downloading the repo - and then should I make a new branch? I'm also happy to send some code through via email for your review first.

Thanks, and I'm looking forward to hearing back from you,

Virginia (from sunny Sydney)

----

References:

[1] https://en.wikipedia.org/wiki/Spectral_concentration_problem

[2] Thompson, "Spectrum estimation and harmonic analysis". Proc. IEEE 70(9), 1982 DOI: 10.1109/PROC.1982.12433

[3] V. Frey et al. "Application of optimal band-limited control protocols to quantum noise sensing", Nature Comms 8, 2189, 2017. DOI: 10.1038/s41467-017-02298-2

[4] L.M. Norris et al. "Optimally band-limited spectroscopy of control noise using a qubit sensor", arXiv preprint at https://arxiv.org/abs/1803.05538

[5] W.H. Press "Numerical Recipes in C++: The Art of Scientific Computing" Cambrige University Press, 3rd Edition