《Delta-Sigma数据转换器从入门到精通》笔记之漫谈均值

《Delta-Sigma数据转换器从入门到精通》笔记之漫谈均值

---

前言

“均值”是δ-σ ADC CIC滤波器的重要概念，深刻地理解“均值”是理解δ-σ ADC 以及其他数据处理方法的前提。《Delta-Sigma数据转换器从入门到精通》反复涉及到这一概念，但是作者并未做深入的介绍，本文尝试着从多个视角对此加以解释。

一、最小方差无偏估计

“均值”应该是所有数据处理方法中最常用的也是最朴素的了。试想以下一个情形，一位工程师被要求设计一个测量电池电压（VBAT）系统，这位工程师很快地交出了答卷——连续取样N次，求样本均值，即：
$$\mathbf{X}\text{=}{\left[\begin{array}{ccccc}x[1]& x[2]& x[3]& ...& x[N]\end{array}\right]}^T$$
$$\hat{x}=\overline{x}=\frac{1}{N}\sum _{\text{n}=1}^N\text{x}[n]\text{\hspace{1em}(1)}$$

问题来了，取样一次不就足够了吗，为什么要取多次？$\hat{x}$为什么不是从$\mathbf{X}$中随机取一个$x[\text{n}]$ ，或者是$x[\text{n}]$的其他线性组合$\hat{x}=\sum _{\text{n}=1}^N{\text{a}}_{\text{n}}x[n]$?糟糕的工程师会这样解释：这是祖传的，经过实践的，运行的好好的，又没出过问题。一般的工程师会这样解释：测量难免有误差，多测几次才准，至于取均值，也很简单， $x[\text{n}]$之间都是平等的，不能厚此薄彼。但是有追求的工程师不满足于这样朴素的回答，我们对此问题进行建模，围绕着模型展开论述。
$$x[\text{n}]\text{=}\theta \text{+}\omega [\text{n}]\begin{array}{cccc}& & & \end{array}\text{n=1,2,}...N\text{\hspace{1em}(2)}$$
其中，$\theta \begin{array}{c}\end{array}\text{(}\theta \in \text{(-}\infty \text{,+}\infty \text{))}$是要估计（通俗的说，要测量的）参数，误差$\omega [\text{n}]$被假设成WGN（高斯白噪声），即$\omega [\text{n}]\sim N\text{(0,}{\sigma }^{\text{2}}\text{)}$，这一假设是符合工程师经验的，工程上遇到的许多误差项都可以近似成WGN，有的难以得知误差的CDF，干脆假设成gauss好了。
这一模型，过于理想，因为随着时间的推移 $\theta$ 总是会减小的。但是，还好在足够短的时间内，模型是够用的。
工程师们常常认为数据是上帝的，但是有的工程师认为，一旦模型建立，论述就要以模型为中心了，这时候，模型就是上帝。我赞同这种观点。
接着式（1），$\theta$ 的一个估计$x[\text{n}]$的均值
$$\hat{\theta }=\frac{1}{N}\sum _{\text{n}=1}^{N}x[n]\text{\hspace{1em}(3)}$$
求$\hat{\theta }$期望
$$\text{E(}\hat{\theta })=\text{E(}\frac{1}{N}\sum _{\text{n}=1}^{N}x[n])$$
$$\text{E(}\hat{\theta })=\frac{1}{N}\sum _{\text{n}=1}^N\text{E}(x[n])$$
$$\text{E(}\hat{\theta })=\frac{1}{N}\sum _{\text{n}=1}^N\theta$$
$$\text{E(}\hat{\theta })=\theta \text{\hspace{1em}(4)}$$

