Enh trigger

.gitignore

@@ -68,5 +68,7 @@ stage/
 
 *.jar
 spft/output/xxx_Test_Right_Only_out-file.yml
-spft/__pycache__/spft.cpython-39.pyc
-spft/__pycache__/spft.cpython-39.pyc
+spft/__pycache__/*
+spft/__pycache__
+spft/testing/output/*.yml
+.ipynb_checkpoints*

.vscode/settings.json

@@ -0,0 +1,3 @@
+{
+    "java.configuration.updateBuildConfiguration": "interactive"
+}

README.md

@@ -8,6 +8,12 @@ This task is modeled on previous work, including the following (most relevant) r
 - https://pubmed.ncbi.nlm.nih.gov/34704176/
 - https://pubmed.ncbi.nlm.nih.gov/33885965/
 
+### Integration (triggering)
+- The session can be triggered to start by keypress values defined in the [configuration file](#configuration-file)
+- Each trial start/stop can also be synchronized with external devices through serial port triggering (or emulated serial port) 
+  - a single trigger message is sent to the external device (e.g, arduino) for both start and stop: `byte[] TRIGGER_MESSAGE = {'1','\n'}`
+  - the generated [output file](#output) will contain a triggersOut field with times and values indicating start (`'0'`) and stop (`'1'`) 
+
 ## Build
 The application is built with Maven 3
 
@@ -22,7 +28,9 @@ This will generate a runnable jar in the `target` folder
 ## Run
 This application uses Java 17. Make sure it is installed before running the application.
 
-- if you are running this on a linux system, you will require explicit access to the incoming data from ports {`/dev/ttyACM0`; `/dev/ttyACM1`}. This can be granted by the _system administrator_ providing the user with access to the `dialout` group with `sudo usermod -a -G dialout theUserNameHere`
+- if you are running this on a linux system, you will require explicit access to the incoming data from ports {`/dev/ttyACM0`; `/dev/ttyACM1`} and to send triggers out from the program via the serial port if desired {`/dev/ttyACM2`}. This can be granted by the _system administrator_ providing the user with access to the `dialout` group with `sudo usermod -a -G dialout theUserNameHere`
+- to ensure that you are recording from the correct devices __always__ plug the devices in to the USB ports in the same order `(left, right)` and you **must** attach the trigger output device last 
+  - this allows you to ensure that you are always aware of which device is which
 
 ### Installing Java17 (linux)
 Java 11 may still be the default for Ubuntu-based OSs, so you must install it explicitly
@@ -39,7 +47,7 @@ Once a JDK or JRE is installed, run
 
 For example, if running a locally built version
 
-`java -jar target/spft-1.0-SNAPSHOT-jar-with-dependencies.jar`
+`java -jar target/spft-1.3-SNAPSHOT-jar-with-dependencies.jar`
 
 This will start a window where you can load a [configuration file](#configuration-file) for a session.
 There is a list of [runtime flags](#runtime-flags) that let you customize some aspects of the application
@@ -72,12 +80,13 @@ If the parameter is missing, the session will start as soon as the experimenter
 `blocks.trials`: A list of references to the sequences that define each trial. The referenced sequence has to exist in
 the top-level pool of sequences  
 `sequences`: A pool of sequences that can be referenced in blocks' trials. Each sequence has the values for the reference
-bar and a frequency in Hz that defines the speed of the sequence 
+bar and a frequency in milliseconds that defines the speed of the sequence. Each value of the sequence is displayed for this length of time.
+  - a sequence must be made up of at least two timepoints (i.e., two height values) to ensure correct timing of feedback at the end of the block
 
 ### Runtime flags
 Runtime flags are JVM System Properties that control some behaviour of the application
 To pass a system property, use the `-Dproperty=value` syntax before the jar name that's standard in java applications. For
-example, to enable debug mode: `java -jar -Ddebug=true target/spft-1.0-SNAPSHOT-jar-with-dependencies.jar`
+example, to enable debug mode: `java -jar -Ddebug=true target/spft-1.3-SNAPSHOT-jar-with-dependencies.jar`
 
 `debug`: Enables extra verbose logging and panel coloured backgrounds to see where each panel ends  
 `spft.forceData.smoothWindowSize`: The number of samples in the averaging window to calculate the height of the force bar. Default: 1  
@@ -117,7 +126,8 @@ relative to the other times in the session
 `blocks.trials`: A list of the presentation values and their actual timestamps using a CPU clock  
 `blocks.endTimestamp`: The end of the block using a CPU clock  
 `devices`: A list of hardware force devices with each element containing the full stream of data starting when one of the triggers is received  
-`triggers`: A lif of triggers received throughout the session
+`triggers`: A list of triggers received throughout the session
+`triggersOut`: If present, indicates serial out hardware device and port name for triggers sent to arduino. Contains the times and values for each trigger
 
 ## Implementation notes
 ### Force bar height

pom.xml

@@ -6,9 +6,9 @@
 
     <groupId>com.github.neuralabc.spft</groupId>
     <artifactId>spft</artifactId>
-    <version>1.1-SNAPSHOT</version>
+    <version>1.4-SNAPSHOT</version>
     <name>Sequence Pinch Force Task</name>
-    <description>A Java implementation of a Sequence Pinch Force task</description>
+    <description>A Java implementation of a Sequence Pinch Force task, with ability to trigger to external serial device</description>
 
     <properties>
         <project.build.sourceEncoding>UTF-8</project.build.sourceEncoding>

spft/2022_06_Trial_by_trial_Bx.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [
     {
@@ -11,6 +11,14 @@
      "text": [
       "Populating the interactive namespace from numpy and matplotlib\n"
      ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/opt/quarantine/miniforge/envs/py3p9/lib/python3.9/site-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.16.5 and <1.23.0 is required for this version of SciPy (detected version 1.23.3\n",
+      "  warnings.warn(f\"A NumPy version >={np_minversion} and <{np_maxversion}\"\n"
+     ]
     }
    ],
    "source": [
@@ -25,7 +33,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [
     {

spft/config_files/force_pyramid.yml

@@ -0,0 +1,25 @@
+sessionName: force_pyramid
+outputSuffix: force_pyramid-file.yml
+interBlockInterval: 2000
+triggers: [ "5", "7" ] # optional. If missing, the task will start as soon as the "Start" button is clicked
+colours: # optional
+  leftReference: 1D8348
+  leftForce: 82E0AA
+forceProportionRange:
+  min: 0.00
+  max: 1.00
+blocks:
+  - name: "Block1"
+    instructions: "Block1 starting"
+    instructionsDuration: 2000
+    feedback: "Block1 finished"
+    feedbackDuration: 2000
+    interTrialInterval: 1000
+    trials:
+      - name: "1_pyramid"
+        sequenceRef: "pyramid"
+
+sequences:
+  pyramid:
+    frequency: 3000
+    valuesLeft: [.05,.2,.4,.6,.6,.4,.2,.05]

spft/config_files/maximum_voluntary_contraction.yml

@@ -0,0 +1,25 @@
+sessionName: maximum_voluntary_contraction
+outputSuffix: mvc-file.yml
+interBlockInterval: 2000
+triggers: [ "5", "7" ] # optional. If missing, the task will start as soon as the "Start" button is clicked
+colours: # optional
+  leftReference: 1D8348
+  leftForce: 82E0AA
+forceProportionRange:
+  min: 0.00
+  max: 1.00
+blocks:
+  - name: "Block1"
+    instructions: "Pinch as hard as you can"
+    instructionsDuration: 5000
+    feedback: "Done, thank you"
+    feedbackDuration: 5000
+    interTrialInterval: 1000
+    trials:
+      - name: "1_mvc"
+        sequenceRef: "mvc"
+
+sequences:
+  mvc:
+    frequency: 1000
+    valuesLeft: [.05,.05,.05,.05,.05]

spft/scripts/compute_CMC.m

@@ -0,0 +1,146 @@
+resample_freq = 500; %need to know the resample frequency, since this determines where our times are
+triggers = [EEG.event.latency];
+triggers = triggers(2:end); %skip the first one, since it appears to be a dummy value
+
+%check if we have the correct number of triggers
+num_triggers = 3*3*5*2; %nine 3 trials for each block, 3 sequences (LRN,SMP,RST), 2 events per trial (start,stop)
+
+%if this is true, then we are good to go!
+numel(triggers) == num_triggers
+
+%
+triggers = triggers./resample_freq*1000; %to put into same time space as "times", divide by freq and convert to ms
+times = EEG.times;
+
+trigger_idxs = [];
+for i = 1:length(triggers)
+    index = find(times >=triggers(i) ,1);
+    trigger_idxs = [trigger_idxs,index];
+end
+
+
+trigger_idxs_on = trigger_idxs(1:2:end);
+trigger_idxs_off = trigger_idxs(2:2:end);
+
+%assuming LRN, SMP, RST, with 3 trials each and 5 blocks
+LRN_idx= [1:3,[1:3]+9*1,[1:3]+9*2,[1:3]+9*3,[1:3]+9*4]; 
+LRN_on = trigger_idxs_on(LRN_idx);
+LRN_off = trigger_idxs_off(LRN_idx);
+SMP_on = trigger_idxs_on(LRN_idx+3);
+SMP_off = trigger_idxs_off(LRN_idx+3);
+RST_on = trigger_idxs_on(LRN_idx+6);
+RST_off = trigger_idxs_off(LRN_idx+6);
+
+%% Concatenate trial data from EEG and EMG before computing coherence, 
+raw_d = EEG.data(5,:); %specify the data that you are going to work with (5 is c3 in our example data)
+for i=1:length(SMP_on)
+  temp_d2 = raw_d(SMP_on(i):SMP_off(i));
+  if i == 1
+    temp_d = temp_d2;
+  else
+    temp_d = horzcat(temp_d,temp_d2);
+end
+temp_d = temp_d./length(SMP_on);
+C3_SMP_vec = temp_d;
+
+for i=1:length(LRN_on)
+  temp_d2 = raw_d(LRN_on(i):LRN_off(i));
+  if i == 1
+    temp_d = temp_d2;
+  else
+    temp_d = horzcat(temp_d,temp_d2);
+end
+temp_d = temp_d./length(LRN_on);
+C3_LRN_vec = temp_d;
+
+for i=1:length(RST_on)
+  temp_d2 = raw_d(RST_on(i):RST_off(i));
+  if i == 1
+    temp_d = temp_d2;
+  else
+    temp_d = horzcat(temp_d,temp_d2);
+end
+temp_d = temp_d./length(RST_on);
+C3_RST_vec = temp_d;
+
+%same thing, but for EMG data
+raw_d = [] %specify the EMG data here
+for i=1:length(SMP_on)
+  temp_d2 = raw_d(SMP_on(i):SMP_off(i));
+  if i == 1
+    temp_d = temp_d2;
+  else
+    temp_d = horzcat(temp_d,temp_d2);
+end
+temp_d = temp_d./length(SMP_on);
+EMG1_SMP_vec = temp_d;
+
+for i=1:length(LRN_on)
+  temp_d2 = raw_d(LRN_on(i):LRN_off(i));
+  if i == 1
+    temp_d = temp_d2;
+  else
+    temp_d = horzcat(temp_d,temp_d2);
+end
+temp_d = temp_d./length(LRN_on);
+EMG1_LRN_vec = temp_d;
+
+for i=1:length(RST_on)
+  temp_d2 = raw_d(RST_on(i):RST_off(i));
+  if i == 1
+    temp_d = temp_d2;
+  else
+    temp_d = horzcat(temp_d,temp_d2);
+end
+temp_d = temp_d./length(RST_on);
+EMG1_RST_vec = temp_d;
+
+%% OLDER STUFF -- Compute average per trial (for testing only)
+%create average for SMP trials
+temp_d = zeros(1,8762); %THIS IS JUST FOR TESTING!
+for i=1:length(SMP_on)
+  temp_d2 = EEG.data(5,SMP_on(i):SMP_off(i));
+  disp(length(temp_d2))
+  temp_d = temp_d + temp_d2(1:8762);
+end
+temp_d = temp_d./length(SMP_on);
+SMP_avg = temp_d;
+
+%now for LRN trials
+temp_d = zeros(1,8762);%THIS IS JUST FOR TESTING
+for i=1:length(LRN_on)
+  temp_d2 = EEG.data(5,LRN_on(i):LRN_off(i));
+  disp(length(temp_d2))
+  temp_d = temp_d + temp_d2(1:8762);
+end
+temp_d = temp_d./length(LRN_on);
+LRN_avg = temp_d;
+
+
+%% EMG
+%initial code for EMG extraction, we assume the timecodes are the same
+% BUT we should likely check
+emg_1 = ALLEEG(2).data(1,:);
+
+%create emg average for SMP trials
+temp_d = zeros(1,8762); %THIS IS JUST FOR TESTING!
+for i=1:length(SMP_on)
+  temp_d2 = emg_1(SMP_on(i):SMP_off(i));
+  disp(length(temp_d2))
+  temp_d = temp_d + temp_d2(1:8762);
+end
+temp_d = temp_d./length(SMP_on);
+emg_1_SMP_avg = temp_d;
+
+
+emg_2 = ALLEEG(2).data(2,:);
+%create emg average for SMP trials
+temp_d = zeros(1,8762); %THIS IS JUST FOR TESTING!
+for i=1:length(SMP_on)
+  temp_d2 = emg_2(SMP_on(i):SMP_off(i));
+  disp(length(temp_d2))
+  temp_d = temp_d + temp_d2(1:8762);
+end
+temp_d = temp_d./length(SMP_on);
+emg_2_SMP_avg = temp_d;
+%%

spft/scripts/compute_maxPower.py

@@ -0,0 +1,29 @@
+'''
+code to calculate the max power of a signal, from here (https://stackoverflow.com/questions/58193209/function-to-determine-the-frequency-of-a-sinusoid)
+this relates to the SPFT sequences, which have max power at 0.52 Hz
+'''
+
+import numpy as np
+from matplotlib import pyplot as plt
+
+
+d = #your data vector [...]
+fs = 1000/12 #sampling frequency (of your signal, of course)
+y_fft = np.fft.fft(d-np.mean(d))           # Original FFT, offset removed
+y_fft = y_fft[:round(len(d)/2)] # First half ( pos freqs )
+y_fft = np.abs(y_fft)           # Absolute value of magnitudes
+y_fft = y_fft/max(y_fft)        # Normalized so max = 1
+
+freq_x_axis = np.linspace(0, fs/2, len(y_fft))
+plt.figure()
+plt.plot(freq_x_axis, y_fft, "o-")
+plt.title("Frequency magnitudes")
+plt.xlabel("Frequency")
+plt.ylabel("Magnitude")
+plt.grid()
+plt.tight_layout()
+plt.show()
+
+f_loc = np.argmax(y_fft) # Finds the index of the max
+f_val = freq_x_axis[f_loc] # The strongest frequency value
+print(f"The strongest frequency is f = {f_val}")