我假设一个专用的DSP可以很容易地做到这一点,但是有没有可能不需要那么复杂的“捷径”?
#1 楼
Jason R的答案有一个缺陷,这一点在Knuth的“计算机编程艺术”卷中进行了讨论。 2.如果您的标准偏差只是平均值的一小部分,就会出现问题:E(x ^ 2)-(E(x)^ 2)的计算对浮点舍入误差非常敏感。 />您甚至可以在Python脚本中自己尝试:
ofs = 1e9
A = [ofs+x for x in [1,-1,2,3,0,4.02,5]]
A2 = [x*x for x in A]
(sum(A2)/len(A))-(sum(A)/len(A))**2
我得到-128.0作为答案,这显然在计算上无效,因为数学预测结果应该为非负数。 />
initialize:
m = 0;
S = 0;
n = 0;
for each incoming sample x:
prev_mean = m;
n = n + 1;
m = m + (x-m)/n;
S = S + (x-m)*(x-prev_mean);
,然后在每个步骤之后,
m
的值是平均值,根据您最喜欢的标准偏差定义,标准偏差可以计算为sqrt(S/n)
或sqrt(S/n-1)
。我上面写的方程式与Knuth中的方程式略有不同,但是在计算上是等效的。用Python编写上述公式的代码,并显示您将得到一个非负的答案(希望接近正确的值)。
test1.py:
import math
def stats(x):
n = 0
S = 0.0
m = 0.0
for x_i in x:
n = n + 1
m_prev = m
m = m + (x_i - m) / n
S = S + (x_i - m) * (x_i - m_prev)
return {'mean': m, 'variance': S/n}
def naive_stats(x):
S1 = sum(x)
n = len(x)
S2 = sum([x_i**2 for x_i in x])
return {'mean': S1/n, 'variance': (S2/n - (S1/n)**2) }
x1 = [1,-1,2,3,0,4.02,5]
x2 = [x+1e9 for x in x1]
print "naive_stats:"
print naive_stats(x1)
print naive_stats(x2)
print "stats:"
print stats(x1)
print stats(x2)
结果:
naive_stats:
{'variance': 4.0114775510204073, 'mean': 2.0028571428571427}
{'variance': -128.0, 'mean': 1000000002.0028572}
stats:
{'variance': 4.0114775510204073, 'mean': 2.0028571428571431}
{'variance': 4.0114775868357446, 'mean': 1000000002.0028571}
您会注意到仍然有一些舍入错误,但还不错,而
naive_stats
只是个呕吐。 评论
$ \ begingroup $
+1获取示例代码的详细答案。当需要浮点实现时,这种方法优于我在答案中指出的方法。
$ \ endgroup $
–Jason R
2012年1月20日下午6:46
$ \ begingroup $
也可能会在C ++实现中进行检查:johndcook.com/standard_deviation.html
$ \ endgroup $
–瑞·马克斯(Rui Marques)
13年5月13日在16:43
$ \ begingroup $
是的,就是这样。他使用Knuth使用的精确方程式。您可以进行一些优化,并且避免使用我的方法检查初始迭代和后续迭代。
$ \ endgroup $
–Jason S
13年5月13日在17:41
$ \ begingroup $
“ Knuth引用了一种计算运行均值的方法(我不记得发明者的名字)” –顺便说一下,这是Welford的方法。
$ \ endgroup $
–Jason S
16 Mar 24 '16 at 17:30
$ \ begingroup $
如果有人可以提供帮助,我已经发布了一个与此相关的问题:dsp.stackexchange.com/questions/31812/…
$ \ endgroup $
–乔纳森
16年6月29日在16:34
#2 楼
对于实时应用来说,找到信号平均值和标准偏差的理想方法是什么。我希望能够在某个时间段内,如果信号偏离
平均值超过3个标准偏差,就可以触发
控制器。
在这种情况下,正确的方法通常是计算指数加权的运行平均值和标准偏差。在指数加权平均值中,均值和方差的估计值偏向最近的样本,从而为您提供过去$ \ tau $秒内的均值和方差的估计值,这可能是您想要的,而不是通常的算术平均值在有史以来的所有样本中。
在频域中,“指数加权移动平均值”只是一个真实的极点。在时域中实现很简单。
时域实现
让
mean
和meansq
是信号均值和均方根的当前估计。在每个周期上,请使用新样本x
更新这些估计值:$ 0 上面表示为命令式程序的内容也可能被描述为信号流程图:
分析
上述算法计算$ y_i = a x_i +(1-a)y_ {i-1} $,其中$ x_i $是样本$ i的输入$,而$ y_i $是输出(即均值的估计值)。这是一个简单的单极IIR滤波器。通过$ z $变换,我们发现传递函数$$ H(z)= \ frac {a} {1-(1-a)z ^ {-1}} $$。
现在将IIR滤波器压缩到自己的块中,该图现在看起来像这样:
为了进入连续域,我们进行替换$ z = e ^ {s T} $其中$ T $是采样时间,$ f_s = 1 / T $是采样率。求解$ 1-(1-a)e ^ {-sT} = 0 $,我们发现连续系统在$ s = \ frac {1} {T} \ log(1-a)$处有一个极点。
选择$ a $:$$ a = 1-\ exp \ left \ {2 \ pi \ frac {T} {\ tau} \ right \} $$
参考文献
Simulink图源可以从https://gist.github.com/1942771
下载
评论
$ \ begingroup $
您能解释一下$ a $如何确定移动平均线的长度吗?以及应该使用$ a $的哪个值?规格0> a> 1是无法满足的。
$ \ endgroup $
– Dilip Sarwate
2012年1月20日在2:02
$ \ begingroup $
这类似于Jason R的方法。此方法将不太准确,但会更快一些并且内存更少。这种方法最终使用了指数窗口。
$ \ endgroup $
– schnarf
2012年1月20日4:44
$ \ begingroup $
糟糕!当然,我的意思是0 $ \ endgroup $
– nibot
2012年2月29日在20:51
$ \ begingroup $
我认为这里可能有错误。单极滤波器在DC处没有0 dB的增益,并且由于您在线性域中应用一个滤波器,在平方域中应用一个滤波器,因此对于E
$ \ endgroup $
–希尔马
2012年3月1日13:22
$ \ begingroup $
在DC处确实具有0 dB的增益。将z = 1(DC)代入H(z)= a /(1-(1-a)/ z),得到1。
$ \ endgroup $
– nibot
2012年1月1日13:27
#3 楼
我以前在嵌入式处理应用程序中使用的一种方法是维护感兴趣信号的和和平方和的累加器:$$
A_ {x,i} = \ sum_ {k = 0} ^ {i} x [k] = A_ {x,i-1} + x [i],A_ {x,-1} = 0
$$
$$
A_ {x ^ 2,i} = \ sum_ {k = 0} ^ {i} x ^ 2 [k] = A_ {x ^ 2,i-1} + x ^ 2 [i],A_ {x ^ 2,-1} = 0
$$
还要在上述等式中跟踪当前时刻$ i $(即,请注意已添加到累加器中的样本数量)。然后,在时间$ i $处的样本均值和标准差为:
$$
\ tilde \ mu = \ frac {A_ {x_i}} {i + 1}
$$
$$
\ tilde \ sigma = \ sqrt {\ frac {A_ {x ^ 2_i}} {i + 1}-\ tilde \ mu ^ 2}
$$
或者您可以使用:
$$
\ tilde \ sigma = \ sqrt {\ frac {A_ {x ^ 2_i}}} {i}-\ tilde \ mu ^ 2}
$$$
取决于您喜欢哪种标准偏差估计方法。这些方程式基于方差的定义:
$$
\ sigma ^ 2 = \ operatorname {E}(X ^ 2)-(\ operatorname {E}(X))^ 2
$$
我过去曾经成功使用过这些方法(尽管我只关注方差估计,而不关注标准差),尽管您必须谨慎使用数字类型如果您要长时间累加,可使用它来保持累加器;您不希望溢出。
编辑:除了上面关于溢出的注释之外,应该注意的是,在浮点算术中实现时,这不是一种数值健壮的算法。在估算的统计数据中。在这种情况下,请查看Jason S的答案以获得更好的方法。
评论
$ \ begingroup $
也许将其重写为$ A_ {x,i} = x [i] + A_ {x,i-1},\ A_ {x,0} = x [0] $会清楚表明您需要在每个步骤中仅添加两个数字,以免有人将其实现为在每个步骤中将$ x $的所有$ i $元素相加。
$ \ endgroup $
–乳香
2011年12月6日15:18
$ \ begingroup $
是的,这样更好。我试图重写以使递归实现更加清晰。
$ \ endgroup $
–Jason R
2011年12月6日在16:12
$ \ begingroup $
-1当我有足够的代表这样做时:这存在数值问题。见Knuth卷。 2
$ \ endgroup $
–Jason S
2012年1月20日,0:23
$ \ begingroup $
这里似乎有一些错别字。为什么在平方根符号下减去$ \ sigma $的平均值?它应该是$ \ mu ^ 2 $才能匹配显示的方程$ \ sigma ^ 2 = E(X ^ 2)-(E(X))^ 2 $,不是吗?另外,尽管我不会拒绝这个答案,但我同意Jason S的观点,认为这种方法可能存在数值问题。
$ \ endgroup $
– Dilip Sarwate
2012年1月20日下午1:33
$ \ begingroup $
@JasonS:我不同意该技术固有的缺陷,尽管我同意您的观点,即在浮点数中实现时,它不是一种数值上可靠的方法。我应该更清楚地知道,我已经在使用整数算术的应用程序中成功使用了它。整数(或分数的定点实现)算术不会遭受您指出的会导致精度损失的问题的困扰。在这种情况下,这是一种适合的方法,每个样品需要较少的操作。
$ \ endgroup $
–Jason R
2012年1月20日6:45
#4 楼
Jason和Nibot的回答在一个重要方面有所不同:Jason的方法计算整个信号的std dev和均值(因为y = 0),而Nibot的计算是“运行中”的计算,即,它权衡比来自样本的样本更强的最近样本遥远的过去。由于该应用程序似乎需要std dev,并且作为时间的函数,因此Nibot的方法可能更合适(对于此特定应用程序)。但是,真正棘手的部分是正确设置时间加权部分。 Nibot的示例使用一个简单的单极点滤波器。
描述此问题的正确方法是,通过过滤$ x [n] $得到$ E [x] $的估计,通过过滤得到$ E [x ^ 2] $的估计$ x [n] ^ 2 $。估计滤波器通常是低通滤波器。这些滤波器应定标为在0 Hz时增益为0dB。否则会有恒定的增益误差。
对信号的了解以及估计所需的时间分辨率可指导选择低通滤波器。较低的截止频率和较高的阶数将导致更好的精度,但响应时间变慢。
要使事情复杂化,在线性域中进一步应用一个滤波器,在平方域中应用另一个。平方会显着改变信号的频谱内容,因此您可能需要在平方域中使用其他滤波器。
下面是一个示例,说明如何根据时间估算均值,均方根和标准差。
%% example
fs = 44100; n = fs; % 44.1 kHz sample rate, 1 second
% signal: white noise plus a low frequency drift at 5 Hz)
x = randn(n,1) + sin(2*pi*(0:n-1)'*5/fs);
% mean estimation filter: since we are looking for effects in the 5 Hz range we use maybe a
% 25 Hz filter, 2nd order so it's not too sluggish
[b,a] = butter(2,25*2/fs);
xmeanEst = filter(b,a,x);
% now we estimate x^2, since most frequency double we use twice the bandwidth
[b,a] = butter(2,50*2/fs);
x2Est = filter(b,a,x.^2);
% std deviation estimate
xstd = sqrt(x2Est)-xmeanEst;
% and plot it
h = plot([x, xmeanEst sqrt(x2Est) xstd]);
grid on;
legend('x','E<x>','sqrt(E<x^2>)','Std dev');
set(h(2:4),'Linewidth',2);
评论
$ \ begingroup $
我的答案中的过滤器对应于y1 = filter(a,[1(1-a)],x);。
$ \ endgroup $
– nibot
2012年3月1日14:05
$ \ begingroup $
区分运行统计数据和总体样本统计数据之间的区别。我的实现可以修改为通过在移动窗口上进行累加来计算运行统计信息,这也可以高效地完成(在每个时间步,从每个累加器中减去刚从窗口中滑出的时间样本)。
$ \ endgroup $
–Jason R
2012年3月1日下午14:13
$ \ begingroup $
nibot,对不起,您是对的,我错了。我会马上纠正
$ \ endgroup $
–希尔马
2012年3月1日17:52
$ \ begingroup $
+1建议对x和x ^ 2进行不同的过滤
$ \ endgroup $
– nibot
2012年3月1日于20:57
#5 楼
与上面的首选答案(Jason S.)类似,并且也从取自Knut的公式(第2卷,第232页)中得出,也可以导出一个公式来替换值,即一步删除并添加值。根据我的测试,替换提供了比两步删除/添加版本更高的精度。和在Jason的帖子中的mean
。值s
是指窗口大小m
。/**
* Replaces the value {@code x} currently present in this sample with the
* new value {@code y}. In a sliding window, {@code x} is the value that
* drops out and {@code y} is the new value entering the window. The sample
* count remains constant with this operation.
*
* @param x
* the value to remove
* @param y
* the value to add
*/
public void replace(double x, double y) {
final double deltaYX = y - x;
final double deltaX = x - mean;
final double deltaY = y - mean;
mean = mean + deltaYX / count;
final double deltaYp = y - mean;
final double countMinus1 = count - 1;
s = s - count / countMinus1 * (deltaX * deltaX - deltaY * deltaYp) - deltaYX * deltaYp / countMinus1;
}
评论
您对信号一无所知吗?它是静止的吗?@Tim假设它是固定的。出于我自己的好奇,非平稳信号的后果是什么?
如果它是固定的,则可以简单地计算运行平均值和标准偏差。如果平均值和标准偏差随时间变化,事情将会变得更加复杂。
密切相关:en.wikipedia.org/wiki/…