详解抑制有限长时间序列中噪声的有效方法

weixin

　　抑制噪声的经典方法是用线性滤波器，使信号通过滤波器后不需要的频率成分被削弱。模拟滤波器是电感电容或有源器件构成的电路网络单元。如果需要用数字信号处理器代替，首先用巴特沃斯，切比雪夫或椭圆函数设计逼近所需的模拟滤波器频率特性，然后再数字化，并考虑实现结构和数据有限字长的影响等。这占据了数字信号处理教科书相当一部分内容。

　　但是，在实际工作中，尤其是基于计算机的离线数据分析，常是没有必要的，还可能费力不讨好。数据时间长度较小时，除非有特别明确的依据和需要，基于经典“线性卷积”滤波的思想和方法是不可取的。

　　有限长时间序列x(n),n=0,1,2,...,N-1，是N维空间中的一点。通过适当的变换处理使噪声与其它有用信号成分可以分离或突显出来，这样就可以抛去或得以修正，然后再重建其它成分构成的信号。或者根据先验知识，对信号建模，用观测数据来估计模型参数，再重构信号。这些是很好的思想方法。

　　以复指数序列exp(j2πkn/N)除以N的平方根(k和n=0，1，2，…,N-1)为第k+1行第n+1列元素的N x N矩阵，其行向量或列向量是N维空间的规范正交基。时间序列x(n)，可以唯一的表达为这组基向量的线性组合。组合系数由信号序列与基向量做内积得到，这即是所谓的FFT。在数字信号处理中，FFT常特指蝶形快速算法，N为2的整数次幂。但在这里，N可以是任意整数，强调在规范正交基上展开和应用，忽略实现傅里叶变换的具体计算机程序结构。

　　连续时间周期函数，可以用傅里叶级数展开，即时域周期化对应频率域离散化。香农采样理论表明，时间离散化时信号频谱周期延拓。FFT则意味着信号在时域和频域的表达都是离散化的，也是周期化的。FFT结果是无穷连续时间信号的频谱周期化后的采样逼近，有明确的物理意义。这可以指导对FFT结果的筛选或修正，以达到抑制噪声的目的。

　　图1 FFT域修正系数降噪的例子

　　图1是FFT域修正系数降噪的例子。蓝线，是动态磁敏感对比核磁共振成像(DSC-MRI)的原始数据。红线，是由幅度最大的8项(若包括直流和对称部分共17项)FFT系数重构的序列。绿线，则是约在数字频段0.078~0.156用高斯函数序列加权从1逐渐减到0，0.078以下保持完全不变，而其余FFT系数全置为0时，逆FFT所重建的信号。两者都很好地抑制了噪声。

　　注意谨慎运用数据补零延拓的方法，虽然很多教科书讲卷积和相关时都用补零后再做周期延拓。如果没有明确理由排斥，使用直接周期化延拓，概念更清晰，应用更方便，效果更佳。对图1中蓝线那样的数据，如果简单补零的话，吉布斯伪影更难对付，结果可能很糟，而应该用均值延拓或者先去除均值后再补零。

　　如上所述，信号可以用非常有限的几个幅度与时间无关的复指数序列的加权和去逼近。那么一般复指数(exp[(a+jb)n]，a,b为实数)序列加权和模型，其模板序列的幅度也是指数函数，逼近效果可以更佳，可以包含FFT模型，但不必涉及周期化离散化和补零延拓等问题。

　　这种复指数和模型，得到了很好的研究，其矩阵代数描述分析显得相当别致。在雷达目标极点估计和核磁共振波谱分析等中都有应用。数值计算，使用矩阵奇异值分解，特征值和伪逆等是矩阵代数教科书的标准内容，有通用程序，Matlab和Scilab也都有简单命令调用函数实现计算。

　　用状态空间算法(《北京理工大学学报》1996年 06期《噪声信号极点估计的一般方法》)估计6个主要复指数，然后用这几个复指数序列的加权和去拟合DSC-MRI数据，得到的降噪结果如图2中红线所示。效果很好。居士偏爱这一方法。

　　图2 6个复指数序列加权和拟合降噪的例子

　　FFT域修正方法和一般复指数和模型逼近方法，都利用了信号在频率域是局域化的的特点。频率域描述是冗余的，冗余部分主要是噪声。如果序列中包含一些时域高度局域化成分，其频率谱很宽，那么降噪效果变差。普通线性滤波器也是基于噪声或干扰在频域可分离设计的，不适用。这时有一个很好的选择是小波包方法(赛特居士新浪博客2011-02-19《八卦图以及小波变换与高维空间的最佳坐标系》)。

　　小波包的多分辨率特点可以适应时域或频域的不同局域化特点，匹配检测出信号中的有用信息。小波分析理论，提供了大量的正交和双正交滤波器序列，使人们可以构造N维空间的各种小波包基序列(《北京理工大学学报》1998-18-1《离散正交小波包及其在噪声抑制中的应用》)。

　　图3 Coiflet小波包处理DSC-MRI数据

　　用长度为30的Coiflet序列构造周期化小波包基序列，对DSC-MRI数据做小波包分解，并根据最小熵准则选择最佳基，最后用最佳基中最大的10项系数重建信号，结果如图3中红线所示。

　　图4 Coiflet小波包处理正弦序列

　　图4包括两张曲线图，是用长度为30的Coiflet序列处理正弦信号的例子，蓝线是已知不含噪声的信号，绿线是噪声污染后的数据，红线是小波包域筛选系数重建的结果，仅用了模值大于3倍噪声标准差的系数。重建信号相对于无污染信号是失真的，失真部分是相关随机序列，随所加白噪声样本的不同而变化。这一特殊信号最适宜使用前面的非小波包方法，失真会减小，但是用小波包依然可以有效地衰减噪声。

　　图5 双正交小波包处理DSC-MRI数据

　　图5显示使用双正交小波包处理DSC-MRI数据的结果。使用的滤波器对为：[0;0;-0.05;0.25;0.6;0.25;-0.05;0;0;0]; [0;-3/280;-3/56;73/280;17/28;73/280;-3/56;-3/280;0;0]。可见正交和双正交处理都是有效的。

　　做“小波变换”，可以先构造并存储变换矩阵和逆变换矩阵。正变换时做矩阵乘法，即与基做内积。在早期Matlab中，用循环语句做普通Mallat卷积下采样的话，速度反而慢很多，处理批量数据甚至可能是灾难。如果不用周期化，而用“无穷时间”零延拓假设，程序更复杂，计算速度问题更明显，而且0时间点和边界效应的控制容易发生错误。

　　有些小波文献还强调预滤波，这里更没必要。着眼于有限时间序列，使用直接周期化，强调构造N维空间的基序列而不是连续时间基函数，是明智的。“小波包”的处理类似。后来Matlab有了小波和小波包工具箱，但是还没机会使用和比较与居士计算结果的差别。这里的例子，是用居士的旧Matlab程序改写到Scilab后计算的结果。

　　居士也写了非周期正交和双正交的小波包程序，没有体会到任何益处。如果保留最佳基上的所有系数，用图3，4和5相应基重建的能量相对误差都小于1.3e-15，完全重建表明程序方案是可行的。降噪效果也不错。

　　在完全重建前提下，程序细节以及滤波器序列的选择，反转和调制等都可根据实际处理效果来调整，这里不必遵从某种标准程序实现或数值结果相同，但是，不失一般性，如果使用Matlab工具箱一类“标准计算”或更佳程序实现降噪，也是有效的或许更好。

　　本文来源于新浪Scite_Jushi的博客，作者：新浪赛特居士SciteJushi。