希尔伯特黄变换算例.doc
【目的】电力工程信号处理应用希尔伯特黄变换1 了解希尔伯特黄变换的理论知识及应用领域2 用Matlab软件仿真，验证希尔伯特黄变换的优点【希尔伯特黄变换】希尔伯特黄变换Hilbert-Huang transform, HHT首先采用EMD方法将信号分解为假设干个IMF分量之与，然后对每个IMF分量进展Hilbert变换得到的瞬时频率与瞬时幅值，从而得到信号的Hilbert谱，Hilbert谱表示了信号完整的时间-频率分布，是具有一定的自适应的时频分析方法。与前面的小波分析方法相比，防止了小波分析基选取的困难。分析非线性、非平稳信号采用基于经历模态分解的HHT方法可以较好地分析信号的局域动态行为与特征。由于HHT方法的种种特点，其在机械振动、生物医学、故障诊断、海洋学科、地震工程学以及经济学各学科中得到了广泛应用。在电力系统领域中，HHT方法可用于谐波分析、同步电机参数辨识、低频震荡分析、电能质量检测、磁铁谐振过电压辨识等方面与超高速方向保护等方面。HHT方法在电力系统中的应用还在进一步的研究与探索中。【EMD分解】对于一个时间序列,其经历模态分解过程如下：（1） 确定原始信号的所有极大值点与极小值点；（2） 采用样条函数求出的上、下包络线，并计算均值；（3） 做差；（4） 是否满足终止条件，假设不满足将作为新的输入信号转至第1步，否那么转为第5步；（5） 令，即为一个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 = -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) < 0); u2 = length(findpeaks(x)+length(findpeaks(-x); if abs(u1-u2) > 1 u = 0;else u = 1;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) < 0);u = find(x(n+1) > 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+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 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(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 = linspace(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 = 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)|'); endendfiguresubplot(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;t=0:Ts:x = sin(2*pi*50*t) + sin(2*pi*100*t) ; %式1imf = 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, 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('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)+sin(2*pi*100*t).*(t>=0.1)得到的波形如下列图3与图4图3图4【学习心得】第 11 页
