mbergmann/gpgrewe2024
1
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

GP_code nach code verschoben

Browse Source
This commit is contained in:
mbergmann 2024-10-22 07:16:25 +00:00
parent b0adb74ae3
commit e8d4ca90fe
1 changed files with 192 additions and 192 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

384
GP_Code.py → code/GP_Code.py
View File

@ -1,192 +1,192 @@
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
""" """
Created on Thu Oct 17 09:23:10 2024 Created on Thu Oct 17 09:23:10 2024
@author: diana @author: diana
""" """
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
import glob import glob
import os import os
import rlxnix as rlx import rlxnix as rlx
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import scipy.signal as sig import scipy.signal as sig
from scipy.integrate import quad from scipy.integrate import quad
### FUNCTIONS ### ### FUNCTIONS ###
def binary_spikes(spike_times, duration, dt): def binary_spikes(spike_times, duration, dt):
""" Converts the spike times to a binary representation. """ Converts the spike times to a binary representation.
Zeros when there is no spike, One when there is. Zeros when there is no spike, One when there is.
Parameters Parameters
---------- ----------
spike_times : np.array spike_times : np.array
The spike times. The spike times.
duration : float duration : float
The trial duration. The trial duration.
dt : float dt : float
the temporal resolution. the temporal resolution.
Returns Returns
------- -------
binary : np.array binary : np.array
The binary representation of the spike times. The binary representation of the spike times.
""" """
binary = np.zeros(int(np.round(duration / dt))) #Vektor, der genauso lang ist wie die stim time binary = np.zeros(int(np.round(duration / dt))) #Vektor, der genauso lang ist wie die stim time
spike_indices = np.asarray(np.round(spike_times / dt), dtype=int) spike_indices = np.asarray(np.round(spike_times / dt), dtype=int)
binary[spike_indices] = 1 binary[spike_indices] = 1
return binary return binary
def firing_rate(binary_spikes, box_width, dt=0.000025): def firing_rate(binary_spikes, box_width, dt=0.000025):
"""Calculate the firing rate from binary spike data. """Calculate the firing rate from binary spike data.
This function computes the firing rate using a boxcar (moving average) This function computes the firing rate using a boxcar (moving average)
filter of a specified width. filter of a specified width.
Parameters Parameters
---------- ----------
binary_spikes : np.array binary_spikes : np.array
A binary array representing spike occurrences. A binary array representing spike occurrences.
box_width : float box_width : float
The width of the box filter in seconds. The width of the box filter in seconds.
dt : float, optional dt : float, optional
The temporal resolution (time step) in seconds. Default is 0.000025 seconds. The temporal resolution (time step) in seconds. Default is 0.000025 seconds.
Returns Returns
------- -------
rate : np.array rate : np.array
An array representing the firing rate at each time step. An array representing the firing rate at each time step.
""" """
box = np.ones(int(box_width // dt)) box = np.ones(int(box_width // dt))
box /= np.sum(box) * dt #Normalization of box kernel to an integral of 1 box /= np.sum(box) * dt #Normalization of box kernel to an integral of 1
rate = np.convolve(binary_spikes, box, mode="same") rate = np.convolve(binary_spikes, box, mode="same")
return rate return rate
def powerspectrum(rate, dt): def powerspectrum(rate, dt):
"""Compute the power spectrum of a given firing rate. """Compute the power spectrum of a given firing rate.
This function calculates the power spectrum using the Welch method. This function calculates the power spectrum using the Welch method.
Parameters Parameters
---------- ----------
rate : np.array rate : np.array
An array of firing rates. An array of firing rates.
dt : float dt : float
The temporal resolution (time step) in seconds. The temporal resolution (time step) in seconds.
Returns Returns
------- -------
frequency : np.array frequency : np.array
An array of frequencies corresponding to the power values. An array of frequencies corresponding to the power values.
power : np.array power : np.array
An array of power spectral density values. An array of power spectral density values.
""" """
frequency, power = sig.welch(rate, fs=1/dt, nperseg=2**15, noverlap=2**14) frequency, power = sig.welch(rate, fs=1/dt, nperseg=2**15, noverlap=2**14)
return frequency, power return frequency, power
def prepare_harmonics(frequencies, categories, num_harmonics, colors): def prepare_harmonics(frequencies, categories, num_harmonics, colors):
points_categories = {} points_categories = {}
for idx, (freq, category) in enumerate(zip(frequencies, categories)): for idx, (freq, category) in enumerate(zip(frequencies, categories)):
points_categories[category] = [freq * (i + 1) for i in range(num_harmonics[idx])] points_categories[category] = [freq * (i + 1) for i in range(num_harmonics[idx])]
points = [p for harmonics in points_categories.values() for p in harmonics] points = [p for harmonics in points_categories.values() for p in harmonics]
color_mapping = {category: colors[idx] for idx, category in enumerate(categories)} color_mapping = {category: colors[idx] for idx, category in enumerate(categories)}
return points, color_mapping, points_categories return points, color_mapping, points_categories
def plot_power_spectrum_with_integrals(frequency, power, points, delta, color_mapping, points_categories): def plot_power_spectrum_with_integrals(frequency, power, points, delta, color_mapping, points_categories):
"""Create a figure of the power spectrum with integrals highlighted around specified points. """Create a figure of the power spectrum with integrals highlighted around specified points.
This function creates a plot of the power spectrum and shades areas around This function creates a plot of the power spectrum and shades areas around
specified harmonic points to indicate the calculated integrals. specified harmonic points to indicate the calculated integrals.
Parameters Parameters
---------- ----------
frequency : np.array frequency : np.array
An array of frequencies corresponding to the power values. An array of frequencies corresponding to the power values.
power : np.array power : np.array
An array of power spectral density values. An array of power spectral density values.
points : list points : list
A list of harmonic frequencies to highlight. A list of harmonic frequencies to highlight.
delta : float delta : float
Half-width of the range for integration around each point. Half-width of the range for integration around each point.
color_mapping : dict color_mapping : dict
A mapping of point categories to colors. A mapping of point categories to colors.
points_categories : dict points_categories : dict
A mapping of categories to lists of points. A mapping of categories to lists of points.
Returns Returns
------- -------
fig : matplotlib.figure.Figure fig : matplotlib.figure.Figure
The created figure object. The created figure object.
""" """
fig, ax = plt.subplots() fig, ax = plt.subplots()
ax.plot(frequency, power) ax.plot(frequency, power)
integrals = [] integrals = []
for point in points: for point in points:
indices = (frequency >= point - delta) & (frequency <= point + delta) indices = (frequency >= point - delta) & (frequency <= point + delta)
integral = np.trapz(power[indices], frequency[indices]) integral = np.trapz(power[indices], frequency[indices])
integrals.append(integral) integrals.append(integral)
# Get color based on point category # Get color based on point category
color = next((c for cat, c in color_mapping.items() if point in points_categories[cat]), 'gray') color = next((c for cat, c in color_mapping.items() if point in points_categories[cat]), 'gray')
ax.axvspan(point - delta, point + delta, color=color, alpha=0.3, label=f'{point:.2f} Hz') ax.axvspan(point - delta, point + delta, color=color, alpha=0.3, label=f'{point:.2f} Hz')
print(f"Integral around {point:.2f} Hz: {integral:.5e}") print(f"Integral around {point:.2f} Hz: {integral:.5e}")
ax.set_xlim([0, 1200]) ax.set_xlim([0, 1200])
ax.set_xlabel('Frequency (Hz)') ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('Power') ax.set_ylabel('Power')
ax.set_title('Power Spectrum with marked Integrals') ax.set_title('Power Spectrum with marked Integrals')
ax.legend() ax.legend()
return fig return fig
### Data retrieval ### ### Data retrieval ###
datafolder = "../data" # Geht in der Hierarchie einen Ordern nach oben (..) und dann in den Ordner 'data' datafolder = "../data" # Geht in der Hierarchie einen Ordern nach oben (..) und dann in den Ordner 'data'
example_file = os.path.join("..", "data", "2024-10-16-ad-invivo-1.nix") example_file = os.path.join("..", "data", "2024-10-16-ad-invivo-1.nix")
dataset = rlx.Dataset(example_file) dataset = rlx.Dataset(example_file)
sams = dataset.repro_runs("SAM") sams = dataset.repro_runs("SAM")
sam = sams[2] sam = sams[2]
## Daten für Funktionen ## Daten für Funktionen
df = sam.metadata["RePro-Info"]["settings"]["deltaf"][0][0] df = sam.metadata["RePro-Info"]["settings"]["deltaf"][0][0]
stim = sam.stimuli[1] stim = sam.stimuli[1]
potential, time = stim.trace_data("V-1") potential, time = stim.trace_data("V-1")
spikes, _ = stim.trace_data("Spikes-1") spikes, _ = stim.trace_data("Spikes-1")
duration = stim.duration duration = stim.duration
dt = stim.trace_info("V-1").sampling_interval dt = stim.trace_info("V-1").sampling_interval
### Anwendung Functionen ### ### Anwendung Functionen ###
b = binary_spikes(spikes, duration, dt) b = binary_spikes(spikes, duration, dt)
rate = firing_rate(b, box_width=0.05, dt=dt) rate = firing_rate(b, box_width=0.05, dt=dt)
frequency, power = powerspectrum(b, dt) frequency, power = powerspectrum(b, dt)
## Important stuff ## Important stuff
eodf = stim.metadata[stim.name]["EODf"][0][0] eodf = stim.metadata[stim.name]["EODf"][0][0]
stimulus_frequency = eodf + df stimulus_frequency = eodf + df
AM = 50 # Hz AM = 50 # Hz
#print(f"EODf: {eodf}, Stimulus Frequency: {stimulus_frequency}, AM: {AM}") #print(f"EODf: {eodf}, Stimulus Frequency: {stimulus_frequency}, AM: {AM}")
frequencies = [AM, eodf, stimulus_frequency] frequencies = [AM, eodf, stimulus_frequency]
categories = ["AM", "EODf", "Stimulus frequency"] categories = ["AM", "EODf", "Stimulus frequency"]
num_harmonics = [4, 2, 2] num_harmonics = [4, 2, 2]
colors = ["green", "orange", "red"] colors = ["green", "orange", "red"]
delta = 2.5 delta = 2.5
### Peaks im Powerspektrum finden ### ### Peaks im Powerspektrum finden ###
points, color_mapping, points_categories = prepare_harmonics(frequencies, categories, num_harmonics, colors) points, color_mapping, points_categories = prepare_harmonics(frequencies, categories, num_harmonics, colors)
fig = plot_power_spectrum_with_integrals(frequency, power, points, delta, color_mapping, points_categories) fig = plot_power_spectrum_with_integrals(frequency, power, points, delta, color_mapping, points_categories)
plt.show() plt.show()