CSN_Skill/MATLAB Scripts/free_move_correlation.m

%% Freely Moving Kinetics 
clear;
%% Import file and extract full signal

[files, path] = uigetfile({'*.MotoTrak;*.ArdyMotor'},...
    'Select Mototrak File', ...
    'multiselect','on');
try 
    cd(path);
catch err
    return
end
if ischar(files)                                                      %If only one file was selected...
    files = {files};                                                    %Convert the string to a cell array.
end
for i = 1:length(files)
    D = dir(files{1,i});
    Bytes(i) = D.bytes;
end
[row, col] = find(Bytes>500);                             
files = files(col);
knob_data.trial_length = 500;
knob_data.num_sessions = length(files);
knob_data.combined_sessions_length = 0;
knob_data.session_length = nan(1, length(files));
knob_data.max_session_trials = 0;
for s = 1:knob_data.num_sessions
    try                                                                     %Try to read in the data file...
        [~,filename,ext] = fileparts(files{s});
        switch ext
            case '.ArdyMotor'
                data = ArdyMotorFileRead(files{s});
            case '.MotoTrak'
                data = MotoTrakFileRead(files{s});
                data = MotoTrak_to_ArdyMotor(data);
                data.rat = data.subject;
        end
        %         temp = ArdyMotorFileRead(files{f});                                 %Read in the data from each file.
    catch err                                                               %If an error occurs...
        warning(['ERROR READING: ' files{s}]);                              %Show which file had a read problem...
        warning(err.message);                                               %Show the actual error message.
        continue
    end
%     data = ArdyMotorFileRead(files{s});
    temp = length(data.trial);
    knob_data.combined_length = temp+knob_data.combined_sessions_length;
    knob_data.session_length(s) = temp;
    
    if knob_data.session_length(s) > knob_data.max_session_trials
        knob_data.max_session_trials = knob_data.session_length(s);
    end
end
knob_data.combined_trials = nan(knob_data.combined_length, 500);
knob_data.trial = nan(knob_data.max_session_trials, 500, knob_data.num_sessions);

% for t = 1:length(data.trial)                                                %Step through each trial.
%         data.trial(t).timestamps = ...
%             double(data.trial(t).sample_times)/86400000 + ...
%             data.trial(t).starttime;                                        %Convert the sample times to serial date numbers.
% end
% 
% time = vertcat(data.trial.timestamps);
for i = 1:length(files)
    signal(:,i) = vertcat(data.trial.signal);
end

%% Find and classify peaks in the force signal

lower = 15;
upper = 18;
three = 3;
nine = 9;

[pks, locs, widths] = findpeaks(signal, 'MinPeakHeight', three, 'MinPeakDistance', 8);  %Finds peaks and their corresponding locations for pulls over 3 grams at least 0.4s apart

peakTimes = locs;
pks = pks(1: length(peakTimes));
NumPeaks = length(pks); 
trial_start = 1:1:NumPeaks;
trial_end = 1:1:NumPeaks;

width_full = widths;

avg_pull = mean(width_full);

for i = 1:NumPeaks
    trial_start(1,i) = peakTimes(i,1) - avg_pull;
    trial_end(1,i) = peakTimes(i,1) + avg_pull;
end

R_trial_start = floor(trial_start);
R_trial_end = ceil(trial_end);

pullTimes = reshape([R_trial_start(:) R_trial_end(:)]', [1, 2* NumPeaks]); 

behavior = cell(1, NumPeaks);

for j = 1:2:length(pullTimes)
    behavior{1,j} = signal(pullTimes(j) : pullTimes(j+1), 1);
end

behavior(cellfun('isempty',behavior)) = [];

[pks, locs, widths] = findpeaks(signal, 'MinPeakHeight', lower, 'MinPeakDistance', 8);

if pks > 0
    for i = 1:length(pks)
        if pks(i,:) > upper
           locs(i,:) = 1;
           widths(i,:) = 1;
        end
    end

    locs(locs == 1) = [];
    widths(widths == 1) = [];

    successPeaks = locs;
end

[~, index_successPeaks, index_peakTimes] = intersect(successPeaks,peakTimes,'rows');
successloc = zeros(1, length(peakTimes)); 

successloc(:,index_peakTimes) = 1;
successloc(1,1) = 0;
  
for j = 1:size(behavior,2) - 1
    temp = corrcoef(behavior{1,j}, behavior{1,j+1});
    pull_corr{1,j} = temp(1,2);
end

avg_pull_corr = mean(abs(cell2mat(pull_corr)));

for i = 1:size(behavior,2) - 1
    if successloc(:,i) == 1
        temp = corrcoef(behavior{1,i}, behavior{1,i-1});
        suc_prevcorr{1,i} = temp(1,2);
    else
        temp = NaN;
        suc_prevcorr{1,i} = NaN;
    end
end

for i = 1:size(behavior,2) - 1
    if successloc(:,i) == 1
        temp = corrcoef(behavior{1,i}, behavior{1,i+1});
        suc_postcorr{1,i} = temp(1,2);
    else
        temp = NaN;
        suc_postcorr{1,i} = NaN;
    end
end

avg_presuccess_corr = nanmean(abs(cell2mat(suc_prevcorr)));
avg_postsuccess_corr = nanmean(abs(cell2mat(suc_postcorr)));

% %%
% numNeurons = size(Cd_N,1);
% activity = cell(numNeurons, length(pullTimes) / 2);
% 
% for i = 1:numNeurons
%     for j = 1:2:length(pullTimes)
%         activity{i,j} = Cd_N(i, pullTimes(j) : pullTimes(j+1));
%     end
% end
% 
% activity(cellfun('isempty',activity)) = [];
% activity = reshape(activity, [numNeurons, length(meanFrames)/2]);
% activity2 = cell2mat(activity);