核心：以下代码计算HHT边际谱及其对应频率
工具包要求：G-Rilling EMD Toolbox，TFTB Toolbox
附：黄锷先生课题组开发的工具包（可以在 这里 找到），这里并未用到。

% Empirical mode decomposition, resulting in intrinc mode functions.
% Without parameter \'MAXMODES\' the process may be seriously delayed by
% decompose original signals into too many IMFs (not necessary, 9 is 
% enough generally)
imfs = emd(oriSig, \'MAXMODES\', 9);
% HHT spectrum: hhtS
[A, f, t] = hhspectrum(imfs);
[hhtS, ~, fCent] = toimage(A, f, t);
% Marginal hilbert spectrum: hhtMS, xf: correspondig frequency
for k = 1 : size(hhtS, 1)
    hhtMS(k) = sum(hhtS(k,  : )) * 1 / fs;
end
xf = fCent(1, :) .* fs;

G-Rilling EMD Toolbox工具包相关配置：了凡春秋

简单来说，设置好路径之后输入 install_emd 即可。

Matlab关于EMD分解后希尔伯特谱分析（相关函数均隶属 G-Rilling EMD Toolbox）

• hhspectrum 函数说明（8楼：老老的学生）

% [A,f,tt] = HHSPECTRUM(imf,t,l,aff)
% Input:
%- imf : matrix with one IMF per row    % 将emd分解得到的IMF代入就可以，就是你的程序中写的c变量，不用加最后一行的趋势项
%- t   : time instants    % 瞬时时间或持续时间 ？？（写[1:信号长度]就可以，真实的时间可以根据采样率转换
%- l   : estimation parameter for instfreq    % 瞬时频率的估计参数 ？？（写1就可以，决定瞬时频率估计时的边界从第几个点开始
%- aff : if 1, displays the computation evolution    % 显示计算进程选项，不想显示写0就可以
%
% Output:
%   - A   : amplitudes of IMF rows
%   - f   : instantaneous frequencies
%   - tt  : truncated time instants    % 截止时间 ？？（截断时间，返回的是瞬时频率对应的时间，要比原来信号的时间按短，由前面的l值决定）
%
% calls:
%   - hilbert  : computes the analytic signal
%   - instfreq : computes the instantaneous frequency    % 瞬时频率
%
% Example：
[A,f,tt] = hhspectrum(c(1:end-1,:),[1:n],1,0);
[im,tt,ff] = toimage(A,f,tt,512);
disp_hhs(im,tt,[],fs);

• 9楼：老老的学生

关于时频图的概念，我认为是与诸如小波，gabor等联合时频分析方法联系在一起的。

小波，gabor等具有多尺度分析的概念，得到的时频分布是一个二维的矩阵（横轴时间，纵轴频率，可以用不同的颜色（光谱图spectrogram）或瀑布图表示不同的幅度）。

对HHT来讲，我觉得气时频图的概念是有些不同的。EMD分解的作用就是把复杂的信号分界为简单的单分量的信号，使其可以应用瞬时频率的概念，hilbert变换的目的就是分析出瞬时频率。所以HHT在每一时刻得到的只有一个值，而不是像小波之类的得到一系列的值（多尺度分析）。所以我们从其时频分布图上看到的是一条线，而不是一幅图。

• 27楼：songzy41
  HHT的重要意义是通过EMD分解后得到IMP，用IMP做hilbert，对每个固有模态函数才有瞬时信息。原信号是有许多固有模态函数之和，所有原信号没有瞬时信息可言。这是HHT的基本出发点。

基于G-Rilling EMD Toolbox的Hilbert谱计算示例代码：nike815

%示例程序
N=1000;
n=1:N;
fs=1000;
t=n/fs;
fx=50;
fy=150;
x=2*sin(2*pi*fx*t);
y=sin(2*pi*fy*t);
z=x+y;
data=z;
imf=emd(data);  % 对输入信号进行EMD分解    
[A,f,t]=hhspectrum(imf);    % 对IMF分量求取瞬时频率与振幅：A：是每个IMF的振幅向量,f:每个IMF对应的瞬时频率，t：时间序列号
[E,t,Cenf]=toimage(A,f);    % 将每个IMF信号合成求取Hilbert谱，E：对应的振幅值，Cenf：每个网格对应的中心频率。这里横轴为时间，纵轴为频率。即时频图(用颜色表示第三维值的大小)和三维图（三维坐标系：时间，中心频率，振幅）        
cemd_visu(data,1:length(data),imf); % 显示每个IMF分量及残余信号
disp_hhs(E);    % 希尔伯特谱
% 画出边际谱
colormap(flipud(gray)); % 黑白显示
%N=length(Cenf);    % 设置频率点数。完全从理论公式出发。网格化后中心频率很重要，大家从连续数据变为离散的角度去思考，相信应该很容易理解
for k=1:size(E,1)
    bjp(k)=sum(E(k,:))*1/fs;
end
figure(3);
plot(Cenf(1,:)*fs,bjp);  % 作边际谱图   进行求取Hilbert谱时频率已经被抽样成具有一定窗长的离散频率，所以此时的频率轴已经是中心频率
xlabel(\'频率 / Hz\');
ylabel(\'幅值\');

% 绘制瞬时包络图和瞬时频率图
figure;
subplot(221),plot(t/N,A(1,:));xlabel(\'时间 t/s\');ylabel(\'幅值\');title(\'imf1分量瞬时包络\');
subplot(222),plot(t/N,f(1,:)*fs);xlabel(\'时间 t/s\');ylabel(\'频率\');title(\'imf1分量瞬时频率\');
subplot(223),plot(t/N,A(2,:));xlabel(\'时间 t/s\');ylabel(\'幅值\');title(\'imf2分量瞬时包络\');
subplot(224),plot(t/N,f(2,:)*fs);xlabel(\'时间 t/s\');ylabel(\'频率\');title(\'imf2分量瞬时频率\');

Hilbert变换-Hilbert谱、Hilbert边际谱和包络谱：于青民

Hilbert谱：信号的希尔伯特变换后做fft，表示信号幅值在整个频率段上随时间和频率的变化规律；
Hilbert边际谱：对hilbert谱做积分，表示信号幅值在整个频率段上随频率的变化情况，它相当于傅里叶谱，但比傅里叶谱具有更高的频率分辨率。Hilbert边际谱是通过对Hilbert谱（在时间上）积分得到的；
Hilbert包络谱：希尔伯特变换后做包络后再fft，不同于Hilbert谱和Hilbert边际谱，是直接对信号进行Hilbert变换后构造解析函数，然后依据解析函数求模值，求的模值即为包络，然后对信号包络进行FFT后得到的即为Hilbert包络谱。

学习笔记：rogen

边际谱从统计意义上表示了整组数据每个频率点的积累幅值分布，而傅立叶谱的某点幅值表示在整个信号里有一个含有此频率的三角函数组分，而且幅值越大只是说明在整个数据段上，局部存在的可能性越大。
再看得到的图形，FFT 表示的是整个数据中，能量在一个频率上分布的可能性地描述，而边际谱表示在在每一个频率上幅值的积累，如果想知道具体时间那么就看HHT谱，这个时间-幅值-频率的三维谱。说到瞬时频率，傅立叶变换不强调局部性，而是强调全局性。咱们的HHT才提出一个唯一的瞬时频率的定义。因此拿瞬时频率来衡量傅立叶变换也是不公平的。
对于边际谱，就是hilbert谱对时间的积分，从积分的角度来讲，就相当于对任意一阶频率把所有的时间上的幅值都加起来了，这就反映这阶频率在所有时间上的幅值积累，而对于他们对于频率所对应的幅值，描述为该频率在整个时间上出现的可能，我个人认为是既然出现了，那就是存在的，不能说是一种可能，只能说他出现的次数的多少，经过累积以后就变成了他的幅值，而对于fourier来讲，他只能说明在用正余弦函数拟合这个信号的时候需要这一阶频率，幅值大，则说明他对拟合这个信号的贡献大，而不是一种出现的可能性的度量。
能量谱，可以理解为对边际谱的平方，这个只是具有能量的形式，而具体是否能代表能量，需要更进一步的探讨，还有一个瞬时能量，这些名词确实很诱人，但究竟如何，还需要大家的努力。

另一个比较详细的HHT资料，附matlab代码：cwjy