ENH: signal: implement signal.planck (Planck-taper window)

[wavexx]

This is my tentative implementation of the Planck-taper window (added as signal.planck). See http://arxiv.org/pdf/1003.2939

Window and resulting spectra match the paper as expected, although I need to double-check the boundary conditions for very low samples to check my understanding that the first and last sample are always 0. If nobody chimes in, I'll mail the author.

Meanwhile, anything missing?

[e-q]

Some tests should be added. You can take a look at what was was done for the Tukey window in gh-4584.

[wavexx]

I'll add some tests tomorrow.

[codecov-io]

@@            master   #6012   diff @@
======================================
  Files          238     238       
  Stmts        43822   43840    +18
  Branches      8215    8219     +4
  Methods          0       0       
======================================
+ Hit          34253   34269    +16
- Partial       2602    2604     +2
  Missed        6967    6967       

Review entire Coverage Diff as of 5c338ca

Powered by Codecov. Updated on successful CI builds.

[wavexx]

Should be good now. I indeed found a typo while cross-checking the formula again. If you look at the paper, I tried to reproduce the same plot here (see page 4):

With the following:

from scipy import signal
from scipy.fftpack import fft, fftshift, fftfreq
import matplotlib.pyplot as plt

smp = 2**16-1
ftlen = 2048**2
fig, axes = plt.subplots(1, 2, figsize=(12, 5.5))
fig.suptitle("Planck-taper window")

window = signal.boxcar(smp)
axes[0].plot(np.linspace(-1, 1, smp), window, label='square')
A = fft(window, ftlen) / (smp/2)
freq = fftshift(fftfreq(ftlen, 1/(smp/2)))
resp = fftshift(np.abs(A))
axes[1].plot(freq, resp, label='square')

for f in [100, 30, 10]:
    e = 1/f
    l = '$\epsilon = 1/{}$'.format(f)
    window = signal.planck(smp, e)
    axes[0].plot(np.linspace(-1, 1, smp), window, label=l)
    A = fft(window, ftlen) / (smp/2)
    freq = fftshift(fftfreq(ftlen, 1/(smp/2)))
    resp = fftshift(np.abs(A))
    axes[1].plot(freq, resp, label=l)

axes[0].legend(loc='lower left')
axes[0].axis([0.79, 1.01, -0.1, 1.1])
axes[0].set_xlabel('Time /s')
axes[1].set_xscale('log')
axes[1].set_yscale('log')
axes[1].axis([1, 1e4, 1e-10, 1])
axes[1].set_xlabel("Frequency /Hz")
axes[1].legend(loc='upper right')

I actually find this visualization much better suited to see the frequency-response compared to the current examples in other window functions. It's a bit more expensive to compute, but I'd love to have that for all functions.

I tried to use subplots inside the example itself, but somehow run into problems and gave up, as regenerating the docs takes some time...

[e-q]

By "this visualization", do you mean logarithmic x-axis? I definitely agree on that point. It could also be appropriate to show responses of related windows. (Though, if you were to undertake this for other windows, you should open a different PR)

For this window, I was curious how it compares to the Tukey window, which is similar in concept. Since the Planck window was designed for GW chirp signals, I wasn't sure how much utility it really would bring to the growing zoo of window functions already in scipy.

The authors in the reference you linked to only compare this window to the rectangular window, which is a bit disingenuous, since pretty much any tapered window would show less spectral leakage than the rectangular window, so it doesn't really argue for the merit of their particular design. In any case, I did the comparison to Tukey myself:

