3軸加速度センサーを用いたシンプルな歩数推定(階段の上り下り)

phyphox(Physical Phone Experiments)で階段を上り下りしたときの3軸加速度センサーのデータを記録し、データ記録中に歩いた歩数を推定したときのメモになります。

先日、こちらのページに記載した方法と同じ方法で歩数を推定しました。使用したJavaScriptも同じです。Google Sheetsにアップロードしたデータの参照先のURLのみ変更しました。階段を上り下りしたときの歩数も同じ方法で大きな誤差なく推定することができていました。

1. こちらの手順で歩行中の3軸加速度センサーのデータをphyphoxを使って記録します。

2. phyphoxで記録した3軸加速度のデータは約0.01秒ごとに記録されています。下図のデータは約43秒間のデータになります。記録したデータを指定した間隔で間引きすることもありますが、今回は間引きしませんでした。

3. 取得したデータの移動平均を計算し、細かな動きやノイズの影響を抑えます。下の図は約0.2秒間の移動平均を取った結果です。過去約0.2秒間のセンサーの値の平均値を各時刻の歩行者の重心の加速度の推定値とします。前回と同様、測定開始直後の約2秒間と測定終了直前の約3秒間のデータは取り除くことにしました。測定する際、測定開始ボタンを押してから4秒ほど経過した後で階段を上り始め、階段を下り終えて3秒以上経過してから測定停止ボタンを押しました。

4. 3軸加速度のデータから重力方向を推定するため、約2秒間の移動平均を取りました。過去約2秒間の3軸加速度センサーの値の平均値ベクトルの向きをその時刻のスマートフォンの向きに対する重力方向の推定値とします。

5. 先の3.で計算した細かな動きやノイズの影響を抑えた各時刻の3軸加速度のセンサー値を、4.で計算した重力方向に射影します。歩行時の重心の動き(加速度)の推定値の重力方向の成分を得ることができます。前回同様、測定を開始してから4秒後以降のデータを対象とすることにしました。

3軸加速度センサーで取得した加速度の3次元ベクトルを重力方向に射影した値は、下記のコードで計算しています。

        const ax = data_short_moving_average.getValue(rowIndex, 1);
        const ay = data_short_moving_average.getValue(rowIndex, 2);
        const az = data_short_moving_average.getValue(rowIndex, 3);

        const gx = data_long_moving_average.getValue(rowIndex, 1);
        const gy = data_long_moving_average.getValue(rowIndex, 2);
        const gz = data_long_moving_average.getValue(rowIndex, 3);
        const absolute_value_of_g = Math.sqrt(gx * gx + gy * gy + gz * gz);

        const accleration_projected_to_garavitational_direction = (ax * gx + ay * gy + az * gz) / absolute_value_of_g;

6. 細かな動きやノイズの影響を抑えるため、5.で計算した重心の動きの重力方向の推定値の約0.2秒間の移動平均を計算します。

7. 重心加速度の重力方向成分は、右足でも左足でも一歩ごとに歩く周期で振動します。重力加速度 9.80665[m/s^2] を中心として、大きくなったり小さくなったりします。先の6.の図にプロットした約36秒間のデータの平均を計算し、プロットしました。meanという名前の直線が平均値になります。比較のため重力加速度 9.80665[m/s^2] もプロットしました。データの平均値(mean)と重力加速度の値は少しずれていますが、重心加速度の重力方向成分はほぼ重力加速度を中心として歩行の周期で振動しています。

8. 歩数を推定する処理を追加してプロットしました。階段は、途中に2箇所の踊り場がある15段の階段です。2箇所の踊り場では上りも下りも2歩ずつ歩きました。17歩歩いて上り、Uターンしてから17歩歩いて下り、元の位置に戻ってきました。下りでは、踊り場での2歩目の上下方向の重心移動は小さく、カウントされていないようでした。上り終えた後、あまり大きく歩かずその場でUターンしたのですが、そのときの歩数は1歩としてカウントされていました。

歩数は下記のようなコードで計算しています。各時刻の重心加速度の重力方向成分と重力加速度の差を計算し、それを加算していきます。加算していった値があるプラスの閾値を超える時刻とマイナスの閾値を下回る時刻が、一歩ごとに一回ずつ見つかるように閾値を調節します。(歩行の周期より短い周期で重心加速度の重力方向成分が小さく揺れた場合はカウントアップせず、歩いたときにはカウントアップできるように閾値の大きさを調節します。)

    let steps = 0;
    let sum = 0;
    let is_positive = false;
    let is_negative = false;

    for (let rowIndex = 0; rowIndex < numRows; rowIndex++) {

        const diff = data.getValue(rowIndex, 1) - 9.80665;

        if ( (diff > 0 && is_positive == false) || (diff < 0 && is_negative == false) ) {
            sum += diff;
        }

        if (sum > positive_threshold) {
            steps++;
            is_positive = true;
            is_negative = false;
            sum = 0;
        }

        if (sum < negative_threshold) {
            is_positive = false;
            is_negative = true;
            sum = 0;
        }

        data.setCell(rowIndex, 4, steps);
    }

