【现代信号处理】17 - 基于滤波器组的谱估计

基于滤波器组的谱估计–Filter-Bank Method

1. 滤波器组

  首先我们介绍一下，什么是滤波器组

1.1 分段处理谱估计问题

  通常随机信号是无限长的，连续的。但是实际信号的经过采样之后的，就变成了离散的，有限长度的。

  现在我们想用有限长的数据估计随机信号的功率谱

Z ( t ) ⇒ { Z ( k ) } k = 1 N − 1 Z(t) \Rightarrow \{ Z(k)\}_{k=1}^{N-1} Z(t)⇒{Z(k)}k=1N−1​

  我们实际是希望得到每一个频点上的功率谱的分量。但是我们只有有限长的数据，这是不可能实现的事情。

  我们只能降低要求。就是不求全部频率点的功率谱，而是划分一个小区间，就用有限的数据，求这一个小区间的功率谱的和

  现在问题就变成了。我们的数据经过采样之后，就划分出了一段工作频带，频带取决于采样率。采样频率的倒数，就是带宽。然后将频宽归一化到(-pi,pi)

  在这个频带内可以划分出很多的小段。只要划分出的小段的数量，比采样的数据数目少即可。

1.2 滤波器作为谱估计的手段

  为了能够获得这个小段，我们只需要让信号通过一个滤波器即可，滤波器的通带就是这个小段

  只要我们滤波器是理想的，我们就能够得到这一小段内功率的和

  因此滤波器就变成了解决谱分析问题的工具

1.3 滤波器与滤波器组

  为什么会有滤波器组的概念呢，因为我们一个滤波器只能得到一个频点周围的功率的和。

  如果我们想看很多频点周围的功率的和，我们就要让滤波器滑动起来，就有了滤波器组的概念

  我们得到了多个频点的功率和之后，就能够得到一个相对比较准确的谱估计了

2. 滤波器组与周期图之间的关系

  我们用滤波器组的观点来重新看待周期图

$${\hat{S}}_Z(\omega )=\frac{1}{N}|\sum _{k=1}^{N-1}Z(k)exp(-j\omega k){|}^2$$

  改变一下形式，在模里面加一个相位,因为相位改变对模没有影响，所以变换前后值不变

$${\hat{S}}_Z(\omega )=\frac{1}{N}|\sum _{k=1}^{N-1}Z(k)exp(j\omega (N-k)){|}^2$$
  这样做是为了构造卷积的样子

$$\begin{gathered} \text{Let }{h}_k=e(j\omega k) \\ k=0,1,...,N-1 \end{gathered}$$

$$\text{Then }{\hat{S}}_Z(\omega )=\frac{1}{N}|\sum _{k=0}^{N-1}Z(k)h_{N-k}(\omega ){|}^2$$

  再变形下，让式子更加接近卷积

