Simplification of of bunch of in-lab scripts to apply Time Shift PCA on MEG KIT data

test_tspca.m

% Author: Jean-Remi King <jeanremi.king@gmail.com>
%
% Licence : To-be-determined

addpath('/media/DATA/Pro/Toolbox/NoiseTools/')  % http://audition.ens.fr/adc/NoiseTools/
demean = @(x) x - repmat(mean(x, 1), [size(x, 1), 1]);

%% Read data, prepare info
raw = load('raw.mat');
data_empty = transpose(raw.empty(1:160, :)); % empty room: chan x time
data_subject = transpose(raw.raw(1:160, :)); % subject data: chan x time
data_subject_clean = data_subject; % initialize clean data
[n_samples, n_chans] = size(data_subject);

% parameters specific to our KIT system
chans_sensors = 1:157; % (KIT: axial gradiometers: sensors over the head)
chans_ref = 158:160;  % (KIT: ref magnetometer)
chans_bad = [20, 40, 64, 112, 115, 152] + 1; % python & matlab indices starts at 0 and 1 respectively
chans_good = setdiff((chans_sensors), chans_bad);


%% 1. Cleans with Least Square -----------------------------------------------
% Find Least Square projection of the reference channels in an empty room and
% remove the corresponding component from the recordings of a subject.
%
% adapted from Adeen Flinker 6/2013 (adeen.f@gmail.com) LSdenoise.m
empty_sensors = demean(data_empty(:, chans_sensors));
empty_sensors(:, chans_bad) = 0;  % remove bad channels
empty_ref = demean(data_empty(:, chans_ref));
W = pinv(empty_ref) * empty_sensors;  % learn the least squares weights
clear('empty_sensors', 'empty_ref', 'data_empty')

subject_ref = demean(data_subject(:, chans_ref));
subject_sensors = demean(data_subject(:, chans_sensors));
subject_sensors = subject_sensors - subject_ref * W;
subject_sensors(:, chans_bad) = 0;

new_ref = subject_sensors * pinv(W);
data_subject(:, chans_sensors) = subject_sensors;
data_subject(:, chans_ref) = new_ref;
clear('subject_ref', 'empty_sensors')

%% 2. Cleans with TSPCA -----------------------------------------------
% Adapted from Dan Hertz dbhertz@umd.edu sqdDenoise.m
% To avoid memory overload, we have to compute TSPCA in blocks
block_size = 20000;
shifts = -100:100;
block_start = max(0, -min(shifts));
block_stop = max(0, max(shifts));
block_increment = block_size - block_start - block_stop;
n_blocks = floor((n_samples - 1) / block_increment);
this_block_start = 0;
this_block_stop = 0;

% MAIN
for this_block = 1:n_blocks
    % select block
    disp(this_block / n_blocks)
    this_block_start = (this_block - 1) * block_increment + 1;  % Start at entry 1 because there is no sample 0
    this_block_stop = (this_block - 1) * block_increment + block_size;
    new_index = 1 + block_start:block_size-block_stop;
    insertLow = this_block_start + block_start; % where are we putting this into the sqd file array
    insertHigh = this_block_stop - block_stop;
    data_block = data_subject(this_block_start:this_block_stop, :);
    data_block_clean = zeros(size(new_index, 2), 157);  % init clean data

    % demean the data channels here (but not ref or trigger)
    % Devigne's algorithm assumes 0 mean => demean per block
    data_block_sensors = demean(data_block(:, chans_good));
    data_block_ref = data_block(:, chans_ref);

    % Calculate a weight matrix for tspca based on the channels
    % that may be clipping
    % XXX JRK: I don't understand this part
    weightsForTSPCA = deweightClippedChannels(data_block_sensors, chans_good);
    weightsByTimeTSPCA = min(weightsForTSPCA, [], 2);
    weightsForSNS = weightsByTimeTSPCA(new_index, :);
    weightsToRestore = weightsForTSPCA(new_index, :);
    weightsAreZero = weightsToRestore==0;

    %% MAIN FUNCTION from Alain de Cheveigné
    data_block_clean_sensors = tsr_dan(data_block_sensors, data_block_ref, shifts, weightsByTimeTSPCA);

    %% JRK: I don't unstadard this part, and my lab actually doesn't uses it, but just in case
    near_neighbors = 0;
    if near_neighbors > 1
        data_block_clean_sensors = sns_weighted(data_block_clean_sensors, near_neighbors, 0, weightsForSNS);
    end

    % save
    data_block_clean(:, chans_good) = data_block_clean_sensors;
    data_subject_clean(insertLow:insertHigh, chans_sensors) = data_block_clean;
    data_subject_clean(insertLow:insertHigh, chans_ref) = data_block_ref(new_index, : );

end

save('raw_clean.mat', 'data_subject_clean')