calcNeuralActivations.m

% This function is part of the NMSM Pipeline, see file for full license.
%
% This function calculates the neural activations given the
% muscleExcitation signals using backward finite difference approximation
%
% (3D mat of num, 3D mat of num, 2D array of num) -> (3D mat of num)
% returns the neural activations


function neuralActivations = calcNeuralActivations(muscleExcitation, ...
    activationTimeConstants, emgTime, numPaddingFrames)
activationTimeConstants = activationTimeConstants / 100;
deactivationTimeConstants = 4 * activationTimeConstants;
% equation 6 from Meyer 2017
inverseDeactivationTimeConstants = 1 ./ deactivationTimeConstants;
% equation 5 from Meyer 2017
differenceOfTimeConstants = (1 ./ activationTimeConstants) - ...
    inverseDeactivationTimeConstants;
% Each emg trial could have different time interval
trialSpecificTimeInterval = mean(diff(emgTime, 1, 2), 2);
neuralActivations = zeros(size(muscleExcitation));
% equation 7 from Meyer 2017
for j = 3:size(muscleExcitation, 3)
    finiteDifferenceTimeConstants = 2 * trialSpecificTimeInterval .* ...
        (differenceOfTimeConstants .* muscleExcitation(:, :, j) + ...
        inverseDeactivationTimeConstants);
    finiteDifferenceDenominator = finiteDifferenceTimeConstants + 3;
    finiteDifferenceNumerator = finiteDifferenceTimeConstants .* ...
        muscleExcitation(:, :, j) + 4 * neuralActivations(:, :, j-1) - ...
        neuralActivations(:, :, j-2);
    neuralActivations(:, :, j) = finiteDifferenceNumerator ./ ...
        finiteDifferenceDenominator;
end
neuralActivations = neuralActivations(:, :, numPaddingFrames + ...
    1:size(emgTime, 2) - numPaddingFrames);
neuralActivations(neuralActivations < 0) = 0; %remove negative neural activ
end