% 1サンプルずつfor文で角度を計算しようとした例
for i = 1:length(hipX)
    thigh = [hipX(i)-kneeX(i), hipY(i)-kneeY(i)];
    shank = [kneeX(i)-ankleX(i), kneeY(i)-ankleY(i)];
    kneeAngle(i) = acosd(dot(thigh,shank)/(norm(thigh)*norm(shank)));
end
% → 動くが遅く、データの向き（行/列）がずれるとエラーになる

% すべてのサンプルをまとめて列ベクトルとして読み込む
thigh = [hipX - kneeX, hipY - kneeY];      % Nx2行列
shank = [kneeX - ankleX, kneeY - ankleY];  % Nx2行列

% 各行（=各サンプル）ごとにドット積とノルムをまとめて計算
dotProd = sum(thigh .* shank, 2);
normProd = vecnorm(thigh,2,2) .* vecnorm(shank,2,2);
kneeAngle = acosd(dotProd ./ normProd);    % for文なしで一気に計算

% ノイズ除去（4次バターワースローパスフィルタ、遮断周波数6Hz）
fs = 100; % サンプリング周波数[Hz]
[b,a] = butter(4, 6/(fs/2), 'low');
kneeAngleSmooth = filtfilt(b,a,kneeAngle);

plot(kneeAngleSmooth); xlabel('サンプル'); ylabel('膝関節角度 [deg]');

clc; clear; close all;   % コマンド履歴・変数・図をすべてクリアしてから開始
T = readtable('data.csv');  % データ読み込み
disp(T.Properties.VariableNames); % 列名を確認してから解析に進む

function angle = calcJointAngle(proximal, distal, joint)
% proximal, distal, joint はそれぞれ [x y] または [x y z] のNx2/Nx3行列
    v1 = proximal - joint;
    v2 = distal - joint;
    cosTheta = sum(v1.*v2,2) ./ (vecnorm(v1,2,2).*vecnorm(v2,2,2));
    angle = acosd(cosTheta);
end

function dataSmooth = smoothSignal(dataRaw, fs, cutoff)
    [b,a] = butter(4, cutoff/(fs/2), 'low');
    dataSmooth = filtfilt(b, a, dataRaw);
end
% 使用例： emgSmooth = smoothSignal(emgRaw, 1000, 20);

function compareBeforeAfter(pre, post, unitLabel)
    [~, p, ~, stats] = ttest(pre, post);
    figure;
    boxplot([pre(:) post(:)], 'Labels', {'訓練前','訓練後'});
    ylabel(unitLabel);
    title(sprintf('訓練前後の比較 (p = %.3f, t = %.2f)', p, stats.tstat));
end
% 使用例： compareBeforeAfter(preROM, postROM, '可動域 [deg]');

% 個人情報を含んだファイル名や変数名をそのまま使ってしまう例
patientData_Yamada_Taro_19450101 = readtable('yamada.csv');
% → ファイル名・変数名・保存フォルダに個人が特定できる情報が残ってしまう

% 事前に匿名IDへの対応表を、解析用フォルダとは別の場所で厳重に管理する
patientData_ID0231 = readtable('ID0231.csv');
% 対応表（氏名⇔ID）は、院内の個人情報管理ルールに沿った場所にのみ保管する