希尔伯特黄变换算例2.docVIP

电力工程信号处理应用 希尔伯特黄变换 【目的】 了解希尔伯特黄变换的理论知识及应用领域 用Matlab软件仿真，验证希尔伯特黄变换的优点 【希尔伯特黄变换】 希尔伯特黄变换（Hilbert-Huang transform, HHT）首先采用EMD方法将信号分解为若干个IMF分量之和，然后对每个IMF分量进行Hilbert变换得到的瞬时频率和瞬时幅值，从而得到信号的Hilbert谱，Hilbert谱表示了信号完整的时间-频率分布，是具有一定的自适应的时频分析方法。与前面的小波分析方法相比，避免了小波分析基选取的困难。分析非线性、非平稳信号采用基于经验模态分解的HHT方法可以较好地分析信号的局域动态行为和特征。 由于HHT方法的种种特点，其在机械振动、生物医学、故障诊断、海洋学科、地震工程学以及经济学各学科中得到了广泛应用。在电力系统领域中，HHT方法可用于谐波分析、同步电机参数辨识、低频震荡分析、电能质量检测、磁铁谐振过电压辨识等方面和超高速方向保护等方面。HHT方法在电力系统中的应用还在进一步的研究和探索中。 【EMD分解】 对于一个时间序列,其经验模态分解过程如下： 确定原始信号的所有极大值点和极小值点； 采用样条函数求出的上、下包络线，并计算均值； 做差； 是否满足终止条件，若不满足将作为新的输入信号转至第（1）步，否则转为第（5）步； 令，即为一个IMF分量，做差； 是否满足终止条件，若不满足则将作为新的输入信号转至第（1）步，若满足则EMD分解过程结束，不能提取的为残余量。具体流程如图1所示。 图1 EMD分解流程图 对于分解总阶数为的时间序列，最后可以表示成 式中，为残余函数，它是以单调函数。 【算例2】 考察两个函数 （1） （2） 编程： [EMD分解程序] 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); 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; 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);