リアルタイムに平均値を求める

平均値

入力信号 \(x_n\) の平均値 \begin{eqnarray} \bar{x}_N &=& \frac{1}{N+1}\sum_{n=0}^N x_n\label{barxN} \end{eqnarray}

実現方法

その1

総和部分を \begin{eqnarray} z_N &=& \sum_{n=0}^N x_n \end{eqnarray} とすると \begin{eqnarray} z_N &=& \left(\sum_{n=0}^{N-1} x_n\right) + x_N \\ &=& z_{N-1} + x_N \end{eqnarray} と書けますので、この部分はフィードバック・ループを含む次のような回路で実現できます。

図1 : \(z_N=\displaystyle\sum_{n=0}^N x_n\) を計算する回路

図2 : 平均値 \(\bar{x}_N=\displaystyle\frac{1}{N+1}\sum_{n=0}^N x_n\) を計算する回路

	float Average1(const float x)
	{
		static int N=0;
		static float z=0;
		return((z+=x)/(++N));
	}
				

図3 : Average1( ) で正弦波の平均値を求めてみる
(横軸:サンプル番号 n, 縦軸:入出力値)

その2

式(\ref{barxN})から 1 サンプル前までの平均値は \begin{eqnarray} \bar{x}_{N-1} &=& \frac{1}{N}\sum_{n=0}^{N-1} x_n \end{eqnarray} と表せますので、両辺を N 倍して左右を入れ替えれば \begin{eqnarray} \sum_{n=0}^{N-1} x_n &=& N \cdot\bar{x}_{N-1} \end{eqnarray} です。これを使うと \begin{eqnarray} \bar{x}_{N} &=& \frac{1}{N+1}\sum_{n=0}^N x_n \\ &=& \frac{1}{N+1}\left(x_N + \sum_{n=0}^{N-1} x_n\right) \\ &=& \frac{1}{N+1}\left(x_N + N \cdot\bar{x}_{N-1}\right) \\ \end{eqnarray} と書き換えることができます。浮動小数点演算には誤差がありますので、完全に等価とはいえませんが、次の回路でも同様の結果が得られます。

図4 : 等価回路

C++ で書けば次のようになります。

	float Average2(const float x)
	{
		static int N = 0;
		static float z = 0;
		z = x + (N++)*z;
		return(z/=N); // 上で N++ してあるので /=N でよい
	}
			

注意点

浮動小数点には「情報落ち」や「丸め」による誤差があるため、入力信号に直流成分があると、図2の方法でも図4の方法でも、 時間が経つと急に誤差が大きくなります。

わかりやすい例として、以下のコードにより図2の方法で作られた Average1( ) に 32767 を与えて繰り返しコールした時の出力を図5に示します。

	ofstream o("out.txt", ios::out);
	o.precision(12);
	for (int n = 0; n <= 32768; n++)
	{
		float a = Average1(32767);
		o << n << "	" << a << endl;
	}
			

図5 : Average1(32767) を繰り返しコールした時の返り値
(横軸:サンプル番号 n, 縦軸:返り値)

引数に与えているのは定数値 32767 ですから、何度コールしても結果(平均値)は 32767 になってほしいところですが、 Average1(32767) の返り値が正しいのは n=511 までで、n=512 では 32767.0019531 が返って来てしまい、その後は誤差が増大しています (誤差の出方は入力する値によって大きく異なります)。

32767 は 2 進法では \(111\ 1111\ 1111\ 1111_2\) と 15bit ですから、512 (2 進法で \(10\ 0000\ 0000_2\) と 10bit) 回足し続けると、結果が 25bit になってしまい、 単精度浮動小数点の仮数部 24bit を超えてしまうため、誤差が出始めるのです。

図6 : フィードバック・ループ内の遅延レジスタの値が大きくなりすぎる

図4の方法で作られた Average2( ) でも同様です。 こちらは遅延レジスタに記憶される数値は大きくならないのですが、N 倍して入力と足し合わせるため、やはり仮数部の桁が足りなくなってしまうのです。

図7 : フィードバックが N 倍されるので、入力との加算で情報落ちが発生する

何百万回もコールして誤差が出るというのなら仕方ないと諦めもつきますが、単精度浮動小数点とはいえ、たった 513 回で計算が合わなくなるのは、ちょっと意外ではないでしょうか。

【注意】直流成分がある信号の総和や平均値を浮動小数点で求める時は、仮数部の桁あふれに注意する必要があります。 特に単精度浮動小数点の場合、サンプリング周波数が高いと、あっという間に桁あふれします。 固定小数点の場合も同様です。 長時間平均を計算し続ける場合は、倍精度浮動小数点の使用をお勧め致します。