求$\hat{\theta }$方差
$$\text{var(}\hat{\theta })=\mathrm{var}\text{(}\frac{1}{N}\sum _{\text{n}=1}^{N}x[n])$$
$$\text{var(}\hat{\theta })=\frac{1}{N^2}\sum _{\text{n}=1}^N\text{var}(x[n])$$
$$\text{var(}\hat{\theta })=\frac{1}{N^2}\sum _{\text{n}=1}^N{\sigma }^2$$
$$\text{var(}\hat{\theta })=\frac{{\sigma }^2}{N}\text{\hspace{1em}(5)}$$
显然的，$\hat{\theta }$是$\theta$ 的一个无偏估计量，并且方差仅为单次取样的1/N，当$\text{N}\to \infty$，$\text{var(}\hat{\theta })\to 0$。$\theta$ 的无偏估计量当然还有很多，任意的
$$\tilde{\theta }=\sum _{\text{n}=1}^N{\text{a}}_{\text{n}}x[n]\begin{array}{cc}& \end{array}(\sum _{\text{n}=1}^N{\text{a}}_{\text{n}}\text{=}1)$$
都是无偏的。那么$\hat{\theta }$是最优的吗？如果以最小方差作为判断依据，答案是肯定的，因为它达到了Cramer-Rao下界。
首先求联合PDF
$$p(\mathbf{X};\theta )=\prod _{n=1}^N\frac{1}{\sqrt{2\pi {\sigma }^2}}{e}^{-\frac{{\left(x[n]-\theta \right)}^2}{2{\sigma }^2}}$$
$$p(\mathbf{X};\theta )=\frac{1}{{\left(2\pi {\sigma }^2\right)}^{\frac{N}{2}}}{e}^{-\frac{\sum _{\text{n=}1}^{\text{N}}{\left(x[n]-\theta \right)}^2}{2{\sigma }^2}}$$
取对数求一阶导
$$\frac{\partial \ln p(\mathbf{X};\theta )}{\partial \theta }=\frac{\partial }{\partial \theta }\left(\text{-}\ln \left[{\left(2\pi {\sigma }^2\right)}^{\frac{N}{2}}\right]-\frac{\sum _{\text{n=}1}^{\text{N}}{\left(x[n]-\theta \right)}^2}{2{\sigma }^2}\right)$$
$$\frac{\partial \ln p(\mathbf{X};\theta )}{\partial \theta }=\frac{\sum _{\text{n=}1}^{\text{N}}\left(x[n]-\theta \right)}{{\sigma }^2}$$
$$\frac{\partial \ln p(\mathbf{X};\theta )}{\partial \theta }=\frac{N\left(\overline{x}-\theta \right)}{{\sigma }^2}$$
$$\frac{\partial \ln p(\mathbf{X};\theta )}{\partial \theta }=\frac{N\overline{x}}{{\sigma }^2}\text{-}\frac{\text{N}\theta }{{\sigma }^2}$$

求二阶导
$$\frac{{\partial }^2\ln p(\mathbf{X};\theta )}{\partial {\theta }^2}=\frac{\partial }{\partial \theta }\left(\frac{N\overline{x}}{{\sigma }^2}\text{-}\frac{\text{N}\theta }{{\sigma }^2}\right)$$
$$\frac{{\partial }^2\ln p(\mathbf{X};\theta )}{\partial {\theta }^2}=\text{-}\frac{\text{N}}{{\sigma }^2}$$
即
$$\frac{1}{\text{-E }[\frac{{\partial }^2\ln p(\mathbf{X};\theta )}{\partial {\theta }^2}]}=\frac{{\sigma }^2}{\text{N}}$$
据此，均值作为该模型参数$\theta$的最佳估计得到了充分的证明，对于一名工程师而言，可以放心地使用了。

二、FIR滤波器

方差减小，意味着直观上，输出更为平滑了，暗示着我们均值是一个低通滤波器。审视式（2）给出的模型，可以将噪声$\omega [\text{n}]$看成交流分量，而$\theta$ 是直流分量，显然，我们需要一个低通滤波器提取$\theta$。仔细审视式（1），可以发现它同时就是FIR滤波器，其Z变换为
$$\text{H(Z)=}\sum _{n=0}^{N-1}\frac{1}{\text{N}}{Z}^{-n}$$
$$\text{H(Z)=}\frac{1}{\text{N}}\sum _{n=0}^{N-1}{Z}^{-n}\text{\hspace{1em}(6)}$$
当N确定时，利用MATLAB可以画出幅频/相频特性曲线。

工程师们如果对你说，均值可以对数据进行滤波，输出更为平滑，请不要怀疑，他们是在用通俗的语言在向你解释FIR滤波器。
和FIR相对应的，还有一个叫IIR的滤波器，其差分方程如下式。
$$y(n)+\sum _{p=1}^P{a}_{p}y(n-p)=\sum _{q=0}^Q{b}_{q}x(n-q)$$
$$y(n)=\sum _{p=1}^P\text{-}{a}_{p}y(n-p)\text{+}\sum _{q=0}^Q{b}_{q}x(n-q)\text{\hspace{1em}(7)}$$
Z变换可得
$$\text{H(Z)}=\frac{\sum _{q=0}^Q{b}_q{\text{Z}}^{\text{-}q}}{1\text{-}\sum _{p=1}^P{a}_p{\text{Z}}^{\text{-}p}}$$

