added data for example cells
This commit is contained in:
383
spectral.py
Normal file
383
spectral.py
Normal file
@@ -0,0 +1,383 @@
|
||||
"""
|
||||
Spectral analysis of neuronal responses.
|
||||
|
||||
## Functions
|
||||
|
||||
- `whitenoise()`: band-limited white noise.
|
||||
- `rate()`: firing rate computed by kernel convolution.
|
||||
- `spectra()`: stimulus- and response power spectra, and cross spectrum.
|
||||
- `susceptibilities()`: stimulus- and response spectra up to second order.
|
||||
- `diag_projection()`: projection of the chi2 matrix onto its diagonal.
|
||||
- `hor_projection()`: horizontal projection of the chi2 matrix.
|
||||
- `peakedness()`: normalized size of a peak expected around a specific frequency.
|
||||
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
from scipy.signal import welch, csd
|
||||
from scipy.stats import norm
|
||||
|
||||
|
||||
def whitenoise(cflow, cfup, dt, duration, rng=np.random.default_rng()):
|
||||
"""Band-limited white noise.
|
||||
|
||||
Generates white noise with a flat power spectrum between `cflow` and
|
||||
`cfup` Hertz, zero mean and unit standard deviation. Note, that in
|
||||
particular for short segments of the generated noise the mean and
|
||||
standard deviation of the returned noise can deviate from zero and
|
||||
one.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
cflow: float
|
||||
Lower cutoff frequency in Hertz.
|
||||
cfup: float
|
||||
Upper cutoff frequency in Hertz.
|
||||
dt: float
|
||||
Time step of the resulting array in seconds.
|
||||
duration: float
|
||||
Total duration of the resulting array in seconds.
|
||||
|
||||
Returns
|
||||
-------
|
||||
noise: 1-D array
|
||||
White noise.
|
||||
"""
|
||||
# number of elements needed for the noise stimulus:
|
||||
n = int(np.ceil((duration+0.5*dt)/dt))
|
||||
# next power of two:
|
||||
nn = int(2**(np.ceil(np.log2(n))))
|
||||
# indices of frequencies with `cflow` and `cfup`:
|
||||
inx0 = int(np.round(dt*nn*cflow))
|
||||
inx1 = int(np.round(dt*nn*cfup))
|
||||
if inx0 < 0:
|
||||
inx0 = 0
|
||||
if inx1 >= nn/2:
|
||||
inx1 = nn/2
|
||||
# draw random numbers in Fourier domain:
|
||||
whitef = np.zeros((nn//2+1), dtype=complex)
|
||||
# zero and nyquist frequency must be real:
|
||||
if inx0 == 0:
|
||||
whitef[0] = 0
|
||||
inx0 = 1
|
||||
if inx1 >= nn//2:
|
||||
whitef[nn//2] = 1
|
||||
inx1 = nn//2-1
|
||||
phases = 2*np.pi*rng.random(size=inx1 - inx0 + 1)
|
||||
whitef[inx0:inx1+1] = np.cos(phases) + 1j*np.sin(phases)
|
||||
# inverse FFT:
|
||||
noise = np.real(np.fft.irfft(whitef))
|
||||
# scaling factor to ensure standard deviation of one:
|
||||
sigma = nn / np.sqrt(2*float(inx1 - inx0))
|
||||
return noise[:n]*sigma
|
||||
|
||||
|
||||
def rate(time, spikes, sigma):
|
||||
""" Firing rate computed by kernel convolution.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
time: ndarray of float
|
||||
Times at which firing rate is evaluated.
|
||||
spikes: list of ndarray of float
|
||||
Spike times.
|
||||
sigma: float
|
||||
Width of the Gaussian kernel as a standard deviation.
|
||||
|
||||
Returns
|
||||
-------
|
||||
rate: ndarray of float
|
||||
Firing rate of convolved spike trains averaged over trials.
|
||||
ratesd: ndarray of float
|
||||
Corresponding standard deviation.
|
||||
"""
|
||||
kernel = norm.pdf(time[time < 8*sigma], loc=4*sigma, scale=sigma)
|
||||
rates = np.zeros((len(spikes), len(time)))
|
||||
xtime = np.append(time, time[-1] + time[1] - time[0])
|
||||
for i, spiket in enumerate(spikes):
|
||||
b, _ = np.histogram(spiket, xtime)
|
||||
rates[i] = np.convolve(b, kernel, 'same')
|
||||
return np.mean(rates, 0), np.std(rates, 0)
|
||||
|
||||
|
||||
def spectra(stimulus, spikes, dt, nfft):
|
||||
"""Stimulus- and response power spectra, and cross spectrum.
|
||||
|
||||
Compute the complex-valued transfer function (first-order
|
||||
susceptibility) and the stimulus-response coherence like this:
|
||||
|
||||
```
|
||||
freqs, pss, prr, prs = spectra(stimulus, spikes, dt, nfft)
|
||||
transfer = prs/pss
|
||||
coherence = np.abs(prs)**2/pss/prr
|
||||
```
|
||||
|
||||
The gain of the transfer function is the absolute value of the
|
||||
transfer function:
|
||||
|
||||
```
|
||||
gain = np.abs(prs)/pss
|
||||
```
|
||||
|
||||
Parameters
|
||||
----------
|
||||
stimulus: ndarray of float
|
||||
Stimulus waveform with sampling interval 'dt'.
|
||||
spikes: list of ndarrays of float
|
||||
Spike times in response to the stimulus.
|
||||
dt: float
|
||||
Sampling interval of stimulus and resolution of the binary spike train.
|
||||
nfft: int
|
||||
Number of samples used for each Fourier transformation.
|
||||
|
||||
Returns
|
||||
-------
|
||||
freqs: ndarray of float
|
||||
The frequencies corresponding to the spectra.
|
||||
pss: ndarray of float
|
||||
Power spectrum of the stimulus.
|
||||
prr: ndarray of float
|
||||
Power spectrum of the response averaged over trials.
|
||||
prs: ndarray of complex
|
||||
Cross spectrum between stimulus and response averaged over trials.
|
||||
|
||||
"""
|
||||
time = np.arange(len(stimulus))*dt
|
||||
freq, pss = welch(stimulus, fs=1/dt, nperseg=nfft, noverlap=nfft//2)
|
||||
prr = np.zeros((len(spikes), len(freq)))
|
||||
prs = np.zeros((len(spikes), len(freq)), dtype=complex)
|
||||
for i, spiket in enumerate(spikes):
|
||||
b, _ = np.histogram(spiket, time)
|
||||
b = b / dt
|
||||
f, rr = welch(b - np.mean(b), fs=1/dt, nperseg=nfft, noverlap=nfft//2)
|
||||
f, rs = csd(b - np.mean(b), stimulus,
|
||||
fs=1/dt, nperseg=nfft, noverlap=nfft//2)
|
||||
prr[i] = rr
|
||||
prs[i] = rs
|
||||
return freq, pss, np.mean(prr, 0), np.mean(prs, 0)
|
||||
|
||||
|
||||
def susceptibilities(stimulus, spikes, dt=0.0005, nfft=2**9):
|
||||
""" Stimulus- and response spectra up to second order.
|
||||
|
||||
Compute the complex-valued transfer function (first-order
|
||||
susceptibility) and the stimulus-response coherence like this:
|
||||
|
||||
```
|
||||
freqs, pss, prr, prs, prss, n = susceptibilities(stimulus, spikes, dt, nfft)
|
||||
transfer = prs/pss
|
||||
coherence = np.abs(prs)**2/pss/prr
|
||||
```
|
||||
|
||||
The gain of the transfer function is the absolute value of the
|
||||
transfer function:
|
||||
|
||||
```
|
||||
gain = np.abs(prs)/pss
|
||||
```
|
||||
|
||||
The second-order susceptibility can be computed like this:
|
||||
|
||||
```
|
||||
chi2 = prss*0.5/np.sqrt(pss.reshape(1, -1)*pss.reshape(-1, 1))
|
||||
```
|
||||
|
||||
Parameters
|
||||
----------
|
||||
stimulus: ndarray of float
|
||||
Stimulus waveform with sampling interval 'dt'.
|
||||
spikes: list of ndarrays of float
|
||||
Spike times in response to the stimulus.
|
||||
dt: float
|
||||
Sampling interval of stimulus and resolution of the binary spike train.
|
||||
nfft: int
|
||||
Number of samples used for each Fourier transformation.
|
||||
|
||||
Returns
|
||||
-------
|
||||
freqs: ndarray of float
|
||||
The frequencies corresponding to the spectra.
|
||||
pss: ndarray of float
|
||||
Power spectrum of the stimulus.
|
||||
prr: ndarray of float
|
||||
Power spectrum of the response averaged over segments.
|
||||
prs: ndarray of complex
|
||||
Cross spectrum between stimulus and response averaged over segments.
|
||||
prss: ndarray of complex
|
||||
Cross bispectrum between stimulus and response averaged over segments.
|
||||
n: int
|
||||
Number of segments.
|
||||
"""
|
||||
freqs = np.fft.fftfreq(nfft, dt)
|
||||
idx = np.argmin(freqs)
|
||||
freqs = np.roll(freqs, -idx)
|
||||
f0 = np.argmin(np.abs(freqs)) # index of zero frequency
|
||||
fidx = np.arange(len(freqs))
|
||||
fsum_idx = fidx.reshape(-1, 1) + fidx.reshape(1, -1) - f0
|
||||
fsum_idx[fsum_idx < 0] = 0
|
||||
fsum_idx[fsum_idx >= len(fidx)] = len(fidx) - 1
|
||||
f0 = len(freqs)//4
|
||||
f1 = 3*len(freqs)//4
|
||||
segments = range(0, len(stimulus) - nfft, nfft)
|
||||
# stimulus:
|
||||
p_ss = np.zeros(len(freqs))
|
||||
fourier_s = np.zeros((len(segments), len(freqs)), complex)
|
||||
n = 0
|
||||
for j, k in enumerate(segments):
|
||||
fourier_s[j] = np.fft.fft(stimulus[k:k + nfft], n=nfft)
|
||||
fourier_s[j] = np.roll(fourier_s[j], -idx)
|
||||
p_ss += np.abs(fourier_s[j]*np.conj(fourier_s[j]))
|
||||
n += 1
|
||||
p_ss /= n
|
||||
# response spectra:
|
||||
time = np.arange(len(stimulus))*dt
|
||||
p_rr = np.zeros(len(freqs))
|
||||
p_rs = np.zeros(len(freqs), complex)
|
||||
p_rss = np.zeros((len(freqs), len(freqs)), complex)
|
||||
n = 0
|
||||
for i, spiket in enumerate(spikes):
|
||||
b, _ = np.histogram(spiket, time)
|
||||
b = b / dt
|
||||
for j, k in enumerate(segments):
|
||||
# stimulus:
|
||||
fourier_s1 = fourier_s[j].reshape(len(fourier_s[j]), 1)
|
||||
fourier_s2 = fourier_s[j].reshape(1, len(fourier_s[j]))
|
||||
fourier_s12 = np.conj(fourier_s1)*np.conj(fourier_s2)
|
||||
# response:
|
||||
fourier_r = np.fft.fft(b[k:k + nfft] - np.mean(b), n=nfft)
|
||||
fourier_r = np.roll(fourier_r, -idx)
|
||||
p_rr += np.abs(fourier_r*np.conj(fourier_r))
|
||||
p_rs += np.conj(fourier_s[j])*fourier_r
|
||||
p_rss += fourier_s12*fourier_r[fsum_idx]
|
||||
n += 1
|
||||
return freqs[f0:f1], p_ss[f0:f1], p_rr[f0:f1]/n, p_rs[f0:f1]/n, p_rss[f0:f1, f0:f1]/n, n
|
||||
|
||||
|
||||
def diag_projection(freqs, chi2, fmax):
|
||||
""" Projection of the chi2 matrix onto its diagonal.
|
||||
|
||||
Adapted from https://stackoverflow.com/questions/71362928/average-values-over-all-offset-diagonals
|
||||
|
||||
Parameters
|
||||
----------
|
||||
freqs: ndarray of float
|
||||
Frequencies of the chi2 matrix.
|
||||
chi2: 2-D ndarray of float
|
||||
Second-order susceptibility matrix.
|
||||
fmax: float
|
||||
Maximum frequency for the projection.
|
||||
|
||||
Returns
|
||||
-------
|
||||
dfreqs: ndarray of float
|
||||
Frequencies of the projection.
|
||||
diagp: ndarray of float
|
||||
Projections of the chi2 matrix onto its diagonal.
|
||||
That is, averages over the anti-diagonals.
|
||||
"""
|
||||
i0 = np.argmin(freqs < 0)
|
||||
i1 = np.argmax(freqs > fmax)
|
||||
if i1 == 0:
|
||||
i1 = len(freqs)
|
||||
chi2 = chi2[i0:i1, i0:i1]
|
||||
n = chi2.shape[0]
|
||||
diagp = np.zeros(n*2-1, dtype=float)
|
||||
for i in range(n):
|
||||
diagp[i:i + n] += chi2[i]
|
||||
diagp[0:n] /= np.arange(1, n+1, 1, dtype=float)
|
||||
diagp[n:] /= np.arange(n-1, 0, -1, dtype=float)
|
||||
dfreqs = np.arange(len(diagp))*(freqs[i0 + 1] - freqs[i0]) + freqs[i0]
|
||||
return dfreqs, diagp
|
||||
|
||||
|
||||
def hor_projection(freqs, chi2, fmax):
|
||||
""" Horizontal projection of the chi2 matrix.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
freqs: ndarray of float
|
||||
Frequencies of the chi2 matrix.
|
||||
chi2: 2-D ndarray of float
|
||||
Second-order susceptibility matrix.
|
||||
fmax: float
|
||||
Maximum frequency for the projection.
|
||||
|
||||
Returns
|
||||
-------
|
||||
hfreqs: ndarray of float
|
||||
Frequencies of the projection.
|
||||
horp: ndarray of float
|
||||
Projections of the chi2 matrix onto its x-axis.
|
||||
That is, averages over columns.
|
||||
"""
|
||||
i0 = np.argmin(freqs < 0)
|
||||
i1 = np.argmax(freqs > fmax)
|
||||
if i1 == 0:
|
||||
i1 = len(freqs)
|
||||
hfreqs = freqs[i0:i1]
|
||||
chi2 = chi2[i0:i1, i0:i1]
|
||||
horp = np.mean(chi2, 1)
|
||||
return hfreqs, horp
|
||||
|
||||
|
||||
def peakedness(freqs, projection, fbase, median=True,
|
||||
searchwin=50, averagewin=10):
|
||||
"""Normalized size of a peak expected around a specific frequency.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
freqs: ndarray of float
|
||||
Frequencies of the projection.
|
||||
projection: ndarray of float
|
||||
Projection of the chi2 matrix.
|
||||
fbase: float
|
||||
The neurons baseline-frequency.
|
||||
That is, the frequency where a peak is expected in the projection.
|
||||
median: bool
|
||||
If True, normalize the peak height by the median of the projection.
|
||||
Otherwise (default), normalize by averaged values of the projection
|
||||
close to the fbase frequency.
|
||||
searchwin: float
|
||||
Search for peak in the projection at fbase plus and
|
||||
minus `searchwin` Hertz.
|
||||
averagewin: float
|
||||
For estimating the reference level, an average is taken in two
|
||||
`averagewin` Hertz wide windows that are located `averagewin`
|
||||
Hertz to the left and right of the detected peak.
|
||||
|
||||
Returns
|
||||
-------
|
||||
p: float
|
||||
The normalized height of the peak close to fbase.
|
||||
fpeak: float
|
||||
The frequency of the detected peak.
|
||||
|
||||
"""
|
||||
sel = (freqs > fbase - searchwin) & (freqs < fbase + searchwin)
|
||||
snippet = projection[sel]
|
||||
if len(snippet) == 0:
|
||||
return np.nan, np.nan
|
||||
peak = np.max(snippet)
|
||||
fpeak = freqs[np.argmax(snippet) + np.argmax(sel)]
|
||||
bleft = np.nan
|
||||
bright = np.nan
|
||||
if median:
|
||||
baseline = np.median(projection)
|
||||
else:
|
||||
mask = (freqs >= fpeak - 2*averagewin) & (freqs <= fpeak - averagewin)
|
||||
bleft = np.mean(projection[mask]) if np.sum(mask) > 0 else np.nan
|
||||
mask = (freqs >= fpeak + averagewin) & (freqs <= fpeak + 2*averagewin)
|
||||
bright = np.mean(projection[mask]) if np.sum(mask) > 0 else np.nan
|
||||
if np.isfinite(bleft) and np.isfinite(bright):
|
||||
baseline = 0.5*(bleft + bright)
|
||||
elif np.isfinite(bleft):
|
||||
baseline = bleft
|
||||
elif np.isfinite(bright):
|
||||
baseline = bright
|
||||
else:
|
||||
baseline = np.nan
|
||||
if np.isnan(peak/baseline):
|
||||
print(peak, fpeak, fbase, baseline, bleft, bright, median)
|
||||
return np.nan, np.nan
|
||||
return peak/baseline, fpeak
|
||||
|
||||
Reference in New Issue
Block a user