Merge pull request #1 from jdpigeon/master

Updated lsl-viewer and minor updates to lab7

lab4/lsl-viewer.py

@@ -1,17 +1,15 @@
 #!/usr/bin/env python
-## code by Alexandre Barachant
 import numpy as np
 import matplotlib.pyplot as plt
-from scipy.signal import butter, filtfilt
+from scipy.signal import butter, lfilter, lfilter_zi
-from time import time, sleep
+from time import sleep
 from pylsl import StreamInlet, resolve_byprop
+from optparse import OptionParser
 import seaborn as sns
 from threading import Thread
 
 sns.set(style="whitegrid")
 
-from optparse import OptionParser
-
 parser = OptionParser()
 
 parser.add_option("-w", "--window",
@@ -55,6 +53,7 @@ class LSLViewer():
         self.dejitter = dejitter
         self.inlet = StreamInlet(stream, max_chunklen=buf)
         self.filt = True
+
         info = self.inlet.info()
         description = info.desc()
 
@@ -105,6 +104,9 @@ class LSLViewer():
 
         self.bf, self.af = butter(4, np.array([1, 40])/(self.sfreq/2.),
                                   'bandpass')
+        zi = lfilter_zi(self.bf, self.af)
+        self.filt_state = np.tile(zi, (self.n_chan, 1)).transpose()
+        self.data_f = np.zeros((self.n_samples, self.n_chan))
 
     def update_plot(self):
         k = 0
@@ -112,33 +114,37 @@ class LSLViewer():
             samples, timestamps = self.inlet.pull_chunk(timeout=1.0,
                                                         max_samples=12)
             if timestamps:
-                self.data = np.vstack([self.data, samples])
                 if self.dejitter:
                     timestamps = np.float64(np.arange(len(timestamps)))
                     timestamps /= self.sfreq
                     timestamps += self.times[-1] + 1./self.sfreq
                 self.times = np.concatenate([self.times, timestamps])
-
                 self.n_samples = int(self.sfreq * self.window)
-                self.data = self.data[-self.n_samples:]
                 self.times = self.times[-self.n_samples:]
+                self.data = np.vstack([self.data, samples])
+                self.data = self.data[-self.n_samples:]
+                filt_samples, self.filt_state = lfilter(
+                    self.bf, self.af,
+                    samples,
+                    axis=0, zi=self.filt_state)
+                self.data_f = np.vstack([self.data_f, filt_samples])
+                self.data_f = self.data_f[-self.n_samples:]
                 k += 1
                 if k == self.display_every:
-                    if self.filt:
-                        data_f = filtfilt(self.bf, self.af, self.data, axis=0)
-                    else:
-                        data_f = self.data
-                        data_f -= data_f.mean(axis=0)
 
+                    if self.filt:
+                        plot_data = self.data_f
+                    elif not self.filt:
+                        plot_data = self.data - self.data.mean(axis=0)
                     for ii in range(self.n_chan):
                         self.lines[ii].set_xdata(self.times[::subsample] -
                                                  self.times[-1])
-                        self.lines[ii].set_ydata(data_f[::subsample, ii] /
+                        self.lines[ii].set_ydata(plot_data[::subsample, ii] /
                                                  self.scale - ii)
+                        impedances = np.std(plot_data, axis=0)
 
-                    impedances = np.std(data_f, axis=0)
+                    ticks_labels = ['%s - %.2f' % (self.ch_names[ii],
-                    ticks_labels = ['%s - %.2f' %
+                                                   impedances[ii])
-                                    (self.ch_names[ii], impedances[ii])
                                     for ii in range(self.n_chan)]
                     self.axes.set_yticklabels(ticks_labels)
                     self.axes.set_xlim(-self.window, 0)

lab5/lsl-viewer.py

@@ -1,17 +1,15 @@
 #!/usr/bin/env python
-## code by Alexandre Barachant
 import numpy as np
 import matplotlib.pyplot as plt
-from scipy.signal import butter, filtfilt
+from scipy.signal import butter, lfilter, lfilter_zi
-from time import time, sleep
+from time import sleep
 from pylsl import StreamInlet, resolve_byprop
+from optparse import OptionParser
 import seaborn as sns
 from threading import Thread
 
 sns.set(style="whitegrid")
 
-from optparse import OptionParser
-
 parser = OptionParser()
 
 parser.add_option("-w", "--window",
@@ -55,6 +53,7 @@ class LSLViewer():
         self.dejitter = dejitter
         self.inlet = StreamInlet(stream, max_chunklen=buf)
         self.filt = True
+
         info = self.inlet.info()
         description = info.desc()
 
@@ -105,6 +104,9 @@ class LSLViewer():
 
         self.bf, self.af = butter(4, np.array([1, 40])/(self.sfreq/2.),
                                   'bandpass')
+        zi = lfilter_zi(self.bf, self.af)
+        self.filt_state = np.tile(zi, (self.n_chan, 1)).transpose()
+        self.data_f = np.zeros((self.n_samples, self.n_chan))
 
     def update_plot(self):
         k = 0
@@ -112,33 +114,37 @@ class LSLViewer():
             samples, timestamps = self.inlet.pull_chunk(timeout=1.0,
                                                         max_samples=12)
             if timestamps:
-                self.data = np.vstack([self.data, samples])
                 if self.dejitter:
                     timestamps = np.float64(np.arange(len(timestamps)))
                     timestamps /= self.sfreq
                     timestamps += self.times[-1] + 1./self.sfreq
                 self.times = np.concatenate([self.times, timestamps])
-
                 self.n_samples = int(self.sfreq * self.window)
-                self.data = self.data[-self.n_samples:]
                 self.times = self.times[-self.n_samples:]
+                self.data = np.vstack([self.data, samples])
+                self.data = self.data[-self.n_samples:]
+                filt_samples, self.filt_state = lfilter(
+                    self.bf, self.af,
+                    samples,
+                    axis=0, zi=self.filt_state)
+                self.data_f = np.vstack([self.data_f, filt_samples])
+                self.data_f = self.data_f[-self.n_samples:]
                 k += 1
                 if k == self.display_every:
-                    if self.filt:
-                        data_f = filtfilt(self.bf, self.af, self.data, axis=0)
-                    else:
-                        data_f = self.data
-                        data_f -= data_f.mean(axis=0)
 
+                    if self.filt:
+                        plot_data = self.data_f
+                    elif not self.filt:
+                        plot_data = self.data - self.data.mean(axis=0)
                     for ii in range(self.n_chan):
                         self.lines[ii].set_xdata(self.times[::subsample] -
                                                  self.times[-1])
-                        self.lines[ii].set_ydata(data_f[::subsample, ii] /
+                        self.lines[ii].set_ydata(plot_data[::subsample, ii] /
                                                  self.scale - ii)
+                        impedances = np.std(plot_data, axis=0)
 
-                    impedances = np.std(data_f, axis=0)
+                    ticks_labels = ['%s - %.2f' % (self.ch_names[ii],
-                    ticks_labels = ['%s - %.2f' %
+                                                   impedances[ii])
-                                    (self.ch_names[ii], impedances[ii])
                                     for ii in range(self.n_chan)]
                     self.axes.set_yticklabels(ticks_labels)
                     self.axes.set_xlim(-self.window, 0)

lab6/lsl-viewer.py

@@ -1,17 +1,15 @@
 #!/usr/bin/env python
-## code by Alexandre Barachant
 import numpy as np
 import matplotlib.pyplot as plt
-from scipy.signal import butter, filtfilt
+from scipy.signal import butter, lfilter, lfilter_zi
-from time import time, sleep
+from time import sleep
 from pylsl import StreamInlet, resolve_byprop
+from optparse import OptionParser
 import seaborn as sns
 from threading import Thread
 
 sns.set(style="whitegrid")
 
-from optparse import OptionParser
-
 parser = OptionParser()
 
 parser.add_option("-w", "--window",
@@ -55,6 +53,7 @@ class LSLViewer():
         self.dejitter = dejitter
         self.inlet = StreamInlet(stream, max_chunklen=buf)
         self.filt = True
+
         info = self.inlet.info()
         description = info.desc()
 
@@ -105,6 +104,9 @@ class LSLViewer():
 
         self.bf, self.af = butter(4, np.array([1, 40])/(self.sfreq/2.),
                                   'bandpass')
+        zi = lfilter_zi(self.bf, self.af)
+        self.filt_state = np.tile(zi, (self.n_chan, 1)).transpose()
+        self.data_f = np.zeros((self.n_samples, self.n_chan))
 
     def update_plot(self):
         k = 0
@@ -112,33 +114,37 @@ class LSLViewer():
             samples, timestamps = self.inlet.pull_chunk(timeout=1.0,
                                                         max_samples=12)
             if timestamps:
-                self.data = np.vstack([self.data, samples])
                 if self.dejitter:
                     timestamps = np.float64(np.arange(len(timestamps)))
                     timestamps /= self.sfreq
                     timestamps += self.times[-1] + 1./self.sfreq
                 self.times = np.concatenate([self.times, timestamps])
-
                 self.n_samples = int(self.sfreq * self.window)
-                self.data = self.data[-self.n_samples:]
                 self.times = self.times[-self.n_samples:]
+                self.data = np.vstack([self.data, samples])
+                self.data = self.data[-self.n_samples:]
+                filt_samples, self.filt_state = lfilter(
+                    self.bf, self.af,
+                    samples,
+                    axis=0, zi=self.filt_state)
+                self.data_f = np.vstack([self.data_f, filt_samples])
+                self.data_f = self.data_f[-self.n_samples:]
                 k += 1
                 if k == self.display_every:
-                    if self.filt:
-                        data_f = filtfilt(self.bf, self.af, self.data, axis=0)
-                    else:
-                        data_f = self.data
-                        data_f -= data_f.mean(axis=0)
 
+                    if self.filt:
+                        plot_data = self.data_f
+                    elif not self.filt:
+                        plot_data = self.data - self.data.mean(axis=0)
                     for ii in range(self.n_chan):
                         self.lines[ii].set_xdata(self.times[::subsample] -
                                                  self.times[-1])
-                        self.lines[ii].set_ydata(data_f[::subsample, ii] /
+                        self.lines[ii].set_ydata(plot_data[::subsample, ii] /
                                                  self.scale - ii)
+                        impedances = np.std(plot_data, axis=0)
 
-                    impedances = np.std(data_f, axis=0)
+                    ticks_labels = ['%s - %.2f' % (self.ch_names[ii],
-                    ticks_labels = ['%s - %.2f' %
+                                                   impedances[ii])
-                                    (self.ch_names[ii], impedances[ii])
                                     for ii in range(self.n_chan)]
                     self.axes.set_yticklabels(ticks_labels)
                     self.axes.set_xlim(-self.window, 0)

lab7/analysis/analysis.py

@@ -0,0 +1,67 @@
+from scipy import signal, stats
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import argparse
+import sys
+
+## sampling rate of ganglion is 200Hz
+fs = 200.0
+
+def extract_data(fname):
+    d = pd.read_csv(fname)
+    eeg = np.array(d.ix[:, 1:5])
+    tag = np.array(d.Marker)
+    start = np.where(tag == 1)[0][0]
+    end = np.where(tag == 2)[0][0]
+    eeg = eeg[start:end]
+    return eeg
+
+## get data
+eeg1 = extract_data(sys.argv[1])
+eeg2 = extract_data(sys.argv[2])
+
+## data may be off by a few samples, so need to align
+N_samp = min(len(eeg1), len(eeg2))
+eeg1 = eeg1[:N_samp]
+eeg2 = eeg2[:N_samp]
+
+## filter signal to remove noise
+b, a = signal.butter(2, (2/(fs/2), 20/(fs/2)), btype='bandpass')
+eeg1 = signal.filtfilt(b, a, eeg1, axis=0)
+eeg2 = signal.filtfilt(b, a, eeg2, axis=0)
+
+## advance window of 200 samples
+## take correlation between signals across each sample
+window = 200
+step = 25
+corr = []
+times = []
+
+for start in np.arange(0, N_samp, step):
+    end = start + window
+    w1 = eeg1[start:end]
+    w2 = eeg2[start:end]
+
+    ## average the correlation across each channel
+    r = 0
+    for c in range(w1.shape[1]):
+        r += stats.pearsonr(w1[:, c], w2[:, c])[0]
+    r /= w1.shape[1]
+
+    mid = (start+end)/2 # middle sample
+    t = mid / fs # convert middle sample to time
+
+
+    times.append(t)
+    corr.append(r)
+
+times = np.array(times)
+corr = np.array(corr)
+
+plt.figure(figsize=(14,4))
+plt.plot(times, np.abs(corr))
+plt.xlabel('Time (s)')
+plt.ylabel('Correlation')
+plt.title('Correlation across brains')
+plt.show()

lab7/analysis/cross_brain_correlation.ipynb

@@ -76,7 +76,9 @@
   {
    "cell_type": "code",
    "execution_count": 21,
-   "metadata": {},
+   "metadata": {
+    "collapsed": true
+   },
    "outputs": [],
    "source": [
     "## filter signal to remove noise\n",
@@ -181,7 +183,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.0"
+   "version": "3.6.1"
   }
  },
  "nbformat": 4,

lab7/recording/lsl-viewer.py

@@ -1,17 +1,15 @@
 #!/usr/bin/env python
-## code by Alexandre Barachant
 import numpy as np
 import matplotlib.pyplot as plt
-from scipy.signal import butter, filtfilt
+from scipy.signal import butter, lfilter, lfilter_zi
-from time import time, sleep
+from time import sleep
 from pylsl import StreamInlet, resolve_byprop
+from optparse import OptionParser
 import seaborn as sns
 from threading import Thread
 
 sns.set(style="whitegrid")
 
-from optparse import OptionParser
-
 parser = OptionParser()
 
 parser.add_option("-w", "--window",
@@ -55,6 +53,7 @@ class LSLViewer():
         self.dejitter = dejitter
         self.inlet = StreamInlet(stream, max_chunklen=buf)
         self.filt = True
+
         info = self.inlet.info()
         description = info.desc()
 
@@ -105,6 +104,9 @@ class LSLViewer():
 
         self.bf, self.af = butter(4, np.array([1, 40])/(self.sfreq/2.),
                                   'bandpass')
+        zi = lfilter_zi(self.bf, self.af)
+        self.filt_state = np.tile(zi, (self.n_chan, 1)).transpose()
+        self.data_f = np.zeros((self.n_samples, self.n_chan))
 
     def update_plot(self):
         k = 0
@@ -112,33 +114,37 @@ class LSLViewer():
             samples, timestamps = self.inlet.pull_chunk(timeout=1.0,
                                                         max_samples=12)
             if timestamps:
-                self.data = np.vstack([self.data, samples])
                 if self.dejitter:
                     timestamps = np.float64(np.arange(len(timestamps)))
                     timestamps /= self.sfreq
                     timestamps += self.times[-1] + 1./self.sfreq
                 self.times = np.concatenate([self.times, timestamps])
-
                 self.n_samples = int(self.sfreq * self.window)
-                self.data = self.data[-self.n_samples:]
                 self.times = self.times[-self.n_samples:]
+                self.data = np.vstack([self.data, samples])
+                self.data = self.data[-self.n_samples:]
+                filt_samples, self.filt_state = lfilter(
+                    self.bf, self.af,
+                    samples,
+                    axis=0, zi=self.filt_state)
+                self.data_f = np.vstack([self.data_f, filt_samples])
+                self.data_f = self.data_f[-self.n_samples:]
                 k += 1
                 if k == self.display_every:
-                    if self.filt:
-                        data_f = filtfilt(self.bf, self.af, self.data, axis=0)
-                    else:
-                        data_f = self.data
-                        data_f -= data_f.mean(axis=0)
 
+                    if self.filt:
+                        plot_data = self.data_f
+                    elif not self.filt:
+                        plot_data = self.data - self.data.mean(axis=0)
                     for ii in range(self.n_chan):
                         self.lines[ii].set_xdata(self.times[::subsample] -
                                                  self.times[-1])
-                        self.lines[ii].set_ydata(data_f[::subsample, ii] /
+                        self.lines[ii].set_ydata(plot_data[::subsample, ii] /
                                                  self.scale - ii)
+                        impedances = np.std(plot_data, axis=0)
 
-                    impedances = np.std(data_f, axis=0)
+                    ticks_labels = ['%s - %.2f' % (self.ch_names[ii],
-                    ticks_labels = ['%s - %.2f' %
+                                                   impedances[ii])
-                                    (self.ch_names[ii], impedances[ii])
                                     for ii in range(self.n_chan)]
                     self.axes.set_yticklabels(ticks_labels)
                     self.axes.set_xlim(-self.window, 0)

lab9/lsl-viewer.py

@@ -1,17 +1,15 @@
 #!/usr/bin/env python
-## code by Alexandre Barachant
 import numpy as np
 import matplotlib.pyplot as plt
-from scipy.signal import butter, filtfilt
+from scipy.signal import butter, lfilter, lfilter_zi
-from time import time, sleep
+from time import sleep
 from pylsl import StreamInlet, resolve_byprop
+from optparse import OptionParser
 import seaborn as sns
 from threading import Thread
 
 sns.set(style="whitegrid")
 
-from optparse import OptionParser
-
 parser = OptionParser()
 
 parser.add_option("-w", "--window",
@@ -55,6 +53,7 @@ class LSLViewer():
         self.dejitter = dejitter
         self.inlet = StreamInlet(stream, max_chunklen=buf)
         self.filt = True
+
         info = self.inlet.info()
         description = info.desc()
 
@@ -105,6 +104,9 @@ class LSLViewer():
 
         self.bf, self.af = butter(4, np.array([1, 40])/(self.sfreq/2.),
                                   'bandpass')
+        zi = lfilter_zi(self.bf, self.af)
+        self.filt_state = np.tile(zi, (self.n_chan, 1)).transpose()
+        self.data_f = np.zeros((self.n_samples, self.n_chan))
 
     def update_plot(self):
         k = 0
@@ -112,33 +114,37 @@ class LSLViewer():
             samples, timestamps = self.inlet.pull_chunk(timeout=1.0,
                                                         max_samples=12)
             if timestamps:
-                self.data = np.vstack([self.data, samples])
                 if self.dejitter:
                     timestamps = np.float64(np.arange(len(timestamps)))
                     timestamps /= self.sfreq
                     timestamps += self.times[-1] + 1./self.sfreq
                 self.times = np.concatenate([self.times, timestamps])
-
                 self.n_samples = int(self.sfreq * self.window)
-                self.data = self.data[-self.n_samples:]
                 self.times = self.times[-self.n_samples:]
+                self.data = np.vstack([self.data, samples])
+                self.data = self.data[-self.n_samples:]
+                filt_samples, self.filt_state = lfilter(
+                    self.bf, self.af,
+                    samples,
+                    axis=0, zi=self.filt_state)
+                self.data_f = np.vstack([self.data_f, filt_samples])
+                self.data_f = self.data_f[-self.n_samples:]
                 k += 1
                 if k == self.display_every:
-                    if self.filt:
-                        data_f = filtfilt(self.bf, self.af, self.data, axis=0)
-                    else:
-                        data_f = self.data
-                        data_f -= data_f.mean(axis=0)
 
+                    if self.filt:
+                        plot_data = self.data_f
+                    elif not self.filt:
+                        plot_data = self.data - self.data.mean(axis=0)
                     for ii in range(self.n_chan):
                         self.lines[ii].set_xdata(self.times[::subsample] -
                                                  self.times[-1])
-                        self.lines[ii].set_ydata(data_f[::subsample, ii] /
+                        self.lines[ii].set_ydata(plot_data[::subsample, ii] /
                                                  self.scale - ii)
+                        impedances = np.std(plot_data, axis=0)
 
-                    impedances = np.std(data_f, axis=0)
+                    ticks_labels = ['%s - %.2f' % (self.ch_names[ii],
-                    ticks_labels = ['%s - %.2f' %
+                                                   impedances[ii])
-                                    (self.ch_names[ii], impedances[ii])
                                     for ii in range(self.n_chan)]
                     self.axes.set_yticklabels(ticks_labels)
                     self.axes.set_xlim(-self.window, 0)