前言

在希爾伯特變換中我們已經知道淤齐，直接對"解析信號"做3瞬屬性的分析時，"瞬時頻率"會出現很多負頻率袜匿，是無意義的更啄！因此為解決這個問題并對非穩(wěn)態(tài)信號做時頻分析，才有了"希爾伯特-黃變換"沉帮！

思路：在對原始信號進行希爾伯特變換之前锈死，要先把信號分解成瞬時頻率具有意義的各個分量/函數！把信號進行分解的方法穆壕，在希爾伯特-黃變換理論中稱為"經驗模態(tài)分解"待牵；分解出來的各個分量，其實是一類函數喇勋，稱為"固有模態(tài)函數"缨该，利用這類函數的局部特性，可以使得函數在任意一點的瞬時頻率都有意義川背。

求瞬時頻率(無意義)
$\downarrow$
經驗模態(tài)分解贰拿，得到一類固有模態(tài)函數
$\downarrow$
對固有模態(tài)函數做希爾伯特變換蛤袒，求瞬時頻率(有意義)

所以，對希爾伯特-黃變換做簡要總結(其實一點都不復雜)：
目標：求原信號的"有物理意義"的瞬時頻率膨更；
步驟：1. 經驗模態(tài)分解妙真；2. 希爾伯特變換/譜分析。

下面對2大步進行詳細的說明荚守，并給出對應的matlab手寫程序珍德。

經驗模態(tài)分解

為什么希爾伯特-黃變換那么好？
個人認為它的關鍵就在于"經驗模態(tài)分解"這一步矗漾。因為這種信號的分解锈候，完全是"方法適應于數據"的，即不同的數據敞贡，方法會根據數據的情況做變換和調整泵琳，而不是始終用固定的那幾種基函數(小波變換、短時傅里葉變換)誊役。
這種自適應思想下的工具获列，一定是非常優(yōu)越的！
什么是固有模態(tài)函數势木？
經驗模態(tài)分解的目標就是將原始信號分解成一類固有模態(tài)函數蛛倦，
固有模態(tài)函數其實只是滿足下面兩個條件的函數而已歌懒，沒啥復雜的：

函數極值點的數目與零點的數目必須相等啦桌，或最多相差1；
極大值和極小值形成的包絡及皂，包絡的均值為0甫男；

其中均值為0的條件對應離散數據不好實現，一般平均值趨近于0（一般和0做差<0.1）即可验烧。

經驗模態(tài)分解何時結束板驳？
經驗模態(tài)分解就是為了分解出一系列的固有模態(tài)函數，所以當剩余數據已經無法再分出來固有模態(tài)函數時碍拆，經驗模態(tài)分解這一步就徹底結束了若治。
如何判斷剩余數據無法再分解了？
剩余數據極大值或極小值點數目有一個為0時(均值條件沒法滿足了)感混，或剩余數據已經是單調時(沒法做包絡線了)端幼，就可以結束了。

分解的流程

設原始信號為 $x(t)$ 弧满，分解步驟如下：

步1：提取輸入信號的所有極大值和極小值點婆跑，并用三次樣條曲線連接出上、下包絡線庭呜，設包絡線均值為 $m(t)$ 滑进，用輸入信號減去均值可得 $h(t)$ ：

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

步2：判斷步1得到的 $h(t)$ 是否滿足固有模態(tài)函數的2個條件犀忱，若不滿足把 $h(t)$ 當作輸入回到步1；若滿足則得到一個固有模態(tài)函數扶关，進入步3:阴汇；
步3：設得到的第k個固有模態(tài)函數為 $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：判斷新的剩余項是否滿足經驗模態(tài)分解的結束條件鲫寄，不滿足將剩余項帶回到步1；滿足則結束經驗模態(tài)分解疯淫。

分解完畢后地来，原始信號可以用n個固有模態(tài)函數 + 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

主調用函數：固有模態(tài)函數用二維矩陣記錄

% 主函數: 經驗模態(tài)分解

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內層循環(huán)
    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;  % 最后一個余項也放進去


% 經驗模態(tài)分解出的函數(最后一個是余項)
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('經驗模態(tài)分解函數');

經驗模態(tài)分解效果：