
function imf = emd(x)
% Empiricial Mode Decomposition (Hilbert-Huang Transform)
% imf = emd(x)
% Func : findpeaks
x= transpose(x(:))%转置基颂为行矩阵
imf = []
while ~ismonotonic(x) %当x不是单调函数,分解终止条件
x1 = x
sd = Inf%均值
%直到x1满足IMF条件,得c1
while (sd >0.1) | ~isimf(x1) %当标准偏差系数sd大于0.1或x1不是固有模态函数时,分量终搏肆郑止条件
s1 = getspline(x1)%上包络线
s2 = -getspline(-x1)%下包络线
x2 = x1-(s1+s2)/2%此处的x2为文章中的h
sd = sum((x1-x2).^2)/sum(x1.^2)
x1 = x2
end
imf{end+1} = x1
x = x-x1
end
imf{end+1} = x
% FUNCTIONS
function u = ismonotonic(x)
%u=0表示x不是单调函数,u=1表示x为单调的
u1 = length(findpeaks(x))*length(findpeaks(-x))
if u1 >0, u = 0
else, u = 1end
function u = isimf(x)
%u=0表示x不是固有模式函数,u=1表示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 = 1end
function s = getspline(x)
%三次样条函数拟合成元数据包络线
N = length(x)
p = findpeaks(x)
s = spline([0 p N+1],[0 x(p) 0],1:N)
-------------------------------------------------------------------------------
--------------------------------------------------------------------------------
function n = findpeaks(x)
% Find peaks.找到极值 ,n为极值点所在位置
% n = findpeaks(x)
n= find(diff(diff(x) >0) <0)
u= find(x(n+1) >x(n))
n(u) = n(u)+1
------------------------------------------------------------------------------------------
----------------------------------------------------------------------------------------
function plot_hht00(x,Ts)
% 双雹首边带调幅信号的EMD分解
% Plot the HHT.
% plot_hht(x,Ts)
%
% :: 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
% Get HHT.
clear all
close all
Ts=0.0005
t=0:Ts:10 % 采样率2000HZ
% 调幅信号
%x=sin(2*pi*t).*sin(40*pi*t)
x=sin(2*pi*t)
s1 = getspline(x)%上包络线
s2 = -getspline(-x)%上包络线
x1 = (s1+s2)/2%此处的x2为文章中的h
figure
plot(t,x)xlabel('Time'), ylabel('Amplitude')title('双边带调幅信号')hold on
plot(t,s1,'-r')
plot(t,s2,'-r')
plot(t,x1,'g')
imf = emd(x)
for k = 1:length(imf)
b(k) = sum(imf{k}.*imf{k})
th = angle(hilbert(imf{k}))
d{k} = diff(th)/Ts/(2*pi)
end
[u,v] = sort(-b)
b = 1-b/max(b)
% Set time-frequency plots.
N = length(x)
c = linspace(0,(N-2)*Ts,N-1)
%
figure
for k = v(1:2)
plot(c,d{k},'k.','Color',b([k k k]),'MarkerSize',3)hold on
set(gca,'FontSize',8,'XLim',[0 c(end)],'YLim',[0 50])%设置x、y轴句柄
xlabel('Time'), ylabel('Frequency')title('原信号时频图')
end
% Set IMF plots.
M = length(imf)
N = length(x)
c = linspace(0,(N-1)*Ts,N)
for k1 = 0:4:M-1
figure
for k2 = 1:min(4,M-k1),
subplot(4,1,k2),
plot(c,imf{k1+k2})
set(gca,'FontSize',8,'XLim',[0 c(end)])
title('EMD分解结果')
end
xlabel('Time')
end
局部均值局禅分解(Local Mean Decomposition) 是由Smith提出的一种新的非线性和非平桐吵尘稳信号分析方法。由于LMD是依据信号本身的信息进行自适应分解的,产生的PF分量具有真实的物理意义,由此得到的时频分布能够清晰准确地反映出信号能量在空间各尺度上的分布规律。编辑本段电子信息技术名词
人生模式数据库(Life Model Database)
编辑本段医学名词
二甲基二硫代氨基甲酸铅(lead dimethyl dithiocarbamate)
编辑本段法国教育制度模式
LMD学制为法国2004年开始实行的新教学制度。已被欧洲普遍接受,碰冲与国际教学体制接轨。 L:Licence (学士) 三年制 M:Masteur(硕士) 二年制 D:Docteur (博士) 三年制 现法国公立学校都已采用LMD教学模式。部分学校在Licence的第二年仍然会颁发原DEUG文凭。 仍有部分私立或培训学校颁发DUT或BTS文凭。 工程师学校一般持续5年,相当于Licence+Masteur。
编辑本段网络释义
欢迎分享,转载请注明来源:内存溢出
微信扫一扫
支付宝扫一扫
评论列表(0条)