NeuroElf wiki

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== GSR data analysis ======

===== Motivation =====
For many experimental tasks (e.g. subjects follwing a set of specific but differing instructions), it is often helpful to collect secondary physiological measurements, such as galvanic skin-response (GSR) over time (both as a manipulation check but also as potential candidates for mediation analysis). The recorded data often requires manual clean-up and inspection, and this page should give you some pointers as to how to conduct some of these steps.

===== Overview =====
The following steps are covered:

  * reading in data files
  * re-sampling data
  * (pre-) filtering data
  * onset information extraction/synchronization from/with GSR data
  * obtaining peak-to-peak measures

===== Reading in data files =====
Currently, the following input formats are supported:
  * ACQ (up until version <= 3.9.7)
  * generic MAT file
  * TXT (use ''object = xff('*.ntt');'' or object = ''xff(filename, 'ntt');'' to read!)

To load the data from an ACQ file, use the following syntax:

<code matlab GSR_data_readacq.m>% use file selector of xff class
gsr = xff('*.acq', 'Please select the GSR-ACQ file to work with...');</code>

To use a variable in a MAT file, simply use the ''load MATFILE'' call in Matlab:
<code matlab GSR_data_loadmat.m>% load a mat file (e.g. an ACQ->MAT converted file)
load HPS1344_session1_GSR.mat;
 
% create new NTT (used for methods on data!)
ntt = xff('new:ntt');
 
% store data from mat file in ntt
ntt.Data = data;</code>

Or if you require to load the data from a text-based (e.g. log) file, you can use the following code:

