武汉理工大学毕业设计（论文）
到的均值趋于零为止，这样就得到了第1个IMF分量c1(t)，它代表信号s(t)中最高频率的分量：
h1(k?1)(t)?m1k(t)?h2k(t) （10）
c1(t)?h1k(t) （11）
（4）将c1(t)从s(t)中分离出来，即得到一个去掉高频分量的差值信号r1(t)，即有
r1(t)?c2(t)?r2(t)
（12）
将r1(t)作为原始数据，重复步骤（1）、（2）（3），得到第二个IMF分量c2(t)，重复n次，得到n个IMF分量。这样就有：
r1(t)?c2(t)?r2(t)???（13） ? rn?1(t)?cn(t)?rn(t)??2.5 Hilbert谱和边际谱
在IMF定义和EMD的基础上，Huang等人系统地提出了一种分析信号的新理论或新方法，命名为希尔伯特一黄变换（Hilbert-Huang Transform简称HHT）它包括两个大组成部分，EMD和与之相应的Hilbeer谱分析方法。即首先用EMD将任意信号s(t)分解成有限个IMF的和。
ns(t)??cj?1j(t)?rn(t) （14）
然后分别对每一个IMF分量用Hilbert变换进行谱分析。最后得到信号的瞬时频率表示
ns(t)?Re?at(t)ei?1j?1(t)j?1(t)dt?Re?ai(t)e?i?1n （15）
这里省略了残余函数rn(t)，Re表示取实部。称式（12）右边为Hilbert时频谱，简称Hilbert谱，记作
nH(?,t)?Re?i?1j?1(t)dtai(t)e? （16）
8
武汉理工大学毕业设计（论文）
它是瞬时振幅在频率．时间平面上的分布。
在式（12）中省略残余函数rn(t)是因为它或者是一个常数，或者是一个单调函数。虽然可以把rn(t)看作一个长周期波的一部分，但考虑到长周期的不确定性，及信号所包含的信息主要在高频分量中，因此做了省略处理。
展开式（12）中，每个分量的幅值和相位都是随时间可变的，而同样信号s(t)的Fourier变换展开式为
?s(t)?Re?aeii?1j?1t （17）
其中ai,?i为常数。这清楚地表明：HHT对信号的瞬时频率表示是Fourier展开的一般化，它不仅提高了信号的效率，而且能够表示可变的频率。因此，新方法突破了傅立叶变换的束缚。
用Hilbert谱可以进一步定义边际谱为
H(?)????? H(?,t)dt （18）
这里由HHT得到的边际谱与Fourier频谱有相似之处，从统计观点上来看，它表示了该频率上振幅（能量）在时间上的累加，能够反映各频率上的能量分布，但因为瞬时频率定义为时间的函数，不同以往Fourier等需要完整的振荡波周期来定义局部的频率值，而且求取的能量值不是全局定义的。因此对信号的局部特征反映更准确，在这方面优于Fourier谱。尤其是在分析非平稳信号时，这种定义对于频率随时间随时变化的信号特征来说，能够反映真实地振动特点。
为了直观的说明Hilbert谱和边际谱，下面给出一个仿真信号来加以说明： 给出的信号是正弦叠加信号，分别是y1?3*sin(2*pi*280*t)和
y2?5*sin(2*pi*93*t)，通过EMD分解程序包emd、hhspectrum和toimage可
以得到信号的Hilbert谱，分析给定的信号我们可以知道其频率分布在280HZ和93HZ处，经过HHT我们观察其Hilbert谱，可以得出同样的结论即280HZ和93HZ处分布着信号的主要能量。
>> clear; >> fs=1000; >> tspan=2; >> t=1/fs:1/fs:tspan; >> N=length（t）;
>> y1=3*sin（2*pi*280*t）;
9
武汉理工大学毕业设计（论文）
>> y2=5*sin（2*pi*93*t）; >> y=y1+y2; >> imf=emd（y）;
>> [A, fa, tt] = hhspectrum（imf）; >> if size（imf,1） > 1
[A,fa,tt] = hhspectrum（imf（1:end-1, :））; else
[A,fa,tt] = hhspectrum（imf）; end
>> [E, tt1] = toimage（A,fa,tt,length（tt））; >> for k = 1:size（E,1）
bjp（k） = sum（E（k,:））*1/fs*1/0.05; end
>> f = （0:N-3）/N*（fs/2）; >> figure（1） >> plot（y）; >> figure（2）
>> imagesc（tt1,[0,0.5*fs],E）; >> set（gca,'YDir','normal'） >> colormap（flipud（gray）） >> xlabel（'时间 T/s'）; >> ylabel（'频率 f/HZ'）;
10
武汉理工大学毕业设计（论文）
50045040035030025020015010050020040060080010001200时间 T/s140016001800频率 f/HZ
图2.1 仿真信号的Hilbert谱图
上例中给出的是仿真信号的Hilbert谱，由于Hilbert谱具有良好的时频聚集性，通过分析信号的Hilbert谱我们可以方便的看出信号的频率变化分布情况。下面给出的信号同样是正弦叠加信号，分别是y1?3*sin(2*pi*241*t)和
y2?5*sin(2*pi*73*t)，通过EMD分解程序包hhspectrum和toimage我们按
照以下程序代码可以得到信号的边际谱，如上图2.2所示。
>> fs=1000; >> tspan=2; >> t=1/fs:1/fs:tspan; >> N=length（t）;
>> y1=3*sin（3*pi*250*t）; >> y2=5*sin（3*pi*95*t）; >> y=[y1;y2;zeros（size（y1））]; >> [A,fa,tt]=hhspectrum（y）;
>> [E,tt1]=toimage（A,fa,tt,length（tt））; >> E=flipud（E）; >> for k=1:size（E,1）
11
武汉理工大学毕业设计（论文）
bjp（k）=sum（E（k,:））*1/fs*1/tspan; end
>> f=（0:N-3）/N*（fs/2）; >> plot（f,bjp）; >> xlabel（'频率 / Hz'）; >> ylabel（'幅值'）;
54.543.53幅值2.521.510.50050100150200250300频率 / Hz350400450500
图2.2 仿真信号的Hilbert边际谱图
注意在程序中E=flipud（E）语句是实现图谱的翻转，这里flipud是一个使矩阵上下翻转的函数，通过翻转，可以使得到的边际谱的频率按从左到右递增。
我们分析信号的边际谱，即为得到个频率分量的能量分布，从而为我们对信号的分析提供帮助。
2.6 HHT的特点
与传统的信号或数据处理方法相比，HHT具有如下特点： （1） HHT能分析非线性非平稳信号。
传统的数据处理方法，如傅立叶变换只能处理线性非平稳的信号，小波变换虽然在理论上能处理非线性非平稳信号，但在实际算法实现中却只能处理线性非
12