《希尔伯特黄变换算例.doc》由会员分享，可在线阅读，更多相关《希尔伯特黄变换算例.doc（12页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、【目的】电力工程信号处理应用希尔伯特黄变换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） 令，即为一个IMF分量，做
3、差；（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)
4、; 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;endfunction s =
5、 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+1);f = Fs/2*lin
6、space(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.% Exam
7、ple 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)*T
8、s,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) s
9、ubplot(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)pl
10、ot(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
11、*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(I
12、MF%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)得到的波形如下列图3与图4图3图4【学习心得】第 11 页