neurodsp-tools / neurodsp Goto Github PK
View Code? Open in Web Editor NEWDigital signal processing for neural time series.
Home Page: https://neurodsp-tools.github.io/
License: Apache License 2.0
Digital signal processing for neural time series.
Home Page: https://neurodsp-tools.github.io/
License: Apache License 2.0
Missing from the spectral module
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.
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.
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
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:
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.
simulation currently only supports constant-rate poisson population, would be useful to have:
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.
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.
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:
amp_by_time()
, phase_by_time()
, freq_by_time()
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.
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.
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
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.
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.
%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()
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)
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.
Since we're making a bunch of API changes, we should make sure all the tutorials still run like they should.
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
.
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 ?
after PR #111, we just need 2 more:
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
Missing test coverage:
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
I'm starting to use neurodsp.sim, and have some possible suggestions:
ToDos:
Open Questions:
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.
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.
Going through PR #80, I noticed a few changes I think should be made to prepare burst.py
for v1.0
detect_bursts()
detect_bursts_dual_amplitude_threshold
to be consistent with detect_bursts_bosc
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 thingsAs 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).
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.
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:
Note: edited above to more clearly reflect ideas.
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
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.
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)
Either a bug or should be noted in documentation. Likely also apparent in other bursty simulator functions.
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.
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
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?
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.
Assure current code follows the pep8 style guide.
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
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:
Add a markdown file that describes expectations for code contributions to this package. This includes:
@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])?
We should algorithmically define t_ker as a function of tau_r and tau_d, rather than requiring users to define it manually.
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.
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.
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.
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.
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
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/
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.