对于式（7），我常常喜欢把它写成
$$y(n)=\sum _{p=1}^P{a}_{p}y(n-p)\text{+}\sum _{q=0}^Q{b}_{q}x(n-q)$$
前者是差分方程的理解，后者意味着输出是历史输入和历史输出的加权组合，ARM公司提供的库也是这样设计的。当然，这只是个人理解，更多的场合，还是使用式（7）更为方便。相比于FIR，IIR意味着什么呢？多了分母，意味着反馈？很多工程师就是这么理解的。我们重新审视式（6），对其乘以（1-Z-1）再除以（1-Z-1）
$$\text{H(Z)=}\frac{1}{\text{N}}\frac{\left(\sum _{n=0}^{N-1}{Z}^{-n}\right)\left(1-Z^{-1}\right)}{1-Z^{-1}}$$
$$\text{H(Z)=}\frac{1}{\text{N}}\frac{\left(1+Z^{-1}+Z^{-1}+...+Z^{-(N-1)}\right)\left(1-Z^{-1}\right)}{1-Z^{-1}}$$
$$\text{H(Z)=}\frac{1}{\text{N}}\frac{1-Z^{-N}}{1-Z^{-1}}\text{\hspace{1em}(8)}$$
FIR便写成了IIR的形式。用MATLAB绘制一下N=10的幅频/相频特性曲线

完全没有变化，只是matlab也误以为它是IIR

由此工程师们应该明白，判断一个滤波器到底是FIR还是IIR，其Z变换是否有分母并不是一个判据。
式（6）是均值滤波（具体的讲，是滑动均值滤波），和它等价的式（8）可以理解成输入累加（分母）减去N项前的输入累加（分子），即数据先经过一个累加器，再经过一个减去N阶差分的减法器，也即式（8）进一步变换成下式
$$\text{H(Z)=}\frac{1}{\text{N}}\times \frac{1}{1-Z^{-1}}\times \left(1-Z^{-N}\right)$$
$$\text{H(Z)=}\frac{1}{\text{N}}\times {\text{H}}_{\text{1}}\left(Z\right)\times {\text{H}}_{\text{2}}\left(Z\right)$$
当然，也可以反过来理解，先经过减法器后经过累加器，即
$$\text{H(Z)=}\frac{1}{\text{N}}\times {\text{H}}_{\text{2}}\left(Z\right)\times {\text{H}}_{\text{1}}\left(Z\right)$$
理解了这一点，δ-σ ADC 中的CIC滤波器就可以理解一半了。再次审视式（6），它有一项“$\frac{1}{\text{N}}$”，我们在阅读ADI提供的PDM麦克风驱动源码的时候发现，并没有这一项，这不稀奇，只要把公式乘以N，等效成输出放大N倍即可，这样做的好处是没有了除法，不会出现浮点运算带来的误差累积，也方便后续的各种定点滤波。

三、最小二乘估计

我们再次审视式（2）所给出的模型，概率论学得好的工程师会发现这是一个平稳随机过程，LS（最小二乘）和ML（最大似然）的方法便在我们脑海中产生。本节先介绍LS。
构造目标函数
$$\text{J(}\theta \text{)=}\sum _{\text{n=}1}^{\text{N}}{\left(\text{x }[\text{ n }]\text{ -}\theta \right)}^2$$
$$\text{J(}\theta \text{)=}\sum _{\text{n=}1}^{\text{N}}\left(\text{x }[\text{ n}{]}^{\text{2}}\text{-2x }[\text{ n }]\theta +{\theta }^2\right)$$
可见是一个二阶函数，可以求极值。令
$$\frac{\partial \text{J(}\theta \text{)}}{\partial \theta }\text{=}0$$
即
$$\frac{\partial }{\theta }\sum _{\text{n=}1}^{\text{N}}\left(\text{x }[\text{ n}{]}^{\text{2}}\text{-2x }[\text{ n }]\theta +{\theta }^2\right)\text{=}0$$
$$\sum _{\text{n=}1}^{\text{N}}\left(\text{-2x }[\text{ n }]+2\theta \right)\text{=}0$$
$$\text{N}\theta \text{=}\sum _{\text{n=}1}^{\text{N}}\text{x }[\text{ n }]$$
$$\theta \text{=}\frac{\text{1}}{\text{N}}\sum _{\text{n=}1}^{\text{N}}\text{x }[\text{ n }]$$