この例では記録済みのデータを対象としていますが、こちらに記載した歩数推定方法は計測したデータをリアルタイムで処理する歩数計にも適用することができます。


補足:この投稿には下記のJavaScriptを使用しました。

<script type="text/javascript" src="https://www.gstatic.com/charts/loader.js"></script>
<script type="text/javascript">

google.charts.load('current', {packages:['corechart']});
google.charts.setOnLoadCallback(Spreadsheet);

function Spreadsheet() {
    var query = new google.visualization.Query('https://docs.google.com/spreadsheets/d/1LPud6SRVzwOSsTKrg2OmTvQw3Khau8vz5ofZcRDmNfU/edit?usp=sharing');
    query.send(drawChart);
}

const skip_length = 1;
const time_interval = 0.01;

function drawChart(response) {
    const data = response.getDataTable();
    const numRows = data.getNumberOfRows();
    const numThinnedRows = Math.floor(numRows / skip_length);
    const remainder = numRows % skip_length;

    for (let rowIndex = 0; rowIndex < numThinnedRows; rowIndex++) {
        data.removeRows(rowIndex, skip_length - 1);
    }

    if (remainder != 0) {
        data.removeRows(numThinnedRows, remainder);
    }

    // data thinning
    const options = {title: 'phyphox 3-axis acceleration sensor data (after data thinning)',
                     hAxis: {title: 'time[s]'},
                     vAxis: {title: 'acceleration [m/s^2]'}};
    const chart = new google.visualization.LineChart(document.getElementById('after_data_thinning'));
    chart.draw(data, options);

    // short moving average
    let title = 'phyphox 3-axis acceleration sensor data (short moving average)';
    let term = 0.2; // seconds
    let moving_average_length = Math.floor(term / time_interval / skip_length);
    let opening_term = 2.0; // seconds
    let closing_term = 3.0; // seconds
    const data_short_moving_average = data.clone();
    drawMovingAverage(data_short_moving_average, moving_average_length, 'short_moving_average', title, opening_term, closing_term);

    // long moving average : estimate direction of gravity
    title = 'phyphox 3-axis acceleration sensor data (long moving average)';
    term = 2.0; // seconds
    moving_average_length = Math.floor(term / time_interval / skip_length);
    opening_term = 2.0; // seconds
    closing_term = 3.0; // seconds
    const data_long_moving_average = data.clone();
    drawMovingAverage(data_long_moving_average, moving_average_length, 'long_moving_average', title, opening_term, closing_term);

    // acceleration in the direction of gravity
    title = 'acceleration in the direction of gravity';
    opening_term = 2.0; // seconds
    drawAccelerationProjectedToGravitaionalDirection(data_short_moving_average, data_long_moving_average,
                                                     'acceleration_gravitaional_direction', title, opening_term);

    // acceleration in the direction of gravity (short moving average)
    title = 'acceleration in the direction of gravity (short moving average)';
    term = 0.2; // seconds
    moving_average_length = Math.floor(term / time_interval / skip_length);
    opening_term = 0.2; // seconds
    closing_term = 0; // seconds
    drawMovingAverage(data_short_moving_average, moving_average_length, 'acceleration_gravitaional_direction_moving_average',
                      title, opening_term, closing_term);

    // acceleration in the direction of gravity (short moving average) : with mean line
    title = 'acceleration in the direction of gravity (short moving average with mean line)';
    drawWithMeanLine(data_short_moving_average, 'acceleration_gravitaional_direction_moving_average_with_mean_line', title);

    // acceleration in the direction of gravity (short moving average) : with mean line and step counting
    title = 'step counting';
    const positive_threshold = 0.05 / time_interval / skip_length;
    const negative_threshold = -0.05 / time_interval / skip_length;
    drawWithStepCounting(data_short_moving_average, 'step_counting', title, positive_threshold, negative_threshold);
}



function drawMovingAverage(data, moving_average_length, graph_id, title, opening_term, closing_term) {

    const numColumns = data.getNumberOfColumns();
    const numRows = data.getNumberOfRows();

    for (let columnIndex = 1; columnIndex < numColumns; columnIndex++) {
        const data_array = [];

        for (let rowIndex = 0; rowIndex < numRows; rowIndex++) {
            data_array.push(data.getValue(rowIndex, columnIndex));

            if (data_array.length == moving_average_length) {

                let sum = 0;
                for (let i = 0; i < moving_average_length; i++) {
                    sum += data_array[i];
                }
                const moving_average = sum / moving_average_length;

                // failed to update tooltip values
                // data.setValue(rowIndex, columnIndex, moving_average);

                // set 4th parameter (formattedValue) to update tooltip values
                data.setCell(rowIndex, columnIndex, moving_average, moving_average);

                data_array.shift();
            } else {

                // set 4th parameter (formattedValue) to update tooltip values
                data.setCell(rowIndex, columnIndex, 0, 0);
            }
        }
    }

    // remove data written during first 2 seconds and last 3 seconds
    const opening_rows = Math.floor(opening_term / time_interval / skip_length);
    const closing_rows = Math.floor(closing_term / time_interval / skip_length);

    data = remove_opening_and_closing_data(data, opening_rows, closing_rows);

    const options = {title,
                     hAxis: {title: 'time[s]'},
                     vAxis: {title: 'acceleration [m/s^2]'}};
    const chart = new google.visualization.LineChart(document.getElementById(graph_id));
    chart.draw(data, options);
}



