希尔伯特黄变换算例2

1、页眉内容电力工程信号处理应用希尔伯特黄变换【目的】1.了解希尔伯特黄变换的理论知识及应用领域2.用 Matlab 软件仿真，验证希尔伯特黄变换的优点【希尔伯特黄变换】希尔伯特黄变换(Hilbert-Hua ng tran sform, HHT)首先采用 EMD 方法将信号分解为若干个 IMF 分量之和，然后对每个 IMF 分量进行 Hilbert 变换得到的瞬时 频率和瞬时幅值，从而得到信号的 Hilbert 谱，Hilbert 谱表示了信号完整的时间- 频率分布，是具有一定的自适应的时频分析方法。与前面的小波分析方法相比， 避免了小波分析基选取的困难。分析非线性、非平稳信号采用基于经验模态分

2、解 的 HHT 方法可以较好地分析信号的局域动态行为和特征。由于 HHT 方法的种种特点，其在机械振动、生物医学、故障诊断、海洋学科、地震工程学以及经济学各学科中得到了广泛应用。在电力系统领域中，HHT 方法可用于谐波分析、同步电机参数辨识、低频震荡分析、电能质量检测、磁铁谐振 过电压辨识等方面和超高速方向保护等方面。HHT 方法在电力系统中的应用还在 进一步的研究和探索中。页眉内容【EMD分解】对于一个时间序列 x(t)，其经验模态分解过程如下：(1)确定原始信号 x(t )的所有极大值点和极小值点；(2)采用样条函数求出 x(t )的上、下包络线，并计算均值 m(t) ；(3)做差 h(t

3、) 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 所示。图 1 EMD 分解流程图对于分解总阶数为n的时间序列，最后可以表示成式中，r(t)为残余函数，它是以单调函数。【算例 2】考察两个函数x1A1sin(2 ?50?t) A2sin(2 ?100?t)(1)A1sin(2 ?5

4、0?t), t 0.1x21(2 )2A2sin(2 ?100?t), t 0.1编程：EMD 分解程序1.emd.mfunction 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).A2)/sum(x1.A2);页眉内容x1 = x2;end imfend+1 = x1; x = x-x1;end

5、 imfend+1 = x; function u = ismonotonic(x) u1 = length(findpeaks(x)*length(findpeaks(-x); if u1 0u = 0; elseu = 1;end function u = isimf(x) N = length(x);u1 = sum(x(1:N-1).*x(2:N) 1u = 0; elseu = 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.Findp

6、esks.mfunction n = findpeaks(x) n = find(diff(diff(x) 0) x(n); n(u) =n(u)+1;3.FFTAnalysis.mfunction Y, f = FFTAnalysis(y, Ts) Fs = 1/Ts;L = length(y);NFFT = 2 切 extpow2(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);endfunction yenvelope, yf, yh, yang

7、le = HilbertAnalysis(y, Ts) yh = hilbert(y);yenvelope = abs(yh); yangle = unwrap(angle(yh);页眉内容yf = diff(yangle)/2/pi/Ts;end4.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

8、_hht(x(1:6000),1/Fs);% Func : emd% imf = emd(x);for k = 1:length(imf)b(k) = sum(imfk.*imfk);th = unwrap(angle(hilbert(imfk); %dk = diff(th)/Ts/(2*pi); % endu,v = sort(-b);% Hilbert 瞬时频率图N = length(x);c = linspace(0,(N-2)*Ts,N-1);% 0:Ts:Ts*(N-2)for k = v(1:2)% 显示能量最大的两个 IMF 的瞬时频率figureplot(c,dk);xlim

9、(0 c(end);ylim(0 1/2/Ts);xlabel(Time/s)ylabel(Frequency/Hz);title(sprintf(IMF%d, k)end% 显示各 IMFM = length(imf);N = length(x);c = linspace(0,(N-1)*Ts,N);% 0:Ts:Ts*(N-1)for k1 = 0:4:M-1figurefor k2 = 1:min(4,M-k1) subplot(4,2,2*k2-1) plot(c,imfk1+k2)set(gca,FontSize,8,XLim,0 c(end);title(spri ntf(第 %d

10、f IMF, k1+k2)xlabel(Time/s) ylabel(sprintf(IMF%d, k1+k2);subplot(4,2,2*k2)yf, f = FFTAnalysis(imfk1+k2, Ts); plot(f, yf)相位瞬时频率b = 1-b/max(b); %后面绘图的亮度控制页眉内容title(sprintf(第%df IMF 的频谱,k1+k2)xlabel(f/Hz) ylabel(|IMF(f)|);endendfigure subplot(211) plot(c,x) set(gca,FontSize,8,XLim,0 c(end);title( 原始信号

11、)xlabel(Time/s) ylabel(Origin);subplot(212)Yf, f = FFTAnalysis(x, Ts);plot(f, Yf)title( 原始信号的频谱 ) xlabel(f/Hz) ylabel(|Y(f)|);指令窗口输入的代码 1.函数 x1A1sin(2 ?50?t) A2sin(2 ?100?t)代码1clc clear all close all % x, Fs = wavread(Hum.wav);% Ts = 1/Fs;% x = x(1:6000);Ts = 0.001;Fs = 1/Ts; t=0:Ts:0.2 x = sin(2*pi

12、*50*t) + sin(2*pi*100*t) ;%(式 1 )imf = emd(x);plot_hht(x,imf,1/Fs);k = 4;y = imfk;N = length(y);t = 0:2*Ts:Ts*(N-1);yenvelope, yfreq, yh, yangle = HilbertAnalysis(y, 1/Fs); yModulate = y./yenvelope;YMf, f = FFTAnalysis(yModulate, Ts);Yf = FFTAnalysis(y, Ts); figure subplot(321) plot(t, y) title(spri

13、ntf(IMF%d, k) xlabel(Time/s)ylabel(sprintf(IMF%d, k);subplot(322) plot(f, Yf) title(sprintf(IMF%d 的频谱 , k) xlabel(f/Hz) ylabel(|IMF(f)|);subplot(323) plot(t, yenvelope) title(sprintf(IMF%d 的包络 , k) xlabel(Time/s)ylabel(envelope);subplot(324) plot(t(1:end-1), yfreq) title(sprintf(IMF%d 的瞬时频率 , k) xla

14、bel(Time/s)ylabel(Frequency/Hz);subplot(325) plot(t, yModulate) title(sprintf(IMF%d 的调制信号 , k) xlabel(Time/s)ylabel(modulation);subplot(326) plot(f, YMf) title(sprintf(IMF%d 调制信号的频谱 , k) xlabel(f/Hz)页眉内容ylabel（|YMf|）;仿真得到的图形如下图 1 和图 2代码和代码 1 几乎相同，只需将式一改为下式x = sin( 2*pi*50*t) .*(t=0.1)得到的波形如下图 3 和图 4图 3【学习心得】2.函数 X2般总,t 0.1t 0.1FM50-5；-x 八Rm:;:：jL / / /II 1 f 4vIM1iLIii yW M 1i y VU il艸200400f/Hz第 2 个 IMF 的频谱6002LIM0200400f