四、最大似然估计

ML的方法和上文求Cramer-Rao下界差不多，许多推导步骤因此跳过。由于$\omega [\text{n}]\sim N\text{(0,}{\sigma }^{\text{2}}\text{)}$，因此$x[n]$的联合PDF为
$$p(\mathbf{X};\theta )=\prod _{n=1}^N\frac{1}{\sqrt{2\pi {\sigma }^2}}{e}^{-\frac{{\left(x[n]-\theta \right)}^2}{2{\sigma }^2}}$$
$$p(\mathbf{X};\theta )=\frac{1}{{\left(2\pi {\sigma }^2\right)}^{\frac{N}{2}}}{e}^{-\frac{\sum _{\text{n=}1}^{\text{N}}{\left(x[n]-\theta \right)}^2}{2{\sigma }^2}}$$
取对数并且令其一阶导为0
$$\frac{\partial \ln p(\mathbf{X};\theta )}{\partial \theta }=0$$
$$\frac{N\overline{x}}{{\sigma }^2}\text{-}\frac{\text{N}\theta }{{\sigma }^2}\text{=}0$$
$$\theta \text{=}\overline{\text{x}}\text{=}\frac{\text{1}}{\text{N}}\sum _{\text{n=}1}^{\text{N}}\text{x }[\text{ n }]$$

五、几何

我不是一个聪明的工程师，常常难以接受纯粹的公式推导，总是试图探究其背后的物理意义，特别是其几何意义。也正是这个原因，我更喜欢“均值”、“方差”的另一种充满直观的称呼：“一阶原点矩”、“二阶中心距”，也喜欢将方差看成正方形，将协方差看成矩形。也只有借助几何的概念，我才能理解微分方程的通解，特解，齐次解、理解傅里叶变换、理解Hilbert空间、理解线性代数、理解…。确实是笨了点，这或许注定了我只是一名普通的工程师，为了一口饭四处奔波。
在这里，我也从几何的视角解释均值。考察观测向量
$$\mathbf{X}\text{=}{\left[\begin{array}{ccccc}x[1]& x[2]& x[3]& ...& x[N]\end{array}\right]}^T$$
和目标向量
$$\theta \text{=}{\left[\begin{array}{ccccc}\theta & \theta & \theta & ...& \theta \end{array}\right]}^T$$
$$\theta \text{=}\theta {\left[\begin{array}{ccccc}1& 1& 1& ...& 1\end{array}\right]}^T$$
在Euclidean空间${\text{R}}^{\text{N}}$当中，向量$\theta$张成了1维子空间$\text{span}(\theta )\in {\text{R}}^{\text{N}}$，我们的目标是在该空间当中寻找一个向量$\hat{\theta }$，该向量最接近$\mathbf{X}$。所谓最接近，是指误差向量范数$\Vert \mathbf{e}\Vert \text{=}\Vert \mathbf{X}\text{-}\hat{\theta }\Vert$最小，显然的，只有span($\theta$)和$\mathbf{e}$正交时，$\Vert \mathbf{e}\Vert$最小，这时候$\hat{\theta }$就是$\mathbf{X}$在span($\theta$)上的投影，一个内积公式便可以解决了。
首先求向量内积。
$$<\theta ,\mathbf{X}>={\theta }^T\mathbf{X}$$
$$<\theta ,\mathbf{X}>=\theta \left(\text{x }[\text{ 1 }]\text{ +x }[\text{ 2 }]\text{ +}...\text{+x }[\text{ N }]\right)$$
$$<\theta ,\mathbf{X}>=\theta \sum _{n=1}^N\text{x }[\text{ n }]$$

由于$\Vert \theta \Vert \ne 1$，因此需要做单位化
$$<\theta ,\theta >={\theta }^T\theta \text{=N}{\theta }^2$$
由此可得
$$\hat{\theta }\text{=}\theta \frac{<\theta ,\mathbf{X}>}{<\theta ,\theta >}$$
$$\hat{\theta }\text{=}\theta {\left[\begin{array}{ccccc}1& 1& 1& ...& 1\end{array}\right]}^T\frac{\theta \sum _{n=1}^N\text{x }[\text{ n }]}{\text{N}{\theta }^2}$$
$$\hat{\theta }\text{=}\frac{1}{\text{N}}\sum _{n=1}^N\text{x }[\text{ n }]{\left[\begin{array}{ccccc}1& 1& 1& ...& 1\end{array}\right]}^T$$
嘿，真形象、真直观、真简单、真有趣！

