Standard deviation for RMS comparison of two windows

Just looking for guidance on my math to calculate some standard deviation values. I am trying to compare some EMG motor unit data by calculating the RMS in two time windows around motor unit spikes. I have the code where it calculates the RMS value in the two units and comapres the data for which window has more dominance. In my graph I would like to add standard deviation bars to the points to visualize the data accurately. I am pretty sure I am not doing my math right so any guidance would be appreciated. Below is a picture of one of the motor units that I am using to calculate the two windows out of the total window of 0 to 50ms.

% Add a single indicator for each motor unit based on the RMS percentage
for j = 1:numMuscles
    if rmsPercentage(j) < 40
        text(j, rmsPercentage(j) + 2, 'M Wave', 'HorizontalAlignment', 'center', 'Color', 'b');
    elseif rmsPercentage(j) > 60
        text(j, rmsPercentage(j) + 2, 'H Wave', 'HorizontalAlignment', 'center', 'Color', 'r');
    else
        text(j, rmsPercentage(j) + 2, 'M/H Wave', 'HorizontalAlignment', 'center', 'Color', 'g'); % Use a different color for combination
    end
end

% Define the time windows
window1 = [0, 30];  % Window 1: 0 to 30 ms
window2 = [31, 50];  % Window 2: 31 to 50 ms

% Initialize arrays to store the RMS activity, percentages, and variability measures
rmsActivityWindow1 = zeros(numMuscles, 1);
rmsActivityWindow2 = zeros(numMuscles, 1);
rmsPercentage = zeros(numMuscles, 1); % For percentage calculation
stdRMSPercentage = zeros(numMuscles, 1); % Standard deviation expressed as a percentage

% Calculate the corresponding sample indices for each window
startIdx1 = round(window1(1) * samplingFrequency / 1000) + 1;
endIdx1 = round(window1(2) * samplingFrequency / 1000);
startIdx2 = round(window2(1) * samplingFrequency / 1000) + 1;
endIdx2 = round(window2(2) * samplingFrequency / 1000);

% Ensure that indices do not exceed array bounds
if endIdx1 > size(averagePulseTrainSignals, 2)
    endIdx1 = size(averagePulseTrainSignals, 2);
end
if endIdx2 > size(averagePulseTrainSignals, 2)
    endIdx2 = size(averagePulseTrainSignals, 2);
end

% Loop through each motor unit to calculate RMS activity and standard deviation
for j = 1:numMuscles
    % Calculate RMS for Window 1
    rmsActivity1 = sqrt(mean(averagePulseTrainSignals(j, startIdx1:endIdx1).^2, 'omitnan'));
    rmsActivityWindow1(j) = rmsActivity1;

    % Calculate RMS for Window 2
    rmsActivity2 = sqrt(mean(averagePulseTrainSignals(j, startIdx2:endIdx2).^2, 'omitnan'));
    rmsActivityWindow2(j) = rmsActivity2;

    % Calculate the mean for both windows
    meanWindow1 = mean(averagePulseTrainSignals(j, startIdx1:endIdx1), 'omitnan');
    meanWindow2 = mean(averagePulseTrainSignals(j, startIdx2:endIdx2), 'omitnan');

    % Calculate the percentage for dominance
    if rmsActivity1 + rmsActivity2 ~= 0
        rmsPercentage(j) = (rmsActivity2 / (rmsActivity1 + rmsActivity2)) * 100; % Calculate percentage

        % Calculate standard deviation for each window
        std1 = sqrt(rmsActivity1^2 - meanWindow1^2);
        std2 = sqrt(rmsActivity2^2 - meanWindow2^2);
        
        % Combine standard deviations as a percentage of RMS, then divide by 10
        combinedStd = (std1 + std2) / 2; % Average of standard deviations
        stdRMSPercentage(j) = (combinedStd / mean([rmsActivity1, rmsActivity2])) * 100 / 10; % Express as percentage and divide by 10
    else
        rmsPercentage(j) = NaN; % Avoid division by zero
        stdRMSPercentage(j) = NaN; % Avoid division by zero
    end
end

% Create a figure for the RMS percentage plot
figure;
xPositions = 1:numMuscles; % X positions for motor units
scatter(xPositions, rmsPercentage, 'filled', 'MarkerFaceColor', 'b');