<code matlab GSR_data_readtxt.m>% use file selector of xff class (you might have to change the file type
% in the selector from ''*.txt'' to ''*.*'' in case the extension is not TXT!
gsr = xff('*.ntt', 'Please select the text-based GSR file to work will...');

% the data is then available as
data = gsr.Data(:, CHANNEL);</code>

===== Re-sampling data =====
Very often, the raw data acquisition sampling rate is much higher than the required rate to perform any meaningful analyses on GSR data (e.g. peak-to-peak measurements). Here is an example of ten seconds of raw data acquired at 500Hz (5,000 data points):

<code matlab GSR_rawplot.m>% getting the data
data = gsr.Channel(1).AmplitudeScale .* gsr.RawData(1, :)';

% plotting from second 504 through 514
plot(data(252000:256999));</code>

{{:gsr_example_5secsraw.png|10 seconds of raw GSR data recording}}

As you can see, the very small signal variations are normally quite uninteresting (as the skin response only very sluggishly follows the stimulation/emotional state of the subject). Hence, it is usually a good idea to resample the signal to a more suitable resolution. To keep some of the finer details (and being able to determine latency values with a reasonable resolution), for this example I resampled the data to 100Hz using the [[acq.Resample|ACQ::Resample]] or [[ntt.Resample|NTT::Resample]] method (internally using the [[resampleaa]] function which works a little better than using just ''downsample''):

<code matlab GSR_resampledplot.m>% resampling the data
% the first channel is actual physiological data and needs hi-quality
% resampling (cubic/gaussian average), the stim channels 2:9
% only need indexing into the correct temporal samples
gsr.Resample(100, struct('cubic', [true, false(1, 8)]));

% plotting same time index with stronger line
data = gsr.ChannelData(1);
set(plot(data(50400:51399)), 'LineWidth', 2);</code>

{{:gsr_example_5secsresampled.png|100Hz-resampled GSR data}}

===== (Pre-) Filtering data =====
While a lot of the very high-frequency noise is (naturally) gone after the resampling, there is still a lot of signal variation that is (very likely) unrelated to any task/stimulation the subject could have reacted to. Further the signal pre-processing, I then applied a difference (first discreet derivative, ''diff'')-filter that removes artifacts in signals which are unlikely to be in the desired spectrum using the [[acq.Filter|ACQ::Filter]]/[[ntt.Filter|NTT::Filter]] method (which makes use of the [[prefilter]] function):

<code matlab GSR_prefilter.m>% call obj.Filter with default parameters for first channel
gsr.Filter(1);

% plotting same time index with stronger line
data = gsr.ChannelData(1);
set(plot(data(50400:51399)), 'LineWidth', 2);</code>

{{:gsr_example_5secsfiltered.png|filtered GSR data}}

Now the data looks reasonably clean to undergo further processing (e.g. peak-to-peak analysis).

===== Onset information =====
In our first example, the onsets are available are channels in the ACQ object we use. The condition (type of stimulus) is encoded in three channels (2 through 4) and the within-trial stimulus distinction (fixation, instruction, pre-stimulus, stimulus, rating) is given by channels 5 through 9. These channels either carry a zero value (condition/stimulus not present) or a value greater than 0 (present).

To obtain a list of onsets, we first find the onsets of all conditions and then lookup those values in the vectors required for the different conditions (in our case three), whereas the onset of interest is coded as channel 8 (stimulus):

<code matlab GSR_acq_stimonsets.m>% all stimuli onsets (use 1 + to cope for the fact that the
% first value in diff rather codes the second sample in the data!
all_onsets = 1 + find(diff(gsr.ChannelData(8)) > 0);

% now split into three conditions
neu = 0.01 * all_onsets(gsr.ChannelData(2, all_onsets) > 0);
neg = 0.01 * all_onsets(gsr.ChannelData(3, all_onsets) > 0);
rea = 0.01 * all_onsets(gsr.ChannelData(4, all_onsets) > 0);</code>

As window length for the inspection of the GSR data, we choose 15 seconds, so in the resolution of the resampled data that means 1,500 samples.

In our second example, the data only contains the onsets (not the condition information). This has to be established separately (which can be even done **after** extracting the peak-to-peak values, thus allowing for a blind-analysis of data, as the examiner doesn't know which onset belongs to which specific condition):

<code matlab GSR_plain_onsets.m>% loading the data
% (in this case a MAT file written by Acqknowledge for ACQ Version >= 4)
% which is in 200Hz sampling frequency with
% GSR = channel 1 and
% ECG = channel 2
load 6701.acq.mat

% the stimulus channel in this example is channel 3,
% coded with 0 Volts for no stimulus and 5V for stimulus
% thus, to find the time points (samples) on which a
% stimulus occurs, we look for spots in this channel
% where the value changes from below 2.5V to above 2.5V
stimactive = (data(:, 3) > 2.5);
onsets = find(~stimactive(2:end) & stimactive(1:end-1));

% creating the necessary information for plotcurves
% the onset is shifted by 200 samples (1 second) to show
% a window of [-1 .. 14] (in seconds)
c = cell(numel(onsets), 2);
for cc = 1:numel(onsets)
    c{cc,1} = sprintf('Onset %d, t=%.3f', cc, onsets(cc)/200);
    c{cc, 2} = [1, onsets(cc)-200, onsets(cc)+2800];
end

% calling plotcurves
plotcurves(data, struct('curves', {c}, 'freq', 200));</code>

===== Peak-to-peak analysis =====
Once we are happy with the onsets, the [[plotcurves]] function can be used to show each window from onset to offset and then determine the peak-to-peak response:

<code GSR_peaktopeak.m>% create required options argument for plotcurves
gsrchannel = 1;
neucurves = cell(numel(neu), 2);
negcurves = cell(numel(neg), 2);
reacurves = cell(numel(rea), 2);
for sc = 1:numel(neu)
    neucurves{sc, 1} = sprintf('LookNeu %02d at t=%.3fs', sc, neu(sc));
    neucurves{sc, 2} = [gsrchannel, neu(sc) - 2, neu(sc) + 14];
end
for sc = 1:numel(neg)
    negcurves{sc, 1} = sprintf('LookNeg %02d at t=%.3fs', sc, neg(sc));
    negcurves{sc, 2} = [gsrchannel, neg(sc) - 2, neg(sc) + 14];
end
for sc = 1:numel(rea)
    reacurves{sc, 1} = sprintf('ReaNeg %02d at t=%.3fs', sc, rea(sc));
    reacurves{sc, 2} = [gsrchannel, rea(sc) - 2, rea(sc) + 14];
end

% spots
spot = {[neu; neg; rea]};
spotnames = {'Stim onset'};

% combine curves
curves = cat(1, neucurves, negcurves, reacurves);

% define sets
sets = { ...
    'LookNeu', (1:numel(neu))'; ...
    'LookNeg', (numel(neu)+1:numel(neu)+numel(neg))'; ...
    'ReaNeg',  (numel(neu)+numel(neg)+1:size(curves,1))'};

% use plot curves
plotcurves(gsr, struct( ...
    'curves',    {curves}, ...
    'sets',      {sets}, ...
    'spot',      {spot}, ...
    'spotnames', {spotnames}));</code>