希尔伯特黄变换算例2

电力工程信号处理应用

希尔伯特黄变换

【目的】

1．了解希尔伯特黄变换的理论知识及应用领域

2．用Matlab软件仿真，验证希尔伯特黄变换的优点

【希尔伯特黄变换】

希尔伯特黄变换（Hilbert-Huang transform, HHT）首先采用EMD方法将信号分解为若干个IMF分量之和，然后对每个IMF分量进行Hilbert变换得到的瞬时频率和瞬时幅值，从而得到信号的Hilbert谱，Hilbert谱表示了信号完整的时间-频率分布，是具有一定的自适应的时频分析方法。与前面的小波分析方法相比，避免了小波分析基选取的困难。分析非线性、非平稳信号采用基于经验模态分解的HHT方法可以较好地分析信号的局域动态行为和特征。

由于HHT方法的种种特点，其在机械振动、生物医学、故障诊断、海洋学科、地震工程学以及经济学各学科中得到了广泛应用。在电力系统领域中，HHT方法可用于谐波分析、同步电机参数辨识、低频震荡分析、电能质量检测、磁铁谐振过电压辨识等方面和超高速方向保护等方面。HHT方法在电力系统中的应用还在进一步的研究和探索中。

【EMD 分解】

对于一个时间序列()x t ,其经验模态分解过程如下： （1） 确定原始信号()x t 的所有极大值点和极小值点；

（2） 采用样条函数求出()x t 的上、下包络线，并计算均值()m t ； （3） 做差()()()h t x t m t =-；

（4） ()h t 是否满足终止条件，若不满足将()h t 作为新的输入信号转至第（1）

步，否则转为第（5）步；

（5） 令()c h t =，c 即为一个IMF 分量，做差()r x t c =-；

（6） r 是否满足终止条件，若不满足则将r 作为新的输入信号转至第（1）步，

若满足则EMD 分解过程结束，不能提取的为残余量。具体流程如图1所示。

EMD 分解过程

图1 EMD 分解流程图

对于分解总阶数为n 的时间序列，最后可以表示成

1()()()n

i i x t c t r t ==+∑

式中，()r t 为残余函数，它是以单调函数。

【算例2】

考察两个函数

112sin(250)sin(2100)x A t A t ππ=∙∙+∙∙ （1） 1

22sin(250),0.1sin(2100),

0.1

A t t x A t t ππ∙∙<⎧=⎨∙∙≥⎩ （2）

编程：

[EMD 分解程序]

1. emd.m

function imf = emd(x) x = transpose(x(:)); imf = [];

while ~ismonotonic(x) x1 = x; sd = Inf;

while (sd > 0.1) || ~isimf(x1) s1 = getspline(x1); s2 = -getspline(-x1); x2 = x1-(s1+s2)/2;

sd = sum((x1-x2).^2)/sum(x1.^2); x1 = x2; end

imf{end+1} = x1; x = x-x1; end

imf{end+1} = x;

function u = ismonotonic(x)

u1 = length(findpeaks(x))*length(findpeaks(-x));

if u1 > 0

u = 0;

else

u = 1;

end

function u = isimf(x)

N = length(x);

u1 = sum(x(1:N-1).*x(2:N) < 0);

u2 = length(findpeaks(x))+length(findpeaks(-x));

if abs(u1-u2) > 1

u = 0;

else

u = 1;

end

function s = getspline(x)

N = length(x);

p = findpeaks(x);

s = spline([0 p N+1],[0 x(p) 0],1:N);

2.Findpesks.m

function n = findpeaks(x)

n = find(diff(diff(x) > 0) < 0);

u = find(x(n+1) > x(n));

n(u) = n(u)+1;

3.FFTAnalysis.m

function [Y, f] = FFTAnalysis(y, Ts)

Fs = 1/Ts;

L = length(y);

NFFT = 2^nextpow2(L);

y = y - mean(y);

Y = fft(y, NFFT)/L;

Y = 2*abs(Y(1:NFFT/2+1));

f = Fs/2*linspace(0, 1, NFFT/2+1);

end

function [yenvelope, yf, yh, yangle] = HilbertAnalysis(y, Ts) yh = hilbert(y);

yenvelope = abs(yh);

yangle = unwrap(angle(yh));

yf = diff(yangle)/2/pi/Ts;

end

4.pot_hht.m

function plot_hht(x,imf,Ts)

% Plot the HHT.

% :: Syntax

% The array x is the input signal and Ts is the sampling period.

% Example on use: [x,Fs] = wavread('Hum.wav');

% plot_hht(x(1:6000),1/Fs);

% Func : emd

% imf = emd(x);

for k = 1:length(imf)

b(k) = sum(imf{k}.*imf{k});

th = unwrap(angle(hilbert(imf{k}))); % 相位

d{k} = diff(th)/Ts/(2*pi); % 瞬时频率

end

[u,v] = sort(-b);

b = 1-b/max(b); % 后面绘图的亮度控制

% Hilbert瞬时频率图

N = length(x);

c = linspace(0,(N-2)*Ts,N-1); % 0:Ts:Ts*(N-2)

for k = v(1:2) % 显示能量最大的两个IMF的瞬时频率

figure

plot(c,d{k});

xlim([0 c(end)]);

ylim([0 1/2/Ts]);

xlabel('Time/s')

ylabel('Frequency/Hz');

title(sprintf('IMF%d', k))

end

% 显示各IMF

M = length(imf);

N = length(x);

c = linspace(0,(N-1)*Ts,N); % 0:Ts:Ts*(N-1)

for k1 = 0:4:M-1

figure