言語不問で実装する高速な移動平均処理

はじめに

移動平均処理は色々な場面で出てきます。単純に実装すると、愚直に計算する畳み込み積分と同じ計算量になるので、少し工夫して高速化してみます。

移動平均処理とは

基本的な移動平均処理では1サンプルずつ移動しながら平均を求めていきます。画像処理で言えばぼかし処理の一種で、音声処理で言えばローパスフィルタの一種です。
FIRフィルタとして考えると、全ての係数が$\frac{1}{N}$のカーネルによる畳み込み（convolution）です。

以下は単純な実装例です。

#include <stdio.h>
#include <stdlib.h>

int main(void)
{
    #define WINDOW_SIZE (2000)
    #define NUM(fixed_array) ((int)(sizeof(fixed_array) / sizeof(fixed_array[0])))
    #define NUM_PROCESSED (NUM(samples) - WINDOW_SIZE + 1)
    int samples[20000];
    float processed[NUM_PROCESSED];
    int i, j;

    /* 乱数生成 */
    srand(1);
    for (i = 0; i < NUM(samples); i++) {
        samples[i] = rand() % 10;
    }

    /* 移動平均処理 */
    for (i = 0; i < NUM_PROCESSED; i++) {
        int sum = 0;
        for (j = 0; j < WINDOW_SIZE; j++) {
            sum += samples[i + j];
        }
        processed[i] = (float)sum / WINDOW_SIZE;
    }

    /* 結果を表示 */
    for (i = 0; i < NUM_PROCESSED; i++) {
        printf("%f\n", processed[i]);
    }

    return 0;
}

二重ループになっているところが重いです。
WINDOW_SIZE=2000はやりすぎにしても、samplesの要素数20000は少ない方で、音声データなんかだと１秒48000サンプルだったり、画像データだと1920x1080x3ピクセル/フレームだったりします。

高速化のアイディア

平均値の元となる総和に注目すると、窓をスライドするたび、窓から出ていく一番古いサンプルの値を総和から減じて、新しく窓に入ってくるサンプルの値を加えることで更新できそうなことがわかります。

データ構造としてはキュー（Queue）ですね。

でもキューなんて実装していられないので、単純に最初の移動平均値を求めたら、あとは総和を更新しながら処理していく形にしてみます。

#include <stdio.h>
#include <stdlib.h>

int main(void)
{
    #define WINDOW_SIZE (2000)
    #define NUM(fixed_array) ((int)(sizeof(fixed_array) / sizeof(fixed_array[0])))
    #define NUM_PROCESSED (NUM(samples) - WINDOW_SIZE + 1)
    int samples[20000];
    float processed[NUM_PROCESSED];
    int i, j;
    int sum = 0;

    /* 乱数生成 */
    srand(1);
    for (i = 0; i < NUM(samples); i++) {
        samples[i] = rand() % 10;
    }

    /* 移動平均処理 */
    for (j = 0; j < WINDOW_SIZE; j++) {
        sum += samples[j];
    }
    processed[0] = (float)sum / WINDOW_SIZE;

    for (i = 1; i < NUM_PROCESSED; i++) {
        sum -= samples[i - 1];
        sum += samples[i + WINDOW_SIZE - 1];
        processed[i] = (float)sum / WINDOW_SIZE;
    }

    /* 結果を表示 */
    for (i = 0; i < NUM_PROCESSED; i++) {
        printf("%f\n", processed[i]);
    }

    return 0;
}

これなら、最初の窓処理は窓サイズ分の処理負荷を避けられないものの、それ以降は窓サイズに関係なく加算2回と除算1回で移動平均値を求められるので、計算オーダーはほぼ$O(n)$だと言えそうです。

おわりに

処理負荷がどれだけ改善したかは、面倒なので計測しません。2 $<$ WINDOW_SIZEであれば、後者の実装の方が多分速いです。あと、領域外アクセスやコーナーケースについては精査していません。

移動平均処理を畳み込み積分の特殊な場合と考えると、ここで紹介した高速化方法は一般化した畳み込み積分に応用できません。

それでも、移動平均自体はよく用いられる処理なので、特別に高速化する方法を知っていると、いざという時に役立つかもしれません。