经验模态分解python_EMD经验模态分解
今天看了些EMD信号分解方面的东西,matlab官网上有个Hilbert-Huang
Transform的代码,代码效率极高啊,人家3句语句就解决了一个大问题,很牛啊!还有一个GRilling的EMD工具箱,好多文件,功能应该相当强大。
这里研究了研究matlab官网的代码,加了些注释、功能演示,效果如下
原始信号由3个正弦信号加噪声组成,如下
下面为做EMD分解的结果
第三次分解信号的瞬时频率如下
第四次分解信号的Hilbert分析
具体代码如下
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 =
find(diff(diff(x) > 0) < 0); %
相当于找二阶导小于0的点
u =
find(x(n+1) > x(n));
n(u) =
n(u)+1; % 加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 = x-x1;
end
imf{end+1} = x;
% 是否单调
function u = ismonotonic(x)
u1 = length(findpeaks(x))*length(findpeaks(-x));
if u1 > 0
u = 0;
else
u = 1;
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;
else
u = 1;
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
经验模态分解python_EMD经验模态分解相关推荐
- 小波分解与小波包分解代码_分解的功能参数和代码可维护性
小波分解与小波包分解代码 Code keeps changing, there's no doubt about that. We always do our best to set some roc ...
- 前端学习总结,经验分享,项目经验分享过程
前言 来,和魔王哪吒总结一下,分享一下自己对前端学习总结,经验分享,以及写过的项目经验分享过程. 如果觉得还不错的话,浏览的过程中,需要您: 点赞,分享,评论 有钱的捧个钱场,没钱的捧个人场 技术实践 ...
- c#中的模态对话框和非模态对话框
模态对话框 弹出窗口阻止调用窗口的所有消息响应. 只有在弹出窗口结束后调用窗口才能继续. 在模态窗口"关闭"后,可以读取模态窗口中信息,包括窗口的返回状态,窗口子控件的值. 非模态 ...
- MFC模态窗口与非模态窗口
MFC模态窗口与非模态窗口 开发工具与关键技术:C++.VisualStudio 作者:何任贤 撰写时间:2019年07月25日 模态窗口的意思是指主窗口在打开模态窗口后,没法再操作主窗口,这就是模态 ...
- 多模态理论张德禄_结构动力学中的模态分析(3) —— 模态参数及实验模态分析...
引言 前面的文章介绍了模态相关的数学基础及实模态分析. 蒙特遇见卡罗:结构动力学中的模态分析(1) -- 线性系统和频响函数zhuanlan.zhihu.com 蒙特遇见卡罗:结构动力学中的模态分析 ...
- 模态对话框和非模态对话框的消息循环分析
1.非模态对话框和父窗口共享当前线程的消息循环 2.模态对话框新建一个新的消息循环,并由当前消息循环派发消息,而父窗口.模态对话框屏蔽了用户对它父窗口的操作,但是不是在消息循环里面屏蔽,所以给父窗口发 ...
- QT中的模态对话框及非模态对话框
模态对话框(Modal Dialog)与非模态对话框(Modeless Dialog)的概念不是Qt所独有的,在各种不同的平台下都存在.又有叫法是称为模式对话框,无模式对话框等.所谓模态对话框就是在其 ...
- Qt中的模态对话框和非模态对话框
模态对话框及非模态对话框(详情见课本P51). 模态对话框:在没有被关闭之前,用户不能与同一个应用程序的其他窗口进行交互,直到该对话框关闭. 非模态对话框:当被打开时,用户既可选择和该对话框进行交互, ...
- 模态窗口和非模态窗口
转载自:https://my.oschina.net/u/2425942/blog/882879 模态窗口就是在该窗口关闭之前,其父窗口不可能成为活动窗口的那种窗口. 例如: 窗口A弹出窗口B,如果窗 ...
- 小甲鱼 OllyDbg 教程系列 (十四) : 模态对话框 和 非模态对话框 之 URlegal 和 movgear
小甲鱼 OD 使用教程:https://www.bilibili.com/video/av6889190?p=22 exeScope 下载:https://pan.baidu.com/s/1dSWap ...
最新文章
- linux 解压 目录,linux 中目录、文件的解压缩
- Chrome OS 云里雾里
- java独步寻花,小班语言《江畔独步寻花》
- Java 之单元测试
- Linux命令行打开不了发行光盘RHEL_6.3 i386 Disc 1
- Java中使用正则表达式校验字符串
- C语言课程设计:学生管理系统
- gimp 抠图_gimp软件如何实现抠图?
- 【Labplus 3】Scratch获取角色造型的数量
- 常见的软件测试方法有,常见的几种软件测试方法都有哪些
- 麦克林排名计算机,麦克林9大热门大学专业院校排名出炉!启德为您解读
- linux 进程 线程 优先级,Linux编程-线程优先级的设定
- specification java_Java Specification类代码示例
- 解决浏览器被好123劫持主页的问题。
- css实现简单的电影院选座功能
- MATLAB技术沙龙之如何批量处理图像的大小
- 稍微挖掘一下思维导图XMind潜力以及那个使用XMind的你
- 《天道》之丁元英的30句经典语录
- 海康摄像头实时显示与字符叠加详解
- 04【结构面】 面试之结构面,什么是结构面?结构面的准备,常见问题分析?