第1页 / 共9页
第2页 / 共9页
第3页 / 共9页
第4页 / 共9页
第5页 / 共9页
第6页 / 共9页
第7页 / 共9页
第8页 / 共9页
EMD分解的MATLAB代码.docx
test.m 文件
clc
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:1;
x = sin(2*pi*10*t) + sin(2*pi*50*t) + sin(2*pi*100*t) + 0.1*randn(1, length(t));
imf = emd(x);
plot_hht(x,imf,1/Fs);
k = 4;
y = imf{k};
N = length(y);
t = 0: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(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)|');
findpeaks.m 文件
function n = findpeaks(x)
% Find peaks. 找极大值点，返回对应极大值点的坐标
n
u
= find(diff(diff(x) > 0) < 0); % 相当于找二阶导小于 0 的点
= find(x(n+1) > x(n));
% 加 1 才真正对应极大值
n(u) = n(u)+1;
点
% 图形解释上述过程
% figure
% subplot(611)
% x = x(1:100);
% plot(x, '-o')
% grid on
%
% subplot(612)
% plot(1.5:length(x), diff(x) > 0, '-o')
% grid on
% axis([1,length(x),-0.5,1.5])
%
% subplot(613)
% plot(2:length(x)-1, diff(diff(x) > 0), '-o')
% grid on
% axis([1,length(x),-1.5,1.5])
%
% subplot(614)
% plot(2:length(x)-1, diff(diff(x) > 0)<0, '-o')
% grid on
% axis([1,length(x),-1.5,1.5])
%
% n
= find(diff(diff(x) > 0) < 0);
% subplot(615)
% plot(n, ones(size(n)), 'o')
% grid on
% axis([1,length(x),0,2])
%
% u
= find(x(n+1) > x(n));
% n(u) = n(u)+1;
% subplot(616)
% plot(n, ones(size(n)), 'o')
% grid on
% axis([1,length(x),0,2])
plot_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_hht(x(1:6000),1/Fs);
% Func : emd
% imf = emd(x);
for k = 1:length(imf)
b(k) = sum(imf{k}.*imf{k});
th
= unwrap(angle(hilbert(imf{k})));
% 相位
d{k} = diff(th)/Ts/(2*pi);
% 瞬时频率
end
[u,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,d{k});
xlim([0 c(end)]);
ylim([0 1/2/Ts]);
xlabel('Time/s')
ylabel('Frequency/Hz');
title(sprintf('IMF%d', k))
end
% 显示各 IMF
M = 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,imf{k1+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(imf{k1+k2}, Ts);
plot(f, yf)
title(sprintf('第%d 个 IMF 的频谱', k1+k2))
xlabel('f/Hz')
ylabel('|IMF(f)|');
end
end
figure
subplot(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)|');
emd.m 文件
function imf = emd(x)
% Empiricial Mode Decomposition (Hilbert-Huang Transform)
% EMD 分解或 HHT 变换
% 返回值为 cell 类型，依次为一次 IMF、二次 IMF、...、最后残差
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
imf{end+1} = x1;
x
end
imf{end+1} = x;
% 是否单调
= x-x1;
function u = ismonotonic(x)
u1 = length(findpeaks(x))*length(findpeaks(-x));
if u1 > 0
u = 0;
u = 1;
else
end
% 是否 IMF 分量
function 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;
u = 1;
else
end
% 据极大值点构造样条曲线
function s = getspline(x)
N = length(x);
p = findpeaks(x);
s = spline([0 p N+1],[0 x(p) 0],1:N);
FFTAnalysis.m 文件
% 频谱分析
function [Y, f] = FFTAnalysis(y, Ts)
Fs = 1/Ts;
L = length(y);
NFFT = 2^nextpow2(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);
end
HilbertAnalysis.m 文件
% Hilbert 分析
function [yenvelope, yf, yh, yangle] = HilbertAnalysis(y, Ts)
yh = hilbert(y);
yenvelope = abs(yh);
yangle = unwrap(angle(yh));
yf = diff(yangle)/2/pi/Ts;
end
ＥＥＭＤ分解代码
% 包络
% 相位
% 瞬时频率
Y: Inputted data;1-d data only
Nstd: ratio of the standard deviation of the added noise and that of Y;
NE: Ensemble number for the EEMD
A matrix of N*(m+1) matrix, where N is the length of the input
data Y, and m=fix(log2(N))-1. Column 1 is the original data, columns 2, 3, ...
m are the IMFs from high to low frequency, and comlumn (m+1) is the
residual (over all trend).
%function allmode=eemd(Y,Nstd,NE)
%
% This is an EMD/EEMD program
%
% INPUT:
%
%
%
% OUTPUT:
%
%
%
%
%
% NOTE:
%
%
%
%
% References:
% Wu, Z., and N. E Huang (2008),
% Ensemble Empirical Mode Decomposition: a noise-assisted data analysis method.
%
%
% code writer: Zhaohua Wu.
% footnote:S.C.Su 2009/03/04
It should be noted that when Nstd is set to zero and NE is set to 1, the
program degenerates to a EMD program.(for EMD Nstd=0,NE=1)
This code limited sift number=10 ,the stoppage criteria can't change.
Advances in Adaptive Data Analysis. Vol.1, No.1. 1-41.
TNM2=TNM+2,original data and residual included in TNM2
assign 0 to TNM2 matrix
4.add noise
5.give initial values before sift
6.start to find an IMF------------------------------------------------IMF loop start
7.sift 10 times to get IMF--------------------------sift loop start
and end
8.after 10 times sift --we got IMF
9.subtract IMF from data ,and let the residual to find next IMF by loop
6.after having all the IMFs---------------------------------------------IMF loop end
9.after TNM IMFs ,the residual xend is over all trend
%
% There are three loops in this code coupled together.
% 1.read data, find out standard deviation ,devide all data by std
% 2.evaluate TNM as total IMF number--eq1.
%
%
% 3.Do EEMD NE times-------------------------------------------------------------loop EEMD start
%
%
%
%
%
%
%
%
% 3.Sum up NE decomposition result-------------------------------------------------loop EEMD end
% 10.Devide EEMD summation by NE,std be multiply back to data
%
% Association: no
% this function ususally used for doing 1-D EEMD with fixed
% stoppage criteria independently.
%
% Concerned function: extrema.m
%
above mentioned m file must be put together
function allmode=eemd(Y,Nstd,NE)
%part1.read data, find out standard deviation ,devide all data by std
xsize=length(Y);
dd=1:1:xsize;
Ystd=std(Y);
Y=Y/Ystd;
%part2.evaluate TNM as total IMF number,ssign 0 to TNM2 matrix
TNM=fix(log2(xsize))-1;
TNM2=TNM+2;
for kk=1:1:TNM2,
for ii=1:1:xsize,
allmode(ii,kk)=0.0;
end
end
%part3 Do EEMD -----EEMD loop start
for iii=1:1:NE,
%EEMD loop -NE times EMD sum together
%part4 --Add noise to original data,we have X1
for i=1:xsize,
temp=randn(1,1)*Nstd;
X1(i)=Y(i)+temp;
end
%part4 --assign original data in the first column
for jj=1:1:xsize,
mode(jj,1) = Y(jj);
end
%part5--give initial 0 to xorigin and xend
xorigin = X1;
xend = xorigin;
%part6--start to find an IMF-----IMF loop start
nmode = 1;
while nmode <= TNM,
xstart = xend; %last loop value assign to new iteration loop
%xstart -loop start data
iter = 1;
%loop index initial value
%part7--sift 10 times to get IMF---sift loop start
while iter<=10,
spline ,please see part11.
[spmax, spmin, flag]=extrema(xstart); %call function extrema
%the usage of
upper= spline(spmax(:,1),spmax(:,2),dd); %upper spline bound of this sift
lower= spline(spmin(:,1),spmin(:,2),dd); %lower spline bound of this sift
mean_ul = (upper + lower)/2;%spline mean of upper and lower
xstart = xstart - mean_ul;%extract spline mean from Xstart
iter = iter +1;
end
%part7--sift 10 times to get IMF---sift loop end
%part8--subtract IMF from data ,then let the residual xend to start to find next IMF
xend = xend - xstart;
nmode=nmode+1;
%part9--after sift 10 times,that xstart is this time IMF
for jj=1:1:xsize,
mode(jj,nmode) = xstart(jj);
end
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
