Spectogram fails for epochs of 0

[skjerns]

I realized that if there is an epoch with 0s, the 2*np.log10(mesh) will create an -inf, which in place prevents the correct scaling of the spectogram, making it more or less invisible.

I propose we replace all NaN and -inf with the lowest number:

visuals.py

            mesh = 20 * np.log10(mesh)
            idx_notfinite = np.isfinite(mesh)==False
            mesh[idx_notfinite] = np.min(mesh[idx_notfinite==False])

What do you think?

[EtienneCmb]

GitMate.io thinks a possibly related issue is https://github.com/EtienneCmb/visbrain/issues/2 (vispy version 0.5?).

[EtienneCmb]

Hi @skjerns ,

You're right, may be something like :

idx_notfinite = ~np.isfinite(mesh)
mesh[idx_notfinite] = 20 * np.log10(mesh[idx_notfinite])

What do you think?

[skjerns]

Mhh, not quite what I meant, don't know if the spectogram-function itself will output any NaNs, but the -inf are coming from log(0), which would not be resolved by your proposal, as they are just 0 before.

But even if we did this, we would still get the scaling wrong, as most values will be something like -80 to -120 while the non-logged values would be 0.

But in this case I propose:

is_zero = mesh==zero
mesh[~is_zero]= 20 * np.log10(mesh[~is_zero])
mesh[is_zero] = mesh.min()

However, I don't know if the mesh itself can already contain NaNs or infs

[EtienneCmb]

You're right, sorry. Can you do the PR?

[EtienneCmb]

Why not something like this :

mesh[mesh <= 0] = 1.
mesh = 20. * np.log10(mesh)

Might be the simplest solution

[skjerns]

Mhh, this way we have the same problem again, we have values of -80 to -120 for regular log(mV) values and 0 for log(1) , changing the whole min/max scaling of the spectogram. Anyway, if you insist on this solution I'm fine with it, it's only a visual thing.

[EtienneCmb]

I'm not really sure to understand your point. mesh is the spectrogram which means that unit is V² / Hz and values should not be < 0?

[skjerns]

Yes, the mesh will always be > 0, but mostly be <1, which will create values of -80 to -120 in the conversion to dB if the signal is in V, but contains data in mV-range. If the signal is 0 at longer periods the log10(x) turn to -inf. Now the plotting will scale the values from mesh.min to mesh.max, which will fail. But I now remember that visbrain has it's values usually scaled in mV? But the problem of scaling persists.

I'll make an example:

import numpy as np
from lspopt import spectrogram_lspopt
# random signal with first half 0 second half some mVs from -500 to +500
data= np.hstack([np.zeros(2**15), np.random.rand(2**15)-0.5]) * 10**3

# standard parameters from visbrain
fs = 128
nperseg = int(30*fs)
overlap = int(0.5 * nperseg)
freq, _, mesh = spectrogram_lspopt(data, fs=fs, nperseg=nperseg, c_parameter=20, noverlap=overlap)

# mesh is now 0 or else most values are around 1000

mesh = 20 * np.log10(mesh)
# mesh is now -inf and else from 26 to 75 dB

Later, the plotting bounds will be -inf to 75, which leads to the plot not being displayed. If we only change the -inf to 0, we scale from 0 to 75, while I argue that we should scale 26 to 75.

This get's more extreme if the original values are smaller than 1 (if we have mV presented as V, ie. values smaller than 1, not sure if this is ever the case in visbrain?), as they are converted to ~-80 to -120 dB

Either way, this is a minor thing, it just makes the plot look nicer if the lower/upper limit of the plot is not a drastic outlier.

I guess alternatively we could also set the _mesh.min(), _mesh.max() below in the code to appropiate values?

[EtienneCmb]

Hi @skjerns ,

Ok ! It could be relatively easy to fix this by adding 1 (something like 20. * np.log10(mesh + 1.)). But as you said, this is just the colobar limits. Anyway, can I let you make a pull request?

[EtienneCmb]

For the _max and _min, we could use something smarter like np.percentile(_mesh, 5) and np.percentile(_mesh, 95)

[EtienneCmb]

One remaining point, but in the future I think it might be better to use the multi-taper of MNE instead of the one of lsopt (@raphaelvallat )

[EtienneCmb]

An other solution, if the only problem is the colorbar limits, it could also be possible to only get the min, max of the finite par of the mesh e.g.

mesh = 20. * np.log10(mesh)
# [...]
_mesh = mesh[sls, :]
# [...]
is_finite = np.isfinite(_mesh)
_min = np.percentile(_mesh[is_finite], 5)
_max = np.percentile(_mesh[is_finite], 95)