(Note: here, the frequency response plots are normalized to DFT bins, as the window's response is independent of the specific sampling frequency used.)

So, there does seem to be some advantage to the Planck window in the eventually steeper rolloff, so I'm in favor of including this in scipy. (I suspect the apparently flat floor may be some discretization noise...)

I haven't looked too closely over your code, but will do soon. In any case, please squish this PR into one commit, and title it according to the developer's guide.

[e-q]

Also, as I discovered when running your code, the plotting goes much faster if you don't plot every point. The source for the plot in my previous comment, which runs pretty quickly, is here: https://gist.github.com/e-q/b5069f22c190a6f72d2d322bac940b60#file-planck-py

[wavexx]

Squashed. As for the visualization, I meant both: it would be nice to show only half of the window for greater detail, and maybe put a fixed window for comparison (boxcar is a bit pointless shape-wise, maybe it would be best to use a well-known window such as hann, or a gaussian). Also in the frequency response, aside for the log axis, it's nice to always compare to the same reference window.

As for the noise in planck at 10e-18, it's probably a float artifact. I'd expect the rolloff to continue indefinitely. The actual point changes depending on the sampling parameters.

[ewmoore]

Is this used anywhere outside of gravitational wave research? Has it actually been used there or is this paper you've linked to just suggesting that they should use it? If it isn't used anywhere else, can you provide some motivation for why we should include it? I'm -something right now, since AFAICT, this isn't a window that anyone is actually using.

[wavexx]

On Wed, Mar 30 2016, Eric Moore notifications@github.com wrote:

Is this used anywhere outside of gravitational wave research? Has it
actually been used there or is this paper you've linked to just
suggesting that they should use it? If it isn't used anywhere else,
can you provide some motivation for why we should include it? I'm
-something right now, since AFAICT, this isn't a window that anyone is
actually using.

I incurred into similar problems as the authors in tremor analysis where
I needed an adjustable tapered window. I experimented with Tukey
initially, but later found also this as an alternative which has a
steeper falloff and better suited my purposes.

I actually didn't check the literature much about other uses (I really
don't care if it has been used only once). I expect one can come up with
dozens of specialty windows, so I see your point. But the definition is
pretty simple and has nice characteristics. I'm surprised to see it
defined only in a 2010 paper honestly.

[e-q]

I can confirm that it is actively used in GW research, but I'm not aware of
its use anywhere else.

On Wed, Mar 30, 2016 at 05:14 Eric Moore notifications@github.com wrote:

Is this used anywhere outside of gravitational wave research? Has it
actually been used there or is this paper you've linked to just suggesting
that they should use it? If it isn't used anywhere else, can you provide
some motivation for why we should include it? I'm -something right now,
since AFAICT, this isn't a window that anyone is actually using.

—
You are receiving this because you commented.

Reply to this email directly or view it on GitHub
#6012 (comment)

[wavexx]

Any consensus? In the same spirit (and if there is interest), I'd also like to contribute the Planck-Bessel variant.

[e-q]

The standard procedure for things like this is to send a message to scipy-dev@scipy.org to try and get more responses and try to build consensus.

[rgommers]

scipy-dev thread: https://mail.scipy.org/pipermail/scipy-dev/2017-January/021639.html

[larsoner]

It sounds like @e-q can confirm it is actively used in gravitational research. @ewmoore does this seem sufficient?

If so this needs a rebase. I'll remove the 1.1.0 milestone as it's probably too tight a timeline.

[ewmoore]

Sure. I'm not opposed to adding windows that are actually being used.

[larsoner]

@wavexx when you get a chance can you rebase this? The window should be added to the signal.windows namespace only nowadays.

[wavexx]

On Wed, Apr 11 2018, Eric Larson wrote: @wavexx when you get a chance can you rebase this? The window should be added to the signal.windows namespace only nowadays.

In a week or so would be too late?

[larsoner]

1.1.0 branches tomorrow, but it's okay to get this into 1.2.0. That would give you time to implement the Planck-Bessel variant, too (assuming that is used in the scientific community, too).

[larsoner]

@wavexx any interest in coming back to this now that scipy/signal/windows/ exists and has its own namespace?

[wavexx]

I'll need to work again in the project where I used the planck taper in a few months. I can revise the implementation then