杨 鹏

(1.核工业理化工程研究院，天津 300180；2.国防科技工业核材料技术创新中心，天津 300180)

小波理论被认为是对傅里叶分析的重大突破，已成为当今从应用数学到信号与图像处理等众多领域的研究热点。小波变换是由法国科学家Morlet于1980年在进行地震数据分析时提出的。由于小波分析方法良好的时频分辨能力，近年来在工程上获得了广泛应用。但在实际信号分析过程中发现，小波变换可能引起严重的频率混叠现象。针对以上问题，本文试图从原理上分析小波分析产生频率混叠的原因，并介绍小波分析消除混叠改进算法的Matlab实现过程[1]。

1 小波算法及其频率混叠

Mallat算法在小波分析中的地位相当于快速傅里叶变换(FFT)在傅里叶变换中的地位。Mallat算法又称塔式算法，它由小波滤波器H、G和h、g对信号进行分解和重构，算法如下。

分解算法为：

(1)

式中：t为离散时间序列号，t=1,2,3,…,N；f(t)为原始信号；j为层数，j=1,2,3,…,J，J=log2N；H和G为时域中小波分解滤波器的系数；Aj为信号f(t)在第j层的近似部分(即低频部分)的小波系数；Dj为信号f(t)在第j层的细节部分(即高频部分)的小波系数[2]。

重构算法为：

(2)

式中：j为分解的层数，若分解的最高层即分解的深度为J，则j=J-1,J-2,…,1,0；h和g为时域中小波重构滤波器的系数。

Mallat算法中包括3个关键环节：与正交镜像滤波器卷积、隔点采样及隔点插零。这些环节要求正交镜像滤波器必须具有理想的截止特性，即h和H为理想低通滤波器，g和G为理想高通滤波器。而Mallat算法所用的正交镜像滤波器实际上并非理想滤波器[3]，以Daubechies(dbN)小波为例，若信号的采样频率为800 Hz，N=5和N=20时的dbN小波的频域特性如图1所示[4]。由图1可看出，低通部分及带通部分在交界处有明显重叠，重叠的长度随N的增大而减少。

2 消除频率混叠的小波分解与重构算法及其Matlab实现

小波变换的频域定义相当于用1族带通滤波器对信号进行滤波，随着小波变换尺度的减小，滤波器的中心频率向高频率移动的同时，其带通宽度也随之增加。小波理论提供了包括傅里叶分析所采用的三角基函数以外的多种小波基函数，使小波分析具备多种情况下的分析能力，作为一种信号的时频分析方法，具有多分辨 分析的特点，它将信号分解成低频部分和高频部分，分解过程中，低频部分失去的信息由高频部分捕获[5]。在第j+1层的分解中，只分解第j层的低频部分，保留第j层的高频部分，逐层分解下去，可进行很深层次的分解。在机械设备故障诊断中，当机器发生故障时，由于机器各零部件的结构不同，输出信号在各频率波段中的表现也不同。在故障诊断的研究中，信号采集后，要对信号进行分析处理，提取并考察其特征参量的变化。

图1 N=5和N=20时dbN小波的频域特性Fig.1 Frequency of dbN with N=5 and 20

基于上述原因，亟需一种消除小波分解和重构过程中产生的频率混叠的算法[6]。改进算法的主要思想为：利用傅里叶变换和傅里叶逆变换来去掉多余的频率成分。基于这一思想的改进算法如图2所示。此方法的Matlab实现过程[7]如下。

图2 改进的小波分解与重构快速算法Fig.2 Improvement algorithm for resolvement and reconstitution of wavelet

2.1 消除取得近似部分Aj过程中频率混叠的方法

近似部分是能量的低阶频率部分，其混叠会导致低阶信号混乱，影响信号基础频率的输出，其消除方法如下：

2) 将FFT结果中频率大于1+fs/2j(fs为信号采样频率)部分谱值置零；

3) 对置零后的结果进行快速傅里叶逆变换；

4) 对快速傅里叶逆变换的结果进行隔点采样，采样后的结果作为Aj再进行下一步分解。

