neurodsp-tools/neurodsp

Digital signal processing for neural time series.

neurodsp's Issues

Spectral tutorial and slope function

Missing from the spectral module

Move data for tutorials in the tutorials directory

Currently, the data used in the tutorials is imported from the tests directory, but this will not work if the user first pip-installs neurodsp, as the data is not included in that install.

f_c float/tuple option is really ugly

I think having to juggle between float and tuple is confusing, as oppose to always requiring a tuple and put in 0s for low/high-pass.

Remove waveform shape and PAC functions

Waveform shape analysis should be done using bycycle: https://github.com/voytekresearch/bycycle

PAC should be done using pacpy: https://github.com/voytekresearch/pacpy
or pactools: https://github.com/pactools/pactools

Move plotting functions

There is some plotting scattered around the module - sometimes plotting functions in other module scripts, or plots defined directly inside functions that other things.

This is perhaps not ideal:

• There is no easy way to see, overall what plotting functionality is available in the module
• Plots defined directly in functions break "1 function = 1 task", and also make them harder to apply (can't apply onto existing returned vars, without re-computing)
• Plot functions within specific sub-modules somewhat dilutes the focus of the module (it's not "task" but "task & related plotting").

Recommendation: can we move all plotting things to a plts module (or even, folder with multiple sections), and import as needed? Then everything still works as currently, but plots are more organized, and easier to find / apply / maintain.

Note: If y'all are cool with this - I can do it.

nonhomogeneous poisson process

simulation currently only supports constant-rate poisson population, would be useful to have:

• nonhomogeneous poisson process given a rate vector
• doubly stochastic poisson process (Cox process)

Need helpful error when defined filter length is longer than signal length

If the user tries to define a filter (in neurodsp.filt.filter) that is longer than the signal, then the output of the convolution becomes the length of the filter instead of the length of the original signal, and this will cause an error downstream.

This is why unit tests would be useful.

Add helpful error message when try to simulate too short of a simulated signal

Currently, the user will get an error if they try to simulate a signal that is too short if it relies on a background process that requires filtering (using filtfilt).

If the filter length is longer than the half the length of the signal, then there will be an error.

We should explicitly raise an error telling the user when the simulation duration is too short.

Possibly more ideally, we should just simulate a longer trace and then chop it down.

I don't think this is necessary for v1.0 because there will be no API changes.

Add verbose parameter to burst and timefrequency

We can currently suppress the Transition width printing in filter_signal() by passing in verbose=False.

However, it's currently impossible to suppress this print in some functions that call filter_signal(), namely:

• timefrequency: amp_by_time(), phase_by_time(), freq_by_time()
• burst: detect_bursts_dual_thresh()

These functions should therefore have verbose kwargs to pass through to filter_signal().

This is useful when the user is looping through signals to detect bursts and doesn't want the filter properties to print every single time.

Error when trying to increase signal duration for hilbert transform

There's a phenomena that the hilbert transform function (scipy.signal.hilbert) runs very slowly for some durations of signals. This can be bypassed by specifying the number of points in the FFT. I've noticed that if I make the number of points equal to the next power of 2, that helps.

However, the current code to implement this when calculating an amplitude envelope (_hilbert_ignore_nan in timefrequency.py) forgets to decrease this signal back to the original length.

Remove imports from __init__.py

As discussed with @rdgao and @TomDonoghue , we're going to change the importing so that we will not support imports like from neurodsp import filter_signal but rather force from neurodsp.filt import filter_signal

Unused parameter in pac/compute_pac_comodulogram

In the pac, module the function compute_pac_comodulogram has a parameter 'N_cycles_amp' that is never actually used anywhere in the function.

Without intimately knowing the implementation, I'm not sure if this is just a small miss, perhaps a legacy thing, or indicative of a bigger issue with the function - but we should figure this out, and fix it.

shape.extrema_interpolated_phase() - dealing with boundary artifacts

This function currently deals with the edge artifact of phase estimation by filling in the boundary of the signal with the phase of the extrema/zerocrossing closest to that boundary. This should be replaced with NaNs to make it clear that the phase value is unknown in the boundary of the signal.

The test will also need to be updated.

Cannot Simulate Bursting Behavior with Oscillator Frequency greater than 100 Hz

Code Used:

%load_ext autoreload
#%matplotlib inline
%config InlineBackend.figure_format = 'retina'

import numpy as np
import scipy as sp

from scipy import signal
from scipy import io

from neurodsp import sim
from neurodsp import spectral
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams['figure.figsize'] = 16, 6
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['Tahoma']

n_seconds = 5
fs = 1000
freq = 200
noise_generator = 'synaptic'
noise_args={'n_neurons':1000, 'firing_rate':5, 't_ker':2., 'tau_r':0.002, 'tau_d':0.02}
rdsym= 0.3
ratio_osc_var= 1.
prob_enter_burst=0.7
prob_leave_burst= 0.4
t=np.arange(0, n_seconds, 1/fs)
data= sim.sim_noisy_bursty_oscillator(n_seconds, fs, freq, noise_generator, noise_args, rdsym=rdsym, ratio_osc_var=ratio_osc_var, prob_enter_burst=prob_enter_burst, prob_leave_burst=prob_leave_burst)
plt.plot(t,data)
plt.show()

---

Error Given:

ValueError Traceback (most recent call last)
in
32 ratio_osc_var=ratio_osc_var,
33 prob_enter_burst=prob_enter_burst,
---> 34 prob_leave_burst=prob_leave_burst)
35 plt.plot(t,data)
36 plt.show()

~/anaconda2/envs/BirdSongToolbox/lib/python3.6/site-packages/neurodsp/sim.py in sim_noisy_bursty_oscillator(n_seconds, fs, freq, noise_generator, noise_args, rdsym, ratio_osc_var, prob_enter_burst, prob_leave_burst, cycle_features, return_components, return_cycle_df)
372 prob_leave_burst=prob_leave_burst,
373 cycle_features=cycle_features,
--> 374 return_cycle_df=True)
375
376 # Determine samples of burst so that can compute signal power over only those times

