[SciPy-Dev] Dev suggestion: expand signal.slepian functions

[Virginia Frey]

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 
<https://doi.org/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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20180408/c68d1eb7/attachment.html>