1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
|
import numpy as np
from nitime.lazy import scipy_stats as stats
from nitime import descriptors as desc
from nitime import algorithms as tsa
from nitime import timeseries as ts
from nitime.index_utils import tril_indices_from
from .base import BaseAnalyzer
def signal_noise(response):
"""
Signal and noise as defined in Borst and Theunissen 1999, Figure 2
Parameters
----------
response: nitime TimeSeries object
The data here are individual responses of a single unit to the same
stimulus, with repetitions being the first dimension and time as the
last dimension
"""
signal = np.mean(response.data, 0) # The estimate of the signal is the
# average response
noise = response.data - signal # Noise is the individual
# repetition's deviation from the
# estimate of the signal
# Return TimeSeries objects with the sampling rate of the input:
return (ts.TimeSeries(signal, sampling_rate=response.sampling_rate),
ts.TimeSeries(noise, sampling_rate=response.sampling_rate))
class SNRAnalyzer(BaseAnalyzer):
"""
Calculate SNR for a response to repetitions of the same stimulus, according
to (Borst, 1999) (Figure 2) and (Hsu, 2004).
Hsu A, Borst A and Theunissen, FE (2004) Quantifying variability in neural
responses ans its application for the validation of model
predictions. Network: Comput Neural Syst 15:91-109
Borst A and Theunissen FE (1999) Information theory and neural coding. Nat
Neurosci 2:947-957
"""
def __init__(self, input=None, bandwidth=None, adaptive=False,
low_bias=False):
"""
Initializer for the multi_taper_SNR object
Parameters
----------
input: TimeSeries object
bandwidth: float,
The bandwidth of the windowing function will determine the number
tapers to use. This parameters represents trade-off between
frequency resolution (lower main lobe bandwidth for the taper) and
variance reduction (higher bandwidth and number of averaged
estimates). Per default will be set to 4 times the fundamental
frequency, such that NW=4
adaptive: bool, default to False
Whether to set the weights for the tapered spectra according to the
adaptive algorithm (Thompson, 2007).
low_bias : bool, default to False
Rather than use 2NW tapers, only use the tapers that have better
than 90% spectral concentration within the bandwidth (still using a
maximum of 2NW tapers)
Notes
-----
Thompson, DJ (2007) Jackknifing multitaper spectrum estimates. IEEE
Signal Processing Magazing. 24: 20-30
"""
self.input = input
self.signal, self.noise = signal_noise(input)
self.bandwidth = bandwidth
self.adaptive = adaptive
self.low_bias = low_bias
@desc.setattr_on_read
def mt_frequencies(self):
return np.linspace(0, self.input.sampling_rate / 2,
self.input.data.shape[-1] / 2 + 1)
@desc.setattr_on_read
def mt_signal_psd(self):
_, p, _ = tsa.multi_taper_psd(self.signal.data,
Fs=self.input.sampling_rate,
BW=self.bandwidth,
adaptive=self.adaptive,
low_bias=self.low_bias)
return p
@desc.setattr_on_read
def mt_noise_psd(self):
p = np.empty((self.noise.data.shape[0],
self.noise.data.shape[-1] / 2 + 1))
for i in range(p.shape[0]):
_, p[i], _ = tsa.multi_taper_psd(self.noise.data[i],
Fs=self.input.sampling_rate,
BW=self.bandwidth,
adaptive=self.adaptive,
low_bias=self.low_bias)
return np.mean(p, 0)
@desc.setattr_on_read
def mt_coherence(self):
""" """
return self.mt_signal_psd / (self.mt_signal_psd + self.mt_noise_psd)
@desc.setattr_on_read
def mt_information(self):
df = self.mt_frequencies[1] - self.mt_frequencies[0]
return -1 * np.log2(1 - self.mt_coherence) * df
#These two formulations should be equivalent
#return np.log2(1+self.mt_snr)
@desc.setattr_on_read
def mt_snr(self):
return self.mt_signal_psd / self.mt_noise_psd
@desc.setattr_on_read
def correlation(self):
"""
The correlation between all combinations of trials
Returns
-------
(r,e) : tuple
r is the mean correlation and e is the mean error of the correlation
(with df = n_trials - 1)
"""
c = np.corrcoef(self.input.data)
c = c[tril_indices_from(c, -1)]
return np.mean(c), stats.sem(c)
|