~/anaconda2/envs/BirdSongToolbox/lib/python3.6/site-packages/neurodsp/sim.py in sim_bursty_oscillator(n_seconds, fs, freq, rdsym, prob_enter_burst, prob_leave_burst, cycle_features, return_cycle_df)
260 int(row['period'] / 4))[1:]
261 rise_t = np.cos(rise_pha) * row['amp']
--> 262 sig[-len(rise_t):] = rise_t
263
264 # Add a cycle with rdsym

ValueError: could not broadcast input array from shape (0) into shape (5)

bug - first cycle should never be declared as oscillating

The amp_consistency and period_consistency features used to determine if a cycle is an oscillatory regime for shape analysis are not defined for the first cycle (NaN). However, the algorithm can still mark this cycle as oscillating (is_cycle). Instead, the first cycle should always have is_cycle=False.

Update tutorials for v1.0

Since we're making a bunch of API changes, we should make sure all the tutorials still run like they should.

iir bandstop filtering broken

The following should result in a signal with the 60Hz removed, but instead it returns all NaN.

x2 = neurodsp.filter(x, Fs, 'bandstop', f_lo=58, f_hi=62,
                           iir=True, butterworth_order=3, plot_frequency_response=True,
                           remove_edge_artifacts=False)

Clarification: This is the case if x contains some elements that are NaN. This is not a problem if iir=False.

Add sphinx site

Allow arbitrary background process in `sim_noisy_bursty_oscillator`

Currently sim_noisy_bursty_oscillator() only takes (optionally filtered) colored noise. However, people may want to use the other types of simulated noise.

I think we should try to get this in by v1.0. Because there will be a slight API change (changing kwargs). Can you do this @rdgao ?

Convert tutorials to .py files for sphinx

after PR #111, we just need 2 more:

• filt
• sim

Add real data to burst detection

In addition to demonstrating the burst detection module on simulated data, we should also demonstrate how it works on the example real data ('./data/sample_data_1.npy') in the tutorial

Test Coverage

Missing test coverage:

• spectral.morlet_transform
• spectral.morlet_convolve
• burst.compute_burst_stats
• plts.*

if n_seconds is non-integer, simulation breaks

