离散数据希尔伯特-黄变换

2019-05-16 本文已影响0人胜负55开

前言

在希尔伯特变换中我们已经知道，直接对"解析信号"做3瞬属性的分析时，"瞬时频率"会出现很多负频率，是无意义的！因此为解决这个问题并对非稳态信号做时频分析，才有了"希尔伯特-黄变换"！

思路：在对原始信号进行希尔伯特变换之前，要先把信号分解成瞬时频率具有意义的各个分量/函数！把信号进行分解的方法，在希尔伯特-黄变换理论中称为"经验模态分解"；分解出来的各个分量，其实是一类函数，称为"固有模态函数"，利用这类函数的局部特性，可以使得函数在任意一点的瞬时频率都有意义。

求瞬时频率(无意义)
$\downarrow$
经验模态分解，得到一类固有模态函数
$\downarrow$
对固有模态函数做希尔伯特变换，求瞬时频率(有意义)

所以，对希尔伯特-黄变换做简要总结(其实一点都不复杂)：
目标：求原信号的"有物理意义"的瞬时频率；
步骤：1. 经验模态分解；2. 希尔伯特变换/谱分析。

下面对2大步进行详细的说明，并给出对应的matlab手写程序。

经验模态分解

为什么希尔伯特-黄变换那么好？
个人认为它的关键就在于"经验模态分解"这一步。因为这种信号的分解，完全是"方法适应于数据"的，即不同的数据，方法会根据数据的情况做变换和调整，而不是始终用固定的那几种基函数(小波变换、短时傅里叶变换)。
这种自适应思想下的工具，一定是非常优越的！
什么是固有模态函数？
经验模态分解的目标就是将原始信号分解成一类固有模态函数，
固有模态函数其实只是满足下面两个条件的函数而已，没啥复杂的：

函数极值点的数目与零点的数目必须相等，或最多相差1；
极大值和极小值形成的包络，包络的均值为0；

其中均值为0的条件对应离散数据不好实现，一般平均值趋近于0（一般和0做差<0.1）即可。

经验模态分解何时结束？
经验模态分解就是为了分解出一系列的固有模态函数，所以当剩余数据已经无法再分出来固有模态函数时，经验模态分解这一步就彻底结束了。
如何判断剩余数据无法再分解了？
剩余数据极大值或极小值点数目有一个为0时(均值条件没法满足了)，或剩余数据已经是单调时(没法做包络线了)，就可以结束了。

分解的流程

设原始信号为 $x(t)$ ，分解步骤如下：

步1：提取输入信号的所有极大值和极小值点，并用三次样条曲线连接出上、下包络线，设包络线均值为 $m(t)$ ，用输入信号减去均值可得 $h(t)$ ：

$h(t) = x(t) - m(t)$

步2：判断步1得到的 $h(t)$ 是否满足固有模态函数的2个条件，若不满足把 $h(t)$ 当作输入回到步1；若满足则得到一个固有模态函数，进入步3:；
步3：设得到的第k个固有模态函数为 $h_{k}(t)$ ，将其赋值给 $c_{k}(t)$ ：

$c_{k}(t) = h_{k}(t)$

将 $c_{k}(t)$ 从剩余原始序列中分离出来，得到新的剩余项：

$r_{k}(t) = x_{k}(t) - c_{k}(t)$

步4：判断新的剩余项是否满足经验模态分解的结束条件，不满足将剩余项带回到步1；满足则结束经验模态分解。

分解完毕后，原始信号可以用n个固有模态函数 + 1个剩余项的形式表示：

$x(t) = \sum_{i=1}^{n}c_{i} + r_{n}$

下面给出各步的matlab程序：

功能函数：

% 非主函数: 求极值点个数
% 注意: 若要极大值点数直接输入x，若要极小值则输入-x
function n = peaks(x)

[value,local] = findpeaks(x);

n = length(local);  % 只要极值点个数

% 非主函数: 获得数据的极大值、极小值包络线(三次样条插值后的数据)
function s = getspline(x,t)

N = length(x);

[value,local] = findpeaks(x);

% 原始插值点: 极值点
t1 = t(local);
x1 = x(local);

s = spline(t1,x1,t);

% 非主函数: 判断当前函数是否为imf —— 零点个数与极值点个数对比
function u = isimf(x)

N = length(x);
u1 = sum( x(1:N-1).*x(2:N) < 0 );  % 零点个数
u2 = peaks(x) + peaks(-x);  % 极大值 + 极小值总点数

if abs(u2-u1) > 1
    u = 0;  % 不是imf函数
else
    u = 1;  % 是imf函数
end

% 非主函数: 剩余信号的单调性判断(是否整个程序结束)
function u = ismono(x)

u1 = peaks(x)*peaks(-x);  

if u1 ~= 0
    u = 0;  % 非单调: 不存在极大值点或极小值点
else
    u = 1;  % 单调：有一个极值不是0，就不是单调的
end

主调用函数：固有模态函数用二维矩阵记录

% 主函数: 经验模态分解

clear; clc;

x = xlsread('shuju.xlsx');
x = x(1001:1001+1023)';  % 有效数据长必须是2^n，所以我取1024，最后10几个点是0
N = length(x);
fs = 100;          % 采样频率 = 1/采样间隔
t = (0:N-1)/fs;   % 时间刻度

imf = [];

num = 1;  % 计数器
xx = x;   % 代替原始数据被操作而已
while ~ismono(xx)
    x1 = xx;
    sd = inf;  % 用sd代替包络均值为0(<0.01)的条件
    % 当sd不满足<0.01，或x1不是imf函数时，进入while内层循环
    while (sd > 0.01) || ~isimf(x1)
        s1 = getspline(x1,t);    % 极大值包络插值结果
        s2 = -getspline(-x1,t);  % 极小值包络插值结果
        x2 = x1 - (s1+s2)/2;
        
        sd = sum( (x1-x2).^2 )/sum(x1.^2);
        x1 = x2;
    end
    
    imf(num,:) = x1;
    xx = xx - x1;
    num = num + 1;
end

imf(num,:) = xx;  % 最后一个余项也放进去


% 经验模态分解出的函数(最后一个是余项)
figure(1);
merge = 0;
[row,col] = size(imf);
for n = 1:row
    % 把经验分解的内容再加起来 
    merge = merge + imf(n,:);
    % 每一个画图
    subplot(row,1,n)
    plot(t,imf(n,:));
    axis([min(t) max(t) -inf inf]);
end
suptitle('经验模态分解函数');

经验模态分解效果：