function remove_opening_and_closing_data(data, opening_rows, closing_rows) {

    data.removeRows(0, opening_rows);

    const numRows = data.getNumberOfRows();
    const rowIndex_StartClosing = numRows - closing_rows;

    data.removeRows(rowIndex_StartClosing, closing_rows);

    return data;
}



function drawAccelerationProjectedToGravitaionalDirection(data_short_moving_average, data_long_moving_average, graph_id, title, opening_term) {

    const numRows = data_short_moving_average.getNumberOfRows();

    for (let rowIndex = 0; rowIndex < numRows; rowIndex++) {
        const ax = data_short_moving_average.getValue(rowIndex, 1);
        const ay = data_short_moving_average.getValue(rowIndex, 2);
        const az = data_short_moving_average.getValue(rowIndex, 3);

        const gx = data_long_moving_average.getValue(rowIndex, 1);
        const gy = data_long_moving_average.getValue(rowIndex, 2);
        const gz = data_long_moving_average.getValue(rowIndex, 3);
        const absolute_value_of_g = Math.sqrt(gx * gx + gy * gy + gz * gz);

        const accleration_projected_to_garavitational_direction = (ax * gx + ay * gy + az * gz) / absolute_value_of_g;

        // set 4th parameter (formattedValue) to update tooltip values
        data_short_moving_average.setCell(rowIndex, 1, accleration_projected_to_garavitational_direction,
                                          accleration_projected_to_garavitational_direction);
    }

    // set new column label
    data_short_moving_average.setColumnLabel(1, "Gravitational Direction Component");

    // remove other columns
    data_short_moving_average.removeColumns(2, 3);

    // remove data written during first 2 seconds
    const opening_rows = Math.floor(opening_term / time_interval / skip_length);
    data_short_moving_average = remove_opening_and_closing_data(data_short_moving_average, opening_rows, 0);

    const options = {title,
                     hAxis: {title: 'time[s]'},
                     vAxis: {title: 'acceleration [m/s^2]'}};
    const chart = new google.visualization.LineChart(document.getElementById(graph_id));
    chart.draw(data_short_moving_average, options);
}



function drawWithMeanLine(data, graph_id, title) {

    const numRows = data.getNumberOfRows();

    let sum = 0;
    for (let rowIndex = 0; rowIndex < numRows; rowIndex++) {
        sum += data.getValue(rowIndex, 1);
    }
    const mean = sum / numRows;

    data.addColumn('number', 'mean');
    for (let rowIndex = 0; rowIndex < numRows; rowIndex++) {
        data.setCell(rowIndex, 2, mean);
    }

    data.addColumn('number', 'standard gravity (9.80665)');
    for (let rowIndex = 0; rowIndex < numRows; rowIndex++) {
        data.setCell(rowIndex, 3, 9.80665);
    }

    const options = {title,
                     hAxis: {title: 'time[s]'},
                     vAxis: {title: 'acceleration [m/s^2]'}};
    const chart = new google.visualization.LineChart(document.getElementById(graph_id));
    chart.draw(data, options);
}



function drawWithStepCounting(data, graph_id, title, positive_threshold, negative_threshold) {

    const numRows = data.getNumberOfRows();

    data.addColumn('number', 'steps'); // columnIndex 4

    let steps = 0;
    let sum = 0;
    let is_positive = false;
    let is_negative = false;

    for (let rowIndex = 0; rowIndex < numRows; rowIndex++) {

        const diff = data.getValue(rowIndex, 1) - 9.80665;

        if ( (diff > 0 && is_positive == false) || (diff < 0 && is_negative == false) ) {
            sum += diff;
        }

        if (sum > positive_threshold) {
            steps++;
            is_positive = true;
            is_negative = false;
            sum = 0;
        }

        if (sum < negative_threshold) {
            is_positive = false;
            is_negative = true;
            sum = 0;
        }

        data.setCell(rowIndex, 4, steps);
    }


    const options = {title,
                     hAxis: {title: 'time[s]'},
                     vAxes: {
                         0: {title: 'acceleration [m/s^2]'},
                         1: {title: 'steps'}
                     },
                     series: {
                         0: {targetAxisIndex: 0},
                         1: {targetAxisIndex: 0},
                         2: {targetAxisIndex: 0},
                         3: {targetAxisIndex: 1}
                     }};

    const chart = new google.visualization.LineChart(document.getElementById(graph_id));
    chart.draw(data, options);
}

</script>

返信を残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です

CAPTCHA