六、协方差

我们再把视线放回δ-σ ADC的滑动均值滤波中。正如上文所述，滤波器输出变得平滑了，因为输出序列的每一项方差都为$\frac{{\sigma }^2}{N}$，每一项都更可能地接近$\theta$，我们自然地想到，输出序列项之间相互更为接近了。这当然是朴素的想法，作为一名工程师，应该学会用数学的语言去阐述。
对于任意的s和t，$\omega [\text{n}]$序列的自协方差
$${\gamma }_{\omega }(s,t)=\mathrm{cov}({\omega }_s,{\omega }_t)$$
$${\gamma }_{\omega }(s,t)=\{\begin{array}{c}{\sigma }^2\text{ s=t}\\ 0\text{ s}\ne \text{t}\end{array}$$
滤波后的序列
$${\gamma }_{\overline{\omega }}(s,t)=\mathrm{cov}({\omega }_{\overline{s}},{\omega }_{\overline{\text{t}}})$$
$${\gamma }_{\overline{\omega }}(s,t)=\mathrm{cov}(\frac{1}{\text{N}}\sum _{\text{n=}0}^{\text{N-1}}{\omega }_{\text{s-n}},\frac{1}{\text{N}}\sum _{\text{n=}0}^{\text{N-1}}{\omega }_{\text{t-n}})$$
$${\gamma }_{\overline{\omega }}(s,t)=\frac{1}{{\text{N}}^{\text{2}}}\mathrm{cov}(\sum _{\text{n=}0}^{\text{N-1}}{\omega }_{\text{s-n}},\sum _{\text{n=}0}^{\text{N-1}}{\omega }_{\text{t-n}})$$
$${\gamma }_{\overline{\omega }}(s,t)=\frac{1}{{\text{N}}^{\text{2}}}E\left(\sum _{\text{n=}0}^{\text{N-1}}{\omega }_{\text{s-n}}\times \sum _{\text{n=}0}^{\text{N-1}}{\omega }_{\text{t-n}}\right)$$
$${\gamma }_{\overline{\omega }}(s,t)=\frac{1}{{\text{N}}^{\text{2}}}E\left[\left({\omega }_{\text{s}}\text{+}{\omega }_{\text{s-1}}\text{+}...\text{+}{\omega }_{\text{s-(N-1)}}\right)\left({\omega }_{\text{t}}\text{+}{\omega }_{\text{t-1}}\text{+}...\text{+}{\omega }_{\text{t-(N-1)}}\right)\right]$$
$${\gamma }_{\overline{\omega }}(s,t)=\{\begin{array}{c}\frac{1}{{\text{N}}^{\text{2}}}\sum _{\text{n=}0}^{\text{N-1-}|\text{s-t}|}{\sigma }^2\begin{array}{cc}& \end{array}\left|s-t\right|\le N-1\\ 0\begin{array}{cccc}& & & \end{array}\left|s-t\right|>N-1\end{array}$$
由式(9)可知，对于$\left|s-t\right|\le N-1$，$\overline{\omega }[\text{n}]$序列的自协方差${\gamma }_{\overline{\omega }}(s,t)$＞0，而且$\left|s-t\right|$越小，${\gamma }_{\overline{\omega }}(s,t)$越大，越相关，自然地，$\overline{\omega }[\text{n}]$相邻项之间越接近，直观的来看，输出越平滑。

写在最后

夜深了，该结束了，而“均值”的故事还远没有结束，比如上文提及的$\theta \in \text{(-}\infty \text{,+}\infty \text{)}$，如果我们先验地认为$\theta \in \text{(a,b)}$，会怎么样呢？手机上的锂离子电池电压，它的值就可以认为是$\theta \in [0\text{,4}\text{.2 }]$，这种条件下数据该如何处理？，工程师们常常使用的中值滤波、去除离群值后再取均值，它的统计参数又该发生什么变化？
这些和δ-σ ADC的关系就不太大了，等《Delta-Sigma数据转换器从入门到精通》系列笔记完成后找机会再聊！
晚安，迟哥。