elephant icon indicating copy to clipboard operation
elephant copied to clipboard

[Bug] peak_detection without analogsignals.times

Open Samcrl0 opened this issue 1 year ago • 4 comments

Describe the bug Hello. After acquired data with an Intan, I would like to use Elephant for the first treatment of these.
In this way, I would like to fill the spiketrain object of neo using the spike_extraction method of elephant, which is empty after the acquisition. Then I could realize the spike sorting. However, Elephant require a list "analogsinal.times" which does not exist in the obect analogsignals of neo.

How the method is supposed to work?

To Reproduce 1. 2. Traceback (most recent call last):

File ~\anaconda3\lib\site-packages\spyder_kernels\py3compat.py:356 in compat_exec exec(code, globals, locals)

File c:\users\clement\desktop\analysis_data\testephy_3.py:63 event = elephant.spike_train_generation.peak_detection(signal_highpassed[0],threshold=-40 * uV, sign="below")

File ~\anaconda3\lib\site-packages\elephant\spike_train_generation.py:321 in peak_detection result_st = neo.SpikeTrain(events_base, units=signal.times.units,

AttributeError: 'Quantity' object has no attribute 'times'

Expected behavior I expect that the spike_extraction method put all the spike events in the spiketrains obect.

Environment

  • Windows 11
  • How you installed elephant : pip
  • Python version: 3.10
  • neo python package version: 0.12.0
  • elephant python package version: 1.0.0
  • (Any additional python package you want to include here)

Samcrl0 avatar Jan 10 '24 09:01 Samcrl0

Hi Samcrl0, thank you for reaching out to us !

The Neo AnalogSignal class implements a times property, as one can see by inspecting the source code at https://github.com/NeuralEnsemble/python-neo/blob/4ef280130be875263878d5e55dbf7ee40f950c9b/neo/core/analogsignal.py#L387. So the case where a Neo AnalogSignal has no times attribute should be impossible.

It looks like the issue is related to signal_highpassed, which appears to be a list of quantity arrays based on the traceback.

To address this, you can verify if the type of signal_highpassed[0] is a Neo AnalogSignal with:

print(type(signal_highpassed[0]))

If not, you might try to create a Neo AnalogSignal from a Quantity array by following the instructions outlined in the Neo documentation here: https://neo.readthedocs.io/en/latest/api_reference.html#neo.core.AnalogSignal

The peak_detection function in Elephant is designed to operate as you explained: It accepts a Neo AnalogSignal as input and outputs a Neo SpikeTrain.

Consider the following example:

from neo.core import AnalogSignal
import quantities as pq
from elephant.spike_train_generation import peak_detection

# Create a Neo AnalogSignal
analogsignal = AnalogSignal([1, 2, 3], units='V', sampling_rate=1*pq.Hz)

# Use peak_detection function from Elephant
result = peak_detection(analogsignal)

# Check if result is a SpikeTrain
print(type(result))

For me this returns:

<class 'neo.core.spiketrain.SpikeTrain'>

Moritz-Alexander-Kern avatar Jan 10 '24 15:01 Moritz-Alexander-Kern

Hi Moritz-Alexander-Kern, Thank you for your answers, it helped me to resolve the issue. However, I get new ones:

Applying this:

signals = data[0].segments[0].analogsignals[0]
signal_highpassed=elephant.signal_processing.butter(signals, highpass_frequency=250 *Hz, sampling_frequency=30000)
event = elephant.spike_train_generation.peak_detection(signal_highpassed,threshold=-40 * uV, sign="below")

where event is :

<class 'neo.core.spiketrain.SpikeTrain'>

I get the peak times in the variable event. But it is a 1-dimensional array, which show all the peak times of all my 16 electrode, I think.

print(event)
[2.26666667e-03 2.91333333e-02 6.50066667e-01 ... 2.99976133e+02
 2.99976200e+02 3.00002100e+02] s

What I expected to obtain was a 2-dimensional object with each spike times for each electrode. Is it possible to get it?

In another side, I tryed the function peak_extraction :

event = elephant.spike_train_generation.spike_extraction(signal_highpassed,threshold= -40 *uV, sign = "below")

but here is the error message:


  Cell In[10], line 1
    event = elephant.spike_train_generation.spike_extraction(signal_highpassed,threshold= -40 *uV, sign = "below")

  File ~\anaconda3\lib\site-packages\elephant\spike_train_generation.py:179 in spike_extraction
    waveforms = waveforms[:, np.newaxis, :]

  File ~\anaconda3\lib\site-packages\quantities\quantity.py:394 in __getitem__
    ret = super().__getitem__(key)

IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed

I do not know what to put for the time_stamp argument in order to have the waveform. I tryed a lot of thing, but I never succeeded.

Do you have some advices, please?

Samcrl0 avatar Jan 12 '24 09:01 Samcrl0

Hi @Samcrl0 ,

What I expected to obtain was a 2-dimensional object with each spike times for each electrode. Is it possible to get it?

Yes, for example by using a for loop and apply peak_detection on the signal from every electrode individually:

from neo.core import AnalogSignal
import quantities as pq
from elephant.spike_train_generation import peak_detection

# Create a Neo AnalogSignal with one electrode
analogsignal_1_electrode = AnalogSignal([0, 1, 0, 2, 0, 1, 0], units=pq.V, sampling_rate=1*pq.Hz)

# Use peak_detection function from Elephant
result_1_electrode = peak_detection(analogsignal_1_electrode)

# Check result
print(f"Signal 1 electrode: {analogsignal_1_electrode}")
print(f"Result with 1 electrode: {result_1_electrode}")


# Create a Neo AnalogSignal with 2 electrodes
analogsignal_2_electrodes = AnalogSignal([[0, 1, 0, 1, 0, 1, 0],
                                          [1, 0, 1, 0, 1, 0, 1]], units=pq.V, sampling_rate=1*pq.Hz)

# Use peak_detection function from Elephant
result_2_electrodes = peak_detection(analogsignal_2_electrodes)

# Check result
print(f"Signal with 2 electrodes: {analogsignal_2_electrodes}")
print(f"Result on 2 electrodes: {result_2_electrodes}")

print("Use peak_detection on single electrodes with for-loop")
results = []
for idx, signal in enumerate(analogsignal_2_electrodes):
    print(f"Signal {idx}: {signal}")
    # Create AnalogSignal
    signal = AnalogSignal(signal, sampling_rate=analogsignal_2_electrodes.sampling_rate)
    # Use peak_detection
    result = peak_detection(signal)
    print(f"Result {idx}: {result}")
    # Append result to results list
    results.append(result)

print(f"Final result: {results}")

Which yields:

Signal 1 electrode: [0, 1, 0, 2, 0, 1, 0] V
Result with 1 electrode: [1. 3. 5.] s
Signal with 2 electrodes: [[0 1 0 1 0 1 0]
 [1 0 1 0 1 0 1]] V
Result on 2 electrodes: [1.] s
Use peak_detection on single electrodes with for-loop
Signal 0: [0 1 0 1 0 1 0] V
Result 0: [1. 3. 5.] s
Signal 1: [1 0 1 0 1 0 1] V
Result 1: [0. 2. 4. 6.] s
Final result: [<SpikeTrain(array([1., 3., 5.]) * s, [0.0 s, 7.0 s])>, <SpikeTrain(array([0., 2., 4., 6.]) * s, [0.0 s, 7.0 s])>]

The final result is a list with 2 Neo SpikeTrain objects containing the events for each electrode respectively.

Moritz-Alexander-Kern avatar Jan 12 '24 14:01 Moritz-Alexander-Kern

The spike_extraction function also requires an AnalogSignal as input.

But there seems to be a mistake in the documentation, the time_stamps argument should be a Neo SpikeTrain, see:

from neo.core import AnalogSignal
import quantities as pq
from elephant.spike_train_generation import spike_extraction, peak_detection

# Create a Neo AnalogSignal with one electrode
analogsignal_1_electrode = AnalogSignal([0, 1, 0, 2, 0, 1, 0], units=pq.V, sampling_rate=1*pq.Hz)

# Use spike_extraction function from Elephant
result_1_electrode = spike_extraction(analogsignal_1_electrode, interval=(-1 * pq.s, 2 * pq.s))

# Check result
print(f"Signal electrode: {analogsignal_1_electrode}")
print(f"Events: {result_1_electrode}")
print(f"Waveforms: {result_1_electrode.waveforms}")

# Use the time_stamp argument:

time_stamps = peak_detection(analogsignal_1_electrode)
print(f"time_stamps: {time_stamps}, {type(time_stamps)}")
result = spike_extraction(analogsignal_1_electrode, time_stamps=time_stamps, interval=(-1 * pq.s, 2 * pq.s))
print(f"Events: {result}")
print(f"Waveforms: {result.waveforms}")

Which yields:

Signal electrode: [[0]
 [1]
 [0]
 [2]
 [0]
 [1]
 [0]] V
Events: [1. 3. 5.] s
Waveforms: [[[[0.]
   [1.]
   [0.]]]


 [[[0.]
   [2.]
   [0.]]]


 [[[0.]
   [1.]
   [0.]]]] V
time_stamps: [1. 3. 5.] s, <class 'neo.core.spiketrain.SpikeTrain'>
Events: [1. 3. 5.] s
Waveforms: [[[[0.]
   [1.]
   [0.]]]


 [[[0.]
   [2.]
   [0.]]]


 [[[0.]
   [1.]
   [0.]]]] V

The error seems similar to the first post, here signal_highpassed is possibly also a list of quantity arrays based on the trace back ?

Moritz-Alexander-Kern avatar Jan 12 '24 14:01 Moritz-Alexander-Kern