[raphaelvallat]

One remaining point, but in the future I think it might be better to use the multi-taper of MNE instead of the one of lsopt (@raphaelvallat )

@EtienneCmb I agree that it would be better to remove the annoying lspopt dependency. However, the MNE multitaper spectrogram (tfr_array_multitaper) is about 1000 times slower than lspopt...e.g. 1 sec versus 1 ms for 30-seconds of data sampled at 100 Hz.

[jasmainak]

Why don't you just copy the code from Ispopt and put it in externals?

[EtienneCmb]

aaaaahhh, when talking about MNE @jasmainak is in the place :D

Why don't you just copy the code from Ispopt and put it in externals?

What do you mean by "externals"? Do you mean "copying the file" and put it inside visbrain?

However, the MNE multitaper spectrogram (tfr_array_multitaper) is about 1000 times slower than lspopt...e.g. 1 sec versus 1 ms for 30-seconds of data sampled at 100 Hz.

@raphaelvallat this is a surprising result as the MNE function is probably highly optimized (multi-core is applied across epochs, which in our case is not really useful). I would say that both uses a different definition of time-frequency resolution which could lead to such differences in computing time...

[jasmainak]

What do you mean by "externals"? Do you mean "copying the file" and put it inside visbrain?

yep :) Usually, we make a folder called externals and the copied packages are inside it. You don't enforce pep8 etc. on these external packages. Just copy it as it is from the original package and update it periodically.

this is a surprising result

yes, perhaps a proper benchmark would help understand what is going on :)

[raphaelvallat]

@jasmainak agree with you that copying-pasting the lspopt core code into Visbrain might be a good option!

Regarding the benchmarks, here are some notable differences between the two methods:

• LSPOPT uses a STFT, whereas MNE uses wavelet convolution; therefore the frequency resolution of LSPOPT depends on the window duration (as in scipy.signal.spectrogram), whereas it is user-specified for MNE. Probably a lower frequency resolution, i.e. 1 wavelet every 2 Hz would speed up the computation. Also, the time resolution will be much higher with MNE (every sample) than LSPOPT (every window).

• LSPOPT approximates the optimal multitapers using Hermite functions and the optimal weights are pre-stored (see more here). Since the computation of the optimal tapers can be quite expensive, I think this is also one of the reason why LSPOPT is faster.

• MNE parallelizes across channels and not across epochs. For Visbrain this is not very useful since the spectrogram is computed on one channel.

[jasmainak]

I see, thanks for the overview! So they aren't strictly comparable. Any reason lsopt uses STFT? aren't wavelets recommended for brain signals?

I do remember that computing the tapers were expensive but there have been some recent optimizations in scipy + mne for this part, so it might be worth checking again.

Also, the time-frequency code in MNE is indeed quite slow now that I think of it. I remember struggling to make it work fast at the beginning of my phd. Any help or optimizations in the code would be welcome I think :)

[skjerns]

@EtienneCmb I've plotted 3 examples.

Columns:

1. current scaling fails on zeros
2. replacing -inf to a percentile value seems to work well for all cases
3. just chanign the scaling doesn't work, as we already do some multiplication with mesh.min() which is -inf before
4. Just taking log(mesh+1) alters the values of the mesh, making the plot in some regular cases appear quite dark

Therefore I suggest to use the second solution, but on the _mesh, to save some computation

is_finite = np.isfinite(_mesh)
_mesh[~is_finite] = np.percentile(_mesh[is_finite], 5)

On my machine this operation increases spectogram creation by 3%, which seems reasonable.

[EtienneCmb]

Really nice @skjerns ! I've two questions :

• The subplot comes from Matplotlib or is it screenshots?
• Did you try replace by NaN?

_mesh[~is_finite] = np.nan

By curiosity, I'm not sure what's going to happened with OpenGL...

[EtienneCmb]

To me, log(mesh + 1.) looks better and you don't have to manually set values. For the contrast issue, we already provide a way to control it from the GUI.

[skjerns]

• it is a screenshot of visbrain
• when I replace with NaN the later call to _mesh.min() will return NaN, destroying the clim, so we would need to replace it with np.nanmean(_mesh), which is in the end the same as np.min(_mesh[is_finite]). But it would save us one call to mesh.min, saving a few 15-30ms.

For me thelog(mesh+1) removes the visibility of the spindles too much (especially in example 2). I do agree that in the last example it might appear better, but bear in mind that the case of epochs with zero is rather rare. My point here is that using the percentile (or alternatively: _mesh[isfinite].min()) does not alter the appearance of the spectogram from the current way they are displayed.

[EtienneCmb]

ok. Can you make the PR with your idea. Thanks ;)