平均値

リアルタイムに入ってくる信号の平均値を求めたい

MORE...

実現方法

総和や平均値をフィードバック回路で計算する方法

MORE...

注意点

直流を入力し続ける場合は倍精度に

MORE...

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

平均値

入力信号 \(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\) を計算する回路
この出力を \(N+1\) で割れば平均値 \(\bar{x}_N\) が得られます。


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

C++ コーディング例を以下に示します。

	float Average1(const float x)
	{
		static int N=0;
		static float z=0;
		return((z+=x)/(++N));
	}
				
この回路に 200 サンプルで 1 周期の正弦波を入れたときの出力を図3に示します。

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

波打ちながらも次第に真の平均値 0 に向かってゆきます。

その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 回で計算が合わなくなるのは、ちょっと意外ではないでしょうか。

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