$$\begin{gathered} =\frac{1}{N}|\sum _{k=0}^{\infty }Z(k)h_{N-k}(\omega ){|}^2 \\ {h}_k(\omega )=\{\begin{array}{ll}exp(j\omega k)& \text{k=1,...,N-1 }\\ 0& \text{others }\end{array} \end{gathered}$$

  然后我们计算下滤波器的频率响应

$$\begin{gathered} H({\omega }^{\prime })=\sum _{k=1}^{N-1}{h}_k({\omega }^{\prime })exp(-j{\omega }^{\prime }k) \\ =\sum _{k=1}^{N-1}exp(j(\omega -{\omega }^{\prime })k) \end{gathered}$$

$$|H({\omega }^{\prime })|=\frac{sin(\frac{N}{2}(\omega -{\omega }^{\prime }))}{sin(\frac{1}{2}(\omega -{\omega }^{\prime }))}$$

  上图是ω为0的时候的频率响应。如果要滑动的话，就是让这个图沿着ω轴滑动即可。我们从滤波器角度审视周期图的话，发现因为主瓣模糊和旁瓣泄漏，这个滤波器其实并不是非常好。

  我们可以发现，周期图就是在用一个滤波器在频率轴进行滑动，但是滤波器的质量却离理想的差很远。

  但是周期图的滤波器效果不好，只是因为我们没有对这个滤波器进行单独的设计。我们只是根据功率谱的定义式，然后通过离散化的方式得到了周期图。

  下面我们就要对通过对滤波器进行单独设计，然后得到更好的谱估计。

3. 数据无关的滤波器组方法–Slepian 滤波器

3.1 设计目标

  也就是我们要找到一组h_k，我们对他的频响有要求。我们先在零频上找这个滤波器，然后可以通过修改相位的方法进行滑动，找到滤波器组。

  首先，我们要求在整个工作频带内的通过能力是归一化的

{ h k } k = 0 N − 1 → H ( ω ) 1 2 π ∫ − π π ∣ H ( ω ) ∣ 2 d ω = 1 \{ h_k \}_{k=0}^{N-1} \rightarrow H(\omega) \\ \frac{1}{2\pi} \int_{-\pi}^{\pi} |H(\omega)|^2 d\omega = 1 {hk​}k=0N−1​→H(ω)2π1​∫−ππ​∣H(ω)∣2dω=1

  并且我们要求在我们关注的频带内的通过能力是最大的。因为这个滤波器是放在零频上的，因此要求零频附近的响应能力最大

$$ma{x}_H{\int }_{-\beta \pi }^{\beta \pi }|H(\omega ){|}^{2}d\omega$$

  相当于，我们做最优化，在全局范围内做了一个约束，然后在局部范围内做了目标函数

3.2 计算

3.2.1 约束条件的计算

$$\begin{gathered} H(\omega )=\sum _{k=0}^{N-1}{h}_{k}exp(-j\omega k)=a(\omega {)}^H\ast h \\ h=(h_0,...,h_{N-1}) \\ a(\omega )=(1,exp(j\omega ),...,exp(j(N-1)\omega ) \end{gathered}$$

  H是共轭转置的意思。h和a都是列向量。

$$\begin{gathered} \frac{1}{2\pi }{\int }_{-\pi }^{\pi }|H(\omega ){|}^{2}d\omega \\ =\frac{1}{2\pi }{\int }_{-\pi }^{\pi }(h^{H}a(\omega )a(\omega {)}^{H}h)d\omega \\ =h^H(\frac{1}{2\pi }{\int }_{-\pi }^{\pi }a(\omega )a(\omega {)}^{H}d\omega )h \end{gathered}$$

  被积函数是一列乘一行，是一个矩阵

$$\begin{gathered} \text{Let }A(\omega )=a(\omega )a^H(\omega ) \\ {A}_{kn}(\omega )=exp(j(k-n)\omega ) \end{gathered}$$

  同时

$${\int }_{-\pi }^{\pi }{A}_{kn}(\omega )\omega ={\delta }_{kn}$$

  因此可知

$$A(\omega )=I$$

  所以

$$\begin{gathered} \frac{1}{2\pi }{\int }_{-\pi }^{\pi }|H(\omega ){|}^{2}d\omega \\ =h^H(\frac{1}{2\pi }{\int }_{-\pi }^{\pi }a(\omega )a(\omega {)}^{H}d\omega )h \\ =h^H(\frac{1}{2\pi }{\int }_{-\pi }^{\pi }A(\omega )d\omega )h \\ =h^{H}h=||h|{|}^2=1 \end{gathered}$$

3.2.2 优化目标的计算

$$\frac{1}{2\pi }{\int }_{-\beta \pi }^{\beta \pi }|H(\omega ){|}^{2}d\omega$$

  β限制的是滤波器的通带宽度。β应该是个很小的数，但是不能无限的小下去。2π的区间里面分段数不能够大于N

$$\begin{gathered} \beta <<1 \\ \frac{2\pi }{2\beta \pi }\le N \\ \Rightarrow \beta \ge \frac{1}{N} \end{gathered}$$

  来看约束条件

$$\begin{gathered} \frac{1}{2\pi }{\int }_{-\beta \pi }^{\beta \pi }|H(\omega ){|}^{2}d\omega \\ =h^H(\frac{1}{2\pi }{\int }_{-\beta \pi }^{\beta \pi }A(\omega )d\omega )h \end{gathered}$$
  我们假定

$$\text{Let }\Gamma =\frac{1}{2\pi }{\int }_{-\beta \pi }^{\beta \pi }A(\omega )d\omega$$

$$\begin{gathered} {\Gamma }_{kn}=\frac{1}{2\pi }{\int }_{-\beta \pi }^{\beta \pi }{A}_{kn}(\omega )d\omega \\ =\frac{1}{2\pi }{\int }_{-\beta \pi }^{\beta \pi }exp(j(k-n)\omega )d\omega \end{gathered}$$

  因为积分区间是对称的，exp中虚部的sin是个奇函数，因此在积分中是0，这个积分只会剩下cos部分

  因此

$$\begin{gathered} {\Gamma }_{kn}=\frac{1}{2\pi }{\int }_{-\beta \pi }^{\beta \pi }cos((k-n)\omega )d\omega \\ =\frac{1}{2\pi }\frac{sin((k-n)\omega )}{k-n}{|}_{-\beta \pi }^{\beta \pi } \\ =\frac{sin((k-n)\beta \pi )}{(k-n)\pi } \end{gathered}$$

  分析优化条件

$$ma{x}_h{h}^H\Gamma hs.t.h^{H}h=1$$
  我们用拉格朗日数乘法进行求解

$$L(h,\lambda )=h^H\Gamma h+\lambda (h^{H}h-1)$$

  注意，γ是个对称矩阵

$${\Gamma }^H=\Gamma$$

  求导

$$\nabla L(h,\lambda )=2\Gamma h-2\lambda h=0$$

  可得

$$\Gamma h=\lambda h$$

  说明h是γ矩阵的特征向量，是特征值最大的那个的特征矢量

3.2.3 对γ矩阵的分析

$$\begin{gathered} \Gamma =\frac{1}{2\pi }{\int }_{-\beta \pi }^{\beta \pi }A(\omega )d\omega \\ \approx \frac{1}{2\pi }\sum _{k=-L}^{L}a(\frac{2\pi }{N}k)a^H(\frac{2\pi }{N}k)\ast \frac{2\pi }{N} \end{gathered}$$

  其中

$$L=\frac{\beta \pi }{2\pi /N}=\frac{\beta N}{2}$$

  这个近似如果要做的好的话，βN要充分大，也就是通带β要尽可能宽，但是这样就会损失分辨率

$$\begin{gathered} \beta N>>1 \\ \beta >>\frac{1}{N} \end{gathered}$$

  因此这个近似可能不准

$$\begin{gathered} \Gamma \approx \frac{1}{2\pi }\sum _{k=-L}^{L}a(\frac{2\pi }{N}k)a^H(\frac{2\pi }{N}k)\ast \frac{2\pi }{N} \\ =\frac{1}{N}\sum _{k=-L}^{L}a(\frac{2\pi }{N}k)a^H(\frac{2\pi }{N}k)=\tilde{\Gamma } \end{gathered}$$

  我们发现这是一个2L+1个矩阵相加得到的式子。因为每个矩阵都是一列乘一行，是个秩一矩阵。2L+1个秩一矩阵相加，秩最多就是2L+1

$$r(A+B)\le r(A)+r(B)$$

$$a(\frac{2\pi k}{N})=(1,exp(j\frac{2k\pi }{N}),...,exp(j\frac{N-1}{N}2k\pi ))$$
  我们来看看a的内积

$$a^H(\frac{2\pi k_1}{N})a(\frac{2\pi k_2}{N})=\{\begin{array}{ll}0& k_1\cancel{=}{k}_2\\ N& k_1=k_2\end{array}$$

  我们发现这是个正交的向量，并且，如果我们把N放到a里面，这还是个规范的向量

$$\frac{1}{\sqrt{N}}{a}^H(\frac{2\pi k_1}{N})\frac{1}{\sqrt{N}}a(\frac{2\pi k_2}{N})=\{\begin{array}{ll}0& k_1\cancel{=}{k}_2\\ 1& k_1=k_2\end{array}$$

  继续分析γ矩阵

$$\tilde{\Gamma }=\frac{1}{N}\sum _{k=-L}^{L}a(\frac{2\pi }{N}k)a^H(\frac{2\pi }{N}k)$$

  假设有这样的k
$$k\in [-L,L]$$

$$\begin{gathered} \tilde{\Gamma }\ast a(\frac{2\pi k}{N})=(\frac{1}{N}\sum _{l=-L}^{L}a(\frac{2\pi }{N}l)a^H(\frac{2\pi }{N}l))\ast a(\frac{2\pi k}{N}) \\ =\frac{1}{N}\sum _{l=-L}^{L}a(\frac{2\pi }{N}l)(a^H(\frac{2\pi }{N}l)\ast a(\frac{2\pi k}{N})) \\ =\frac{1}{N}\sum _{l=-L}^{L}a(\frac{2\pi }{N}l)N\ast {\delta }_{kl}=a(\frac{2\pi }{N}k) \end{gathered}$$

  可得

$$\tilde{\Gamma }\ast a(\frac{2\pi k}{N})=a(\frac{2\pi }{N}k)$$

  我们就得到了λ为1时候的特征矢量。并且这样的特征矢量有2L+1个，因此1是2L+1重的特征值。并且这些特征矢量彼此正交

$$a(\frac{2\pi k}{N})$$

  这里近似的有些厉害了，因为得到的特征矢量是复三角函数，但是实际上并不是

  γ矩阵的特征值也有特点，前面的一些非常大，后面的一些几乎为0，这是因为γ矩阵的秩最大为2l+1

3.2.4 slepian window

  这样设计得到的滤波器叫做slepian滤波器

4. Multitaper

4.1 周期图谱中的trade off

  在周期图分析中，我们提到过Bias Variance Trade off，就是我们扩大数据量之后，对方差校正没有影响，对均值校正有影响，为了降低方差，我们只能把数据分成不同的段，求功率谱，然后取平均，这样才能得到方差比较小的功率谱估计，可是这样的话，得到的数据均值偏差就会变大，就有了Bias Variance Trade off

4.2 slepian window与周期图谱估计

4.2.1 Multitaper

  我们可以通过多窗口的方法来解决周期图中的trade-off问题。这种方法叫做multiple windows(Multitaper)

  也就是，我们对同一组数据，使用不同的窗进行谱分析，因为窗是我们特别设计的，所以得到的谱分析结果是不相关的。然后谱分析结果取均值，就能够降低方差。这样的话，就能够同时优化均值和方差了

  怎么设计这样的窗呢?slepian window就有这样的功能

  我们就用γ矩阵产生的不同滤波器来做这个思路的谱估计

$$\begin{gathered} {\Gamma }_{kn}=\frac{1}{2\pi }{\int }_{-\beta \pi }^{\beta \pi }cos((k-n)\omega )d\omega \\ =\frac{1}{2\pi }\frac{sin((k-n)\omega )}{k-n}{|}_{-\beta \pi }^{\beta \pi } \\ =\frac{sin((k-n)\beta \pi )}{(k-n)\pi } \end{gathered}$$

  假设γ矩阵有若干个特征值，然后我们就有了新的谱估计结果

$$\begin{gathered} \Gamma h^{(k)}={\lambda }_k{h}^{(k)} \\ k=1,2,...,L \\ \Rightarrow {\hat{S}}_Z=\frac{1}{NL}\sum _{k=1}^L|\sum _{n=0}^{N-1}{h}_n^{(k)}Z(n)exp(-jn\omega ){|}^2 \end{gathered}$$

  这个谱估计，里面是一个周期图，外面是不同滤波器得到结果的均值。周期图里面是加了窗的，就是L个大特征值对应的特征矢量

4.2.2 证明不同滤波器谱分析结果不相关

  因为我们进行Multitaper的谱分析时，必须要证明每个滤波器得到的功率谱估计是不相关的，所以我们在这里证明一下

$$\begin{gathered} E(\sum _{n=0}^{N-1}{h}_n^{(k_1)}Z(n))(\overline{\sum _{l=0}^{N-1}{h}_l^{(k_2)}Z(l)}) \\ =\sum _{n=0}^{N-1}\sum _{l=0}^{N-1}{h}_n^{(k_1)}\overline{h_l^{(k_2)}}E(Z(n)Z^{\ast }(l)) \\ =\sum _{n=0}^{N-1}\sum _{l=0}^{N-1}{h}_n^{(k_1)}\overline{h_l^{(k_2)}}{R}_Z(n-l) \\ =\sum _{n=0}^{N-1}\sum _{l=0}^{N-1}{h}_n^{(k_1)}\overline{h_l^{(k_2)}}\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{S}_Z(\omega )exp(j\omega (n-l))d\omega \\ =\frac{1}{2\pi }{\int }_{-\pi }^{\pi }(\sum _{n=0}^{N-1}{h}_n^{(k_1)}exp(j\omega n))(\sum _{l=0}^{N-1}\overline{h_l^{(k_2)}}exp(-j\omega l))S_Z(\omega )d\omega \end{gathered}$$

  我们做一个换元，令-ω = ω，为了凑成傅里叶变换。后面的微元变号，前面的积分换方向会抵消。功率谱密度是平方，正负号不影响。然后就exp里面改变符号即可

$$\begin{gathered} =\frac{1}{2\pi }{\int }_{-\pi }^{\pi }(\sum _{n=0}^{N-1}{h}_n^{(k_1)}exp(-j\omega n))(\sum _{l=0}^{N-1}\overline{h_l^{(k_2)}}exp(j\omega l))S_Z(\omega )d\omega \\ =\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{H}^{(k_1)}(\omega )\overline{H^{(k_2)}(\omega )}{S}_Z(\omega )d\omega \end{gathered}$$
  因为对功率谱密度加了两个滤波器，两个滤波器都是在零频附近响应比较强，因此，上式可以近似等于在零频附近的能量。

  然后我们假设功率谱密度几乎在零频附近是个常数
$$\begin{gathered} \approx \frac{1}{2\pi }{\int }_{-\beta \pi }^{\beta \pi }{H}^{(k_1)}(\omega )\overline{H^{(k_2)}(\omega )}{S}_Z(\omega )d\omega \\ \approx S_Z(0)\frac{1}{2\pi }{\int }_{-\beta \pi }^{\beta \pi }{H}^{(k_1)}(\omega )\overline{H^{(k_2)}(\omega )}d\omega \\ =S_Z(0)\frac{1}{2\pi }{\int }_{-\beta \pi }^{\beta \pi }(a^H(\omega )h^{(k_1)})(\overline{(h^{(k_2)}{)}^T\overline{a(\omega )}})d\omega \\ =\frac{S_Z(0)}{2\pi }(h^{(k_2)}{)}^H({\int }_{-\beta \pi }^{\beta \pi }a(\omega )a^H(\omega )d\omega )h^{(k_1)} \\ =\frac{S_Z(0)}{2\pi }(h^{(k_2)}{)}^H\Gamma h^{(k_1)}=0(k_1\cancel{=}{k}_2) \end{gathered}$$

4.2.3 小结

  其实周期图谱估计也是一种滤波器组方法，只不过周期图中用的滤波器组没有进行优化设计，有主瓣模糊和旁瓣泄漏的问题

  而对滤波器通过约束条件和优化目标的设计，就得到了slepian滤波器组方法

  同时，slepian除了最大的特征值对应的特征向量可用来做滤波器组以外，slepian的其他特征向量也可以用来做滤波器组谱估计，虽然性能差一点。

  但是，值得注意的是，这些滤波器组处理过数据之后，得到的功率谱数据是不相关的。因此，对多个滤波器组得到的谱估计结果进行平均，就可以得到方差更小的功率谱估计结果。

  这种multitaper的方法解决了周期图谱估计中矛盾点。即扩大数据长度对方差的误差减小没有帮助，把数据分段处理，最后加权平均又会影响均值的误差。而使用slepian滤波器组就不用对数据分段，利用不同滤波器组处理同一组数据，然后加权平均，同时兼顾了均值和方差的误差的减小。

5. 数据相关的谱估计方法-Capon Spectrum

5.1 问题引入

  slepian在设计滤波器的时候，是数据独立的。因为约束条件和优化条件都是针对h来的，跟数据域没有关系，因此通过slepian做的滤波器组方法是数据无关的。

  但是，我们希望从数据相关的角度来进行谱估计方法设计。比如机器学习，训练的模型都是数据相关的，用猫的图片进行训练，得到的模型肯定识别不出汽车

  接下来引入一种数据相关的谱分析方法

5.2 约束条件与目标函数

  假设我们有一个滤波器，并且要求这个滤波器在给定频率上的响应是个常数，这里就取1

{ h k } k = 0 N − 1 Given  ω 0 H ( ω 0 ) = 1 \{ h_k\}^{N-1}_{k=0} \\ \text{Given } \omega_0 \\ H(\omega_0) = 1 {hk​}k=0N−1​Given ω0​H(ω0​)=1

  我们在这个约束的条件下，要求这个滤波器输出的功率最小

$$\begin{gathered} mi{n}_{h}E|h^{H}Z{|}^2 \\ Z=(Z(0),...,Z(n-1){)}^T \\ h=(h0,...,h_{N-1}{)}^T \end{gathered}$$

  相当于我们要求了滤波器在给定频率点上的响应，并且要求这个滤波器在这个频点周围的面积尽可能窄，这样才能有最小的功率输出

$$mi{n}_{h}E|h^{H}Z{|}^{2}s.t.a^H(\omega )h=1$$

5.3 与slepian的对比

  slepian设计滤波器的时候，是要求滤波总的功率谱是一定的，要求经过滤波器在给定频点处的功率最大。说明滤波器在给定频点的地方占据了最大的空间。是约束好整体条件，优化局部。

$$\begin{gathered} \text{Slepian} \\ ma{x}_h{h}^H\Gamma hs.t.h^{H}h=1 \\ \text{1964 Bell} \end{gathered}$$

  而这种方法是约束好了局部条件，然后优化全局。就是要求在频点处的响应是确定的，并且要求功率输出最小。
$$\begin{gathered} mi{n}_{h}E|h^{H}Z{|}^{2}s.t.a^H(\omega )h=1 \\ \text{1969 MIT} \end{gathered}$$

  这种方法在优化里面叫对偶

  MIT的方法在优化中已经用到了数据域，因此是数据相关的。并且约束条件上，还用到了频率，还是频率相关的

  因此这种方法设计的滤波器组，并不是简单的数据滑动，而是在每一个频点上，都是专门设计的滤波器。而slepian是同一个滤波器在不同频点上滑动

  因此MIT的这种方法在所有非参数谱分析中，是最好的

5.4 计算

$$\begin{gathered} E|h^{H}Z{|}^2=E(h^{H}Z{Z}^{H}h) \\ =h^{H}E(Z{Z}^H)h=h^H{R}_{Z}h \end{gathered}$$

  约束条件和优化目标

$$mi{n}_{h}E|h^{H}Z{|}^{2}s.t.a^H(\omega )h=1$$

  拉格朗日数乘法

$$\begin{gathered} L(h,\lambda )=\frac{1}{2}{h}^H{R}_{Z}h-\lambda (a^H(\omega )h-1) \\ {\nabla }_{h}L=R_{Z}h-\lambda a(\omega ) \\ h=\lambda R_Z^{-1}a(\omega ) \end{gathered}$$

  下面我们就求λ

  代入约束条件
$$a^H(\omega )h=1$$

$$\begin{gathered} a^H(\omega )\lambda R_Z^{-1}a(\omega )=1 \\ \lambda =\frac{1}{a^H(\omega )R_Z^{-1}a(\omega )} \end{gathered}$$

  我们就得到了确定频点处的滤波器

$$h=\frac{R_Z^{-1}a({\omega }_0)}{a^H({\omega }_0)R_Z^{-1}a({\omega }_0)}$$

  我们就可以得到新的功率谱估计了

  这个功率谱就是滤波器输出的功率

$$\begin{gathered} S_Z({\omega }_0)=h^H{R}_{Z}h \\ =\frac{a({\omega }_0{)}^H{R}_Z^{-1}{R}_Z{R}_Z^{-1}a({\omega }_0)}{(a^H({\omega }_0)R_Z^{-1}a({\omega }_0){)}^2} \\ =\frac{a({\omega }_0{)}^H{R}_Z^{-1}a({\omega }_0)}{(a^H({\omega }_0)R_Z^{-1}a({\omega }_0){)}^2} \\ =\frac{1}{a^H({\omega }_0)R_Z^{-1}a({\omega }_0)} \end{gathered}$$

5.5 小结

  这种谱叫做Capon Spectrum。在非参数功率谱估计中做到了极致

  这种方法也叫做MVDR(Minimum Variance Distortionless Response) 最小方差无损响应

  Capon的方法的滤波器谱在频点处是非常尖的

   取法乎上，得之中。就是对capon最好的描述。因为capon把别人追求的目标函数变成了约束条件