% Add error bars for the standard deviation percentage
hold on;
errorbar(xPositions, rmsPercentage, stdRMSPercentage, 'k', 'LineStyle', 'none');

% Add a single indicator for each motor unit based on the RMS percentage
for j = 1:numMuscles
    if rmsPercentage(j) < 40
        text(j, rmsPercentage(j) + 2, 'M Wave', 'HorizontalAlignment', 'center', 'Color', 'b', 'FontSize', 10);
    elseif rmsPercentage(j) > 60
        text(j, rmsPercentage(j) + 2, 'H Wave', 'HorizontalAlignment', 'center', 'Color', 'r', 'FontSize', 10);
    else
        text(j, rmsPercentage(j) + 2, 'M/H Wave', 'HorizontalAlignment', 'center', 'Color', 'g', 'FontSize', 10);
    end
end

% Set labels and title with increased font size
xlabel('Motor Units', 'FontSize', 14);
ylabel('RMS Percentage', 'FontSize', 14);
title('RMS Activity of Motor Units as Percentage with Std Dev Error Bars', 'FontSize', 16);
ylim([0 100]); % Set y-limits from 0 to 100 percent
xticks(1:numMuscles); % Set x-ticks to correspond to motor units
xticklabels(1:numMuscles); % Label the x-ticks with numbers
grid on; % Add a grid for clarity

% Adjust font size for axes
set(gca, 'FontSize', 12);

hold off;

Answer 1

[Star Strider]

Root-mean-square (RMS) calculations are best suited to estimate the D-C amplitude of a signal over long periods or the entire signal. To look at individual MUAPS or other pulses in an ensemble, I characteristically calculate the mean, and then 95% confidence intervals using the standard error of the mean.

t = linspace(0, 1, 101);
y =  sinc(-5:0.1:5);
y = y + randn(50,numel(y))*0.25;

figure
plot(t, y)
grid
xlabel('t')
ylabel('V')

ymean = mean(y)

ymean = 1×101

0.0184 0.0069 0.0878 0.0740 0.0773 0.0240 0.1220 0.0242 0.0791 -0.0267 -0.0022 -0.0699 -0.0063 -0.0649 -0.1165 -0.1017 -0.0698 -0.1258 -0.0553 -0.0812 0.0644 -0.0136 0.0946 0.1068 0.0760 0.1399 0.1491 0.1499 0.0892 0.0781

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

ySEM  = std(y)/sqrt(size(y,1))

ySEM = 1×101

0.0313 0.0366 0.0329 0.0312 0.0368 0.0359 0.0338 0.0352 0.0301 0.0323 0.0357 0.0378 0.0364 0.0362 0.0293 0.0354 0.0393 0.0339 0.0358 0.0369 0.0317 0.0381 0.0362 0.0330 0.0352 0.0326 0.0362 0.0437 0.0369 0.0358

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

cv = tinv([0.025 0.975], size(y,1)-1)

cv = 1×2

-2.0096 2.0096

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

ci95 = ySEM .* cv(:)

ci95 = 2×101

-0.0629 -0.0736 -0.0662 -0.0627 -0.0740 -0.0721 -0.0680 -0.0708 -0.0606 -0.0649 -0.0717 -0.0760 -0.0731 -0.0728 -0.0589 -0.0712 -0.0790 -0.0682 -0.0719 -0.0742 -0.0638 -0.0766 -0.0728 -0.0663 -0.0707 -0.0656 -0.0727 -0.0878 -0.0742 -0.0718 0.0629 0.0736 0.0662 0.0627 0.0740 0.0721 0.0680 0.0708 0.0606 0.0649 0.0717 0.0760 0.0731 0.0728 0.0589 0.0712 0.0790 0.0682 0.0719 0.0742 0.0638 0.0766 0.0728 0.0663 0.0707 0.0656 0.0727 0.0878 0.0742 0.0718

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure
hp1 = plot(t, ymean, '-', 'DisplayName','\mu');
hold on
hp2 = patch([t flip(t)], [ci95(1,:)+ymean flip(ci95(2,:)+ymean)],  'r', 'FaceAlpha',0.5, 'DisplayName','±95% CI', 'EdgeColor','none');
hold off
legend('Location','best')
grid

EDIT — Corrected typographical errorrs.

.

[Austin]

Thank you for the response and guidance, I will implement this into the script!

[Star Strider]

My pleasure!