希尔伯特黄变换算例2(8页).doc

《希尔伯特黄变换算例2(8页).doc》由会员分享，可在线阅读，更多相关《希尔伯特黄变换算例2(8页).doc（8页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

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

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

3、 令，即为一个IMF分量，做差；（6） 是否满足终止条件，若不满足则将作为新的输入信号转至第（1）步，若满足则EMD分解过程结束，不能提取的为残余量。具体流程如图1所示。图1 EMD分解流程图对于分解总阶数为的时间序列，最后可以表示成式中，为残余函数，它是以单调函数。【算例2】考察两个函数 （1） （2）编程：EMD分解程序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 =

4、 -getspline(-x1); x2 = x1-(s1+s2)/2; sd = sum(x1-x2).2)/sum(x1.2); x1 = x2; end imfend+1 = x1; x = x-x1;endimfend+1 = x;function u = ismonotonic(x)u1 = length(findpeaks(x)*length(findpeaks(-x);if u1 0 u = 0;else u = 1;endfunction u = isimf(x)N = length(x);u1 = sum(x(1:N-1).*x(2:N) 1 u = 0;else u = 1

5、;endfunction s = getspline(x)N = length(x);p = findpeaks(x);s = spline(0 p N+1,0 x(p) 0,1:N);1. Ffunction n = findpeaks(x)n = find(diff(diff(x) 0) x(n);n(u) = n(u)+1; function Y, f = FFTAnalysis(y, Ts)Fs = 1/Ts;L = length(y);NFFT = 2nextpow2(L);y = y - mean(y);Y = fft(y, NFFT)/L;Y = 2*abs(Y(1:NFFT/2

6、+1);f = Fs/2*linspace(0, 1, NFFT/2+1);endfunction yenvelope, yf, yh, yangle = HilbertAnalysis(y, Ts)yh = hilbert(y);yenvelope = abs(yh); yangle = unwrap(angle(yh); yf = diff(yangle)/2/pi/Ts; endfunction plot_hht(x,imf,Ts)% Plot the HHT.% : Syntax% The array x is the input signal and Ts is the sampli

7、ng 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(imfk.*imfk); th = unwrap(angle(hilbert(imfk); % 相位 dk = diff(th)/Ts/(2*pi); % 瞬时频率endu,v = sort(-b);b = 1-b/max(b); % 后面绘图的亮度控制% Hilbert瞬时频率图N = length(x);c = li

8、nspace(0,(N-2)*Ts,N-1); % 0:Ts:Ts*(N-2)for k = v(1:2) % 显示能量最大的两个IMF的瞬时频率 figure plot(c,dk); xlim(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-1 figure for k2 =

9、 1:min(4,M-k1) subplot(4,2,2*k2-1) plot(c,imfk1+k2) set(gca,FontSize,8,XLim,0 c(end); title(sprintf(第%d个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) title(sprintf(第%d个IMF的频谱, k1+k2) xlabel(f/Hz) ylabel(|IMF(f)|); endendfigu

10、resubplot(211)plot(c,x)set(gca,FontSize,8,XLim,0 c(end);title(原始信号)xlabel(Time/s)ylabel(Origin);subplot(212)Yf, f = FFTAnalysis(x, Ts);plot(f, Yf)title(原始信号的频谱)xlabel(f/Hz)ylabel(|Y(f)|);指令窗口输入的代码1.函数代码1clcclear allclose all% x, Fs = wavread(Hum.wav);% Ts = 1/Fs;% x = x(1:6000);Ts = 0.001;Fs = 1/Ts;

11、t=0:Ts:x = sin(2*pi*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);figuresubplot(321)plot(t,

12、 y)title(sprintf(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)xlabel(

13、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(f)|);仿真得到的图形如下图1和图2图1图22.函数代码和代码1几乎相同，只需将式一改为下式x = sin(2*pi*50*t) .*(t=0.1)得到的波形如下图3和图4图3图4【学习心得】-第 6 页