Looks like its slicing over n_seconds, so if n_seconds is non-integer (e.g., 2 vs 2.) it breaks.

---

TypeError Traceback (most recent call last)
in ()
18 ratio_osc_var=ratio_osc_var,
19 prob_enter_burst=prob_enter_burst,
---> 20 prob_leave_burst=prob_leave_burst)
21
22 plt.plot(t, data)

~/anaconda/envs/py3k/lib/python3.5/site-packages/neurodsp/sim.py in sim_noisy_bursty_oscillator(n_seconds, fs, freq, noise_generator, noise_args, rdsym, ratio_osc_var, prob_enter_burst, prob_leave_burst, cycle_features, return_components, return_cycle_df)
372 prob_leave_burst=prob_leave_burst,
373 cycle_features=cycle_features,
--> 374 return_cycle_df=True)
375
376 # Determine samples of burst so that can compute signal power over only those times

~/anaconda/envs/py3k/lib/python3.5/site-packages/neurodsp/sim.py in sim_bursty_oscillator(n_seconds, fs, freq, rdsym, prob_enter_burst, prob_leave_burst, cycle_features, return_cycle_df)
280 sig = np.append(sig, cycle_t)
281 last_cycle_oscillating = True
--> 282 sig = sig[:n_samples]
283
284 if return_cycle_df:

TypeError: slice indices must be integers or None or have an index method

sim suggestions

I'm starting to use neurodsp.sim, and have some possible suggestions:

ToDos:

