《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$$
嘿，真形象、真直观、真简单、真有趣！

六、最小二乘的几何解释

最小二乘有对应的几何解释，对于我们这些搞振动信号分析的人来说算是常识了，上文“几何”法是一维的，可以直接列出投影向量的公式，在高维问题中，这么做复杂了。高维问题中常常是利用误差向量和子空间的任意向量正交的特点列方程，求解。这些是线性代数基础知识，基础得很，不做过多解释了，套套公式吧。
令$\mathbf{U}={\left[\begin{array}{cccc}1& 1& ...& 1\end{array}\right]}^T\in R^N$，误差向量$\mathbf{X}-\mathbf{U}\theta$和$\forall \gamma \in span(\mathbf{U})$正交，显然的
$$⟨\mathbf{U},\mathbf{X}-\mathbf{U}\theta ⟩={\mathbf{U}}^T\left(\mathbf{X}-\mathbf{U}\theta \right)=0$$
解得
$$\theta ={\left({\mathbf{U}}^T\mathbf{U}\right)}^{-1}{\mathbf{U}}^T\mathbf{X}$$
$$\theta =\frac{1}{N}\sum _{n=1}^{N}x\left[n\right]$$

七、协方差

我们再把视线放回δ-σ 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数据转换器从入门到精通》系列笔记完成后找机会再聊！
晚安，迟哥。