在Matlab中，此过程的实现通过定义一个函数m文件来实现，设函数名为r_appcoef。首先用[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters(wname)语句获得小波分解所要用到的正交镜像滤波器，其中Lo_D,Hi_D,Lo_R,Hi_R分别对应G、H、g、h；wfilters(wname)表示所选择的小波基函数名[8-9]。

利用fft(convn(A,Lo_D))命令即可实现步骤1。

值得探讨的是步骤2中的频率问题。从滤波的角度看，Mallat算法是将信号f(t)分解到一系列子带的滤波过程。各子带的频率范围与信号f(t)的采样频率有关，设原始信号的采样频率为f(t)，Mallat算法中各子带的频率范围列于表1。

表1 Mallat算法的频率范围Table 1 Limit’s frequency of Mallat algorithm

步骤2中，在得到FFT的结果后，首先应构建连续频率轴，然后再通过隔点插零的形式进行分解。其流程图示于图3。

图3 消除取得近似部分频率混叠的流程图Fig.3 Flow chart of reducing frequency mixed of approximate part

2.2 消除取得细节部分Dj过程中频率混叠的方法

细节部分是能量的高阶频率部分，其混叠会导致高频信号混乱，影响信号细节频率的输出，其消除方法如下：

2) 将FFT结果中f<1+fs/2j部分谱值置零；

3) 对置零后的结果进行快速傅里叶逆变换，结果作为DJ，不进行隔点采样。

消除流程示于图4。

图4 消除取得细节部分频率混叠的流程图Fig.4 Flow chart of reducing frequency mixed of detail part

2.3 消除近似部分Aj重构过程中频率混叠的方法

近似部分是能量的一个低阶重要部分，其重构过程中的混叠会导致低阶信号混乱，影响信号的输出，其消除方法如下：在第2j尺度上，由未被隔点采样的信号开始重构，这样在2j尺度上不进行隔点插零，而是从下1步，即第2j-1尺度上开始隔点插零，由于第1步重构不产生虚假频率成分，而后续的虚假成分一般可被h滤掉，因此就消除了Aj重构过程中的频率混叠[10]。

2.4 消除细节部分Aj重构过程中频率混叠的方法

细节部分是能量的一个高频重要部分，其混叠会导致高频信号混乱，影响信号细节频率的输出，其消除方法如下：

1) 设dj是由Dj重构的结果，对dj进行FFT；

2) 将FFT结果中非[1+fs/2j，fs/2j]部分谱值置零；

3) 对置零后的结果进行快速傅里叶逆变换，其结果即为真实的dj。

3 改进后的算法分析实例

在实现过程中，由理论上的推导转化成实际的编程，还应注意以下问题。

1) 信号在与滤波器函数卷积中的长度问题。在将信号与滤波器函数进行卷积时，使用的是Matlab中的convn函数，x=convn(a,b)计算的长度为size(a)+size(b)-1。在分解与重构过程中都有信号与滤波器卷积的环节，这样最后得到的结果其长度必然与原信号的长度不同，因此需利用Matlab中的wkeep函数来控制信号的长度。设N为原信号长度，A1为第1层近似部分重构后的结果，则运行A1=wkeep(A1,N)命令即可取得A1的长度为N的中心部分。这样就保证了信号长度的一致性。

2) 信号隔点插零和隔点采样的方式。Matlab提供了dyaddown和dyadup两个小波分析中用于隔点采样和隔点插零的函数。dyaddown的用法为Y=dyaddown(X,EVENODD)，当EVENODD取偶数时，Y(k)=X(2k)；当EVENODD取奇数时，Y(k)=X(2k+1)。dyadup的用法为Y=dyadup(X,EVENODD)，当EVENODD取偶数时，Y(2k-1)=X(k)，Y(2k)=0；当EVENODD取奇数时，Y(2k-1)=0，Y(2k)=X(k)。为保证信号长度在经隔点采样和隔点插零后不发生长度上的变化，应取适当的EVENODD值以保持一致。在使用dyaddown时EVENODD取0，在使用dyadup时EVENODD取1。

为更清楚地体现改进算法的效果，对某实验过程中的测量信号进行分析。待分析信号为该实验采集到的微弱突变信号，采样频率为500 Hz，采样点数取5 000，取dbN小波，N取40。利用默认小波分析算法和改进小波分析算法对信号进行3层分解，结果分别示于图5、6。由图5可看出，信号D2和D1、D3和D2的频谱存在明显混叠。而从图6可看出，各层的频率都在理论上的频段范围内，未出现频率混叠现象。

4 结论

本文设计了Mallat小波变换改进算法，利用傅里叶变换和傅里叶逆变换构成Mallat小波变换在分解和重构过程中所需的严格正交镜像滤波器。改进前后的小波算法对实际信号的分析结果表明，本文设计的Mallat小波变换改进算法可消除信号分析中出现的频率混叠现象。

图5 默认小波分析算法分析结果Fig.5 Analysis result of wavelet algorithm

图6 改进小波分析算法分析结果Fig.6 Analysis result of improvement wavelet algorithm