• Add in general colored noise generation (implementation available here: https://github.com/felixpatzelt/colorednoise) [Note: Richard added his own implementation].
• The 'Input Suggestion' docs for sim_oscillator or wrong / misplaced, right? Am I missing something? [Fixed: with more general doc clean ups.]
• I'll probably update some variable names, to update to our API conventions (ex - brownNf -> brown_nf or similar).

Open Questions:

• Is there any reason not to generalize all the functions that currently add 1/f^2 noise, to be able to add 1/f^n noise, with set-able n?
• Mixed conventions: some functions take in a signal length, others a number of samples. We should consolidate on one approach - any suggestions on which is better?

Anyways - I'll use this issue for continuing points of discussion, and open a PR soon-ish with some updates - so let me know of any thoughts about things here.

Wrong polarity of peak and trough voltage in trough-centered cycles

The peak and trough voltages need to be inversed in the cycle by cycle dataframe. This is only the case for trough-centered cycles. This is an issue in the cycle_features.py file in the shape module.

clean bursting module

Going through PR #80, I noticed a few changes I think should be made to prepare burst.py for v1.0

• detect_bursts()
  • Should be renamed to something like detect_bursts_dual_amplitude_threshold to be consistent with detect_bursts_bosc
  • Variable name improvements
    • X → signal
    • Min_osc_periods → min_cycles
    • Devation_type → average_method
  • Threshold parameters should be required
  • Remove algorithm input (only the median method)
  • Remove filter_fn optional input (unnecessary, and only complicates things)
• detect_bursts_bosc() I think this should either be labeled as experimental, private, or taken out since it really hasn't been tested or used much at all. Maybe just taken out cause it's a published method by another group (though there is no public code I can find for it)
• get_stats() - rename to something more descriptive like compute_bursting_statistics and could change a few things

New module: simulating oscillatory neural signals

As discussed with @rdgao.

We should have a function (or set of functions) that allow us to construct signals that contain oscillations with various properties (e.g. stationary or bursty) along with noise (probably brown noise), and control the relative power between these two.

Some current code snippets below can be used as a basis.

Simulate a stationary oscillator with arbitrary rise-decay symmetry:

import numpy as np
from scipy import signal

rdsym = .4
period_ms = 100
Fs = 1000
N_cycles = 50

rise_ms = np.round(period_ms * rdsym)
decay_ms = period_ms - rise_ms

pha_one_cycle = np.hstack([np.linspace(0, np.pi, decay_ms+1), np.linspace(-np.pi, 0, rise_ms+1)[1:-1]])
phase_t = np.tile(pha_one_cycle, N_cycles)
x_sawtooth = np.cos(phase_t)
t = np.arange(0, len(x_sawtooth)/Fs, 1/Fs)

High-pass-filtered brown noise:

def simfiltonef(T, f_range, Fs, N):
    """ Simulate a band-pass filtered signal with 1/f^2 
    Input suggestions: f_range=(2,None), Fs=1000, N=1000
    
    Parameters
    ----------
    T : float
        length of time of simulated oscillation
    Fs : float
        oscillation sampling rate
    f_range : 2-element array (lo,hi)
        frequency range of simulated data
        if None: do not filter
    N : int
        order of filter
    """

    if f_range is None:
        # Do not filter
        # Generate 1/f^2 noise
        brownN = simbrown(int(T*Fs))
        return brownN
    elif f_range[1] is None:
        # Make filter order odd if necessary
        nyq = Fs / 2.
        if N % 2 == 0:
            print('NOTE: Increased high-pass filter order by 1 in order to be odd')
            N += 1
            
        # Generate 1/f^2 noise
        brownN = simbrown(int(T*Fs+N*2))

        # High pass filter
        taps = signal.firwin(N, f_range[0] / nyq, pass_zero=False)
        brownNf = signal.filtfilt(taps, [1], brownN)
        return brownNf[N:-N]

    else:
        # Bandpass filter
        # Generate 1/f^2 noise
        brownN = simbrown(int(T*Fs+N*2))
        # Filter
        nyq = Fs / 2.
        taps = signal.firwin(N, np.array(f_range) / nyq, pass_zero=False)
        brownNf = signal.filtfilt(taps, [1], brownN)
        return brownNf[N:-N]
    

def simbrown(N):
    """Simulate a brown noise signal (power law distribution 1/f^2)
    with N samples"""
    wn = np.random.randn(N)
    return np.cumsum(wn)

T = len(phase_t) / Fs
f_range = (2, None)
N = 1000
x_brown = simfiltonef(T, f_range, Fs, N)

Combine these signals and control relative power:

sawtooth_power = np.mean(x_sawtooth**2)
brown_power = np.mean(x_brown**2)
power_ratio = .2 # Higher power ratio is more noise
x_brown = np.sqrt(x_brown**2 * power_ratio * sawtooth_power / brown_power) * np.sign(x_brown)
x = x_sawtooth + x_brown

This code creates a stationary oscillator (top), noise (middle), and a composite signal (bottom).

different simulations of same params for sim_noisy_bursty_oscillator yield totally different (arbitrary) y-values

instead of oscillating around 0 (as usually happens in raw LFP data, for example), different simulations of the same params in sim_noisy_bursty_oscillator result in oscillations around y-values of anywhere between 300 to 350, but it seems totally arbitrary which of these the signal will oscillate around on any one simulation.

spectral CV

New simulations

Just throwing out a couple ideas of possible things to add to sims that I would be interested in talking about (at some point), and possibly working on implementing to be able to play with:

• simulate transients (ERP-like)
• simulate time-varying 1/f backgrounds

Note: edited above to more clearly reflect ideas.

Split general filter function into specialized functions

We talked about how the current filter function is a bit of a mess to deal with because of the different requirements for bandpass, highpass, lowpass, bandstop.

After discussing, @TomDonoghue and @rdgao and I generally agreed that it would be cleaner if we split them up into different sub-functions that could alternatively be called in addition to the general function

phase-amplitude coupling

Add burst characterization function

Add a function to the burst module that outputs statistics of bursting in a signal (e.g. number of burst, mean burst length, etc.)

Takes as an input the binary array output by the associated burst detection function.

Add defaults for noise process in simulating noisy oscillators

Right now, the user has to specify the noise process and ALL of its parameters when running sim_noisy_oscillator() or sim_noisy_bursty_oscillator().

I think it would be good to have some default noise process and parameters here. Defaults would be useful when just wanting to test something out quick on some simulated data.

I don't have much of a preference of what the default is. How about pink noise?

I think it would be good for the user to be able to just choose the noise method and whatever (if any) parameters they want to specify, and not have to specify all of them.

So, users should be able to run these commands without error:
sim_noisy_bursty_oscillator(n_seconds, fs, freq)
sim_noisy_bursty_oscillator(n_seconds, fs, freq, 'filtered_powerlaw')
sim_noisy_bursty_oscillator(n_seconds, fs, freq, 'filtered_powerlaw', {'exponent': -2})

For the last case of incomplete argument specification, we'll need to call the noise simulation functions in _return_noise_sim() differently, like, for example:

n_seconds = 100
fs = 1000
noise_args = {'exponent': -2}
sim_filtered_noise(n_seconds, fs, **noise_args)

Final burst cycle can end later than final sample in sim_noisy_bursty_oscillator() (sim.py)

Either a bug or should be noted in documentation. Likely also apparent in other bursty simulator functions.

power spectrum

Variable names

We don't have super consistent variable naming. With a breaking 1.0 release, should we go through and systematize it all to snake_case (so, for example, things like 'N_samples' -> 'n_samples').

Also, some specific cases:

• filter overwrites the python filter function, which can be really problematic.
  • Especially because NDSP by default imports all functions directly in local namespace, this squashes filter always, meaning you can never use Python filter, and so has side effects / breaks things
• psd, as a function, strikes me as not a great name. Elsewhere, psd is used as a variable name, and it sounds like data more than a function
  • Heuristic: function names should be 'verb-y', because they do things.

We also have a bunch of single character variables, but this (in some cases), might be fine for the math-y parts, and also some names, like 'Fs', that are non-compliant, but (I think), are inherited from style in scipy.signal, or similar. Do we want to update those names too?

bandstop filter error

If the user defines a bandstop filter with incompatible time and frequency resolution, then the current filter function throws an error when the transition bands cannot be estimated. This is now replaced with a warning.

pep8

Assure current code follows the pep8 style guide.

filter error unhelpful for transition bandwidth

currently, when a filter length is not long enough for computing transition bandwidth (I assume), it throws an unhelpful error, after the bandwidth warning. This occurs in bandstop filtering but presumably generalizes to the other options too for FIR.

From @tunmiseo

=UserWarning: Error computing transition bandwidth of the filter. Defined filter length may be too short.
raise warnings.warn(‘Error computing transition bandwidth of the filter. Defined filter length may be too short.’)
/code/lib/neurodsp/neurodsp/filt.py in filter(x, Fs, pass_type, fc, N_cycles, N_seconds, iir, butterworth_order, plot_frequency_response, return_kernel, verbose, compute_transition_band, remove_edge_artifacts)
238 warnings.warn(‘Transition bandwidth is ’ + str(np.round(transition_bw, 1)) + ' Hz. This is greater than the desired pass/stop bandwidth of ' + str(np.round(pass_bw, 1)) + ' Hz’)
239 except StopIteration:
--> 240 raise warnings.warn(‘Error computing transition bandwidth of the filter. Defined filter length may be too short.’)
241
242 # Remove edge artifacts

TypeError: exceptions must derive from BaseException

Add burst detection method: amplitude threshold computed from background fit

As opposed to choosing the amplitude threshold for burst detection off the variance of the amplitude envelope, some promising methods fit a line to the power spectrum to determine the baseline amplitude level, if no oscillation is present. We should implement this in our burst module.

2 examples are:

• Better OSCillation Detector (BOSC)
  • https://www.sciencedirect.com/science/article/pii/S1053811910011614?via%3Dihub
• Multiple Oscillations Detection ALgorithm (MODAL)
  • https://www.biorxiv.org/content/early/2017/10/12/202374
  • http://www.andrewwatrous.com/uploads/7/8/5/9/78599468/modal.m

SOP documentation

Add a markdown file that describes expectations for code contributions to this package. This includes:

• Documentation Numpy docs
• Coding style. pep8
• Testing. At least with example data to test consistency in results across versions. Should pass TravisCI.
• Tutorial notebook requirement
• Dependencies. Attempt to avoid, especially if outside of Anaconda
• Any standard APIs, e.g. sampling rate denoted as 'Fs', time series denoted as 'x', etc.

lowpass and highpass parameters flipped around?

@srcole
in filt.py:

isn't a lowpass usually defined such that f_lo is 0 and f_hi is the cut off (i.e. [0 f_hi]), and vice versa for highpass ([f_lo, fs/2])?

in noise sim, t_ker should be a function of tau_r and tau_d

We should algorithmically define t_ker as a function of tau_r and tau_d, rather than requiring users to define it manually.

Option whether or not to remove edge artifacts in timefrequency functions

It would be nice if by default the functions in time_frequency.py like amp_by_time by default removed edge artifacts (fill with NaN) and had an option that kept the edge artifacts. Currently, these functions keep edge artifacts.

Simulating filtered noise breaks when n_seconds is too small

sim_filtered_noise() breaks for some simulations that users may want to run and doesn't give a particularly helpful error. We can change the code to get past this error.

Basically, if the signal is not 3x longer than the filter, then an error occurs, e.g.

n_seconds = 10
fs = 1000
exponent = -2
sim_filtered_noise(n_seconds, fs, exponent, f_range=(2, None), filter_order=5001)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-20-584a47fbfc12> in <module>()
      2 fs = 1000
      3 exponent = -2
----> 4 sim_filtered_noise(n_seconds, fs, exponent, f_range=(2, None), filter_order=5001)

/gh/bv/neurodsp/neurodsp/sim.py in sim_filtered_noise(n_seconds, fs, exponent, f_range, filter_order)
    794         # High pass filter
    795         taps = signal.firwin(filter_order, f_range[0] / nyq, pass_zero=False)
--> 796         noise = signal.filtfilt(taps, [1], noise)
    797 
    798     # Band pass filtered

~/anaconda/lib/python3.6/site-packages/scipy/signal/signaltools.py in filtfilt(b, a, x, axis, padtype, padlen, method, irlen)
   3097     # method == "pad"
   3098     edge, ext = _validate_pad(padtype, padlen, x, axis,
-> 3099                               ntaps=max(len(a), len(b)))
   3100 
   3101     # Get the steady state of the filter's step response.

~/anaconda/lib/python3.6/site-packages/scipy/signal/signaltools.py in _validate_pad(padtype, padlen, x, axis, ntaps)
   3147     if x.shape[axis] <= edge:
   3148         raise ValueError("The length of the input vector x must be at least "
-> 3149                          "padlen, which is %d." % edge)
   3150 
   3151     if padtype is not None and edge > 0:

ValueError: The length of the input vector x must be at least padlen, which is 15003.

The signal is 10,000 samples long (n_seconds * fs), and scipy.signal.filtfilt() seems to require at least 3 times the length of the filter (5,001 samples, in this case).

One issue. Is that we should not be using filtfilt. We should be using our own highpass/bandpass filter functionality that uses one-pass filtering, to be consistent.

We also need to deal with edge artifacts of the filter. To simulate 10 seconds of a signal that has a highpass filter of 3 seconds, we should simulate 16 seconds, and then chop off the edge artifacts.

amp_by_time fix to work with signals with NaN

The default filter method replaces regions prone to edge artifact with NaN. However, if we want to estimate the amplitude with amp_by_time, then the amplitude array is all NaN. It should be fixed to properly estimate amplitude.

This is likely also the case for phase_by_time

Clarification:
The amp_by_time and phase_by_time functions do the bandpass filtering, but perhaps the user will have done some low-pass filtering that will have caused NaNs in the signal and want to use that as the input time series.

pypi description formatting - setup.py

As seen at https://pypi.org/project/neurodsp/

The second paragraph is not formatted in bullets as intended.

The final paragraph is supposed to be metadata for pypi but is just in the description.

See bycycle or fooof for functional setup.py files. They look the same to me, so I'm not sure what the difference is that makes the neurodsp pypi page not render correctly.

Burst detection

Add a module for burst detection. The main function should take as input a signal and frequency range of interest, and output the samples at which the oscillation is present

PYPI stuff

I think we can drop the egg-info file on master branch. AFAIK it's not needed - it's actually out of date (it's the info for v0.1 release), so that actually makes it sort of confusing to have hanging around, with the current tagged version (also on PYPI) being 0.2.

Also: we should probably add a description to show up on PYPI:
https://pypi.org/project/neurodsp/