首页 > 编程语言 >基础版EMD分解函数Matlab程序

基础版EMD分解函数Matlab程序

时间:2023-06-07 19:44:08浏览次数:48  
标签:EMD end imf findpeaks length 分解 Matlab x1

不调用matlab自带emd(x)函数,将其内容简化为如下部分
EMD分解基础步骤可以参见:[意念回复:经验模态分解(EMD)](https://blog.csdn.net/weixin_39910711/article/details/124661068?spm=1001.2014.3001.5506)
原始程序为百度搜索,结合ChatGPT后给出相应注释。

% EMD分解程序
% 作者:CSDN:Mist_n  博客园:Drizzly_n
% 日期:2023.06.07
% 注释:调用子程序EMD分解基本结构(原始文件来源未知,根据内容添加注释)
% 参考:https://blog.csdn.net/weixin_39910711/article/details/124661068?spm=1001.2014.3001.5506

clc
clear all

% 定义输入信号
Ts = (1/512);
Fs = 1/Ts;
t=0:Ts:1;
v1 = 0.5*cos(2*20*pi*t);
v2 = cos(2*30*pi*t);
input = v1+v2;
% x = awgn(input,1);
x=input;

%绘制出入信号时域图
figure(1)
plot(t,x);hold on;
grid on
xlabel('t');ylabel('x');title('初始信号');
legend('0.5cos(40x*pi)+cos(60x*pi)');

%求取本证模态函数imf并绘制其时图像及频谱
imf = emd(x);

[M, N]=size(imf);
figure(2);
subplot(M+1,1,1);
plot(x);hold on;
ylabel('x')
set(gca,'xtick',[])
axis tight;
title('合成信号x做EMD分解');

for i=2:M
    subplot(M+1,1,i);
    plot(imf(i-1,:),'r');hold on;xlim([100 400]);
    ylabel (['IMF ' num2str(i-1)]);
    set(gca,'xtick',[])
    xlim([1 length(x)])
end

subplot(M+1,1,M+1)
plot(imf(M,:));hold on;
ylabel(['RS ' num2str(M)])
xlim([1 length(x)])

%计算相关系数并绘制
 for i = 1:1:M
     cc(i)=min(min(corrcoef(imf(i,:),x)));
 end
figure(3)
plot(cc,'-g<','LineWidth',1.5,'MarkerEdgeColor','b','MarkerFaceColor','b','MarkerSize',5);
set(gca,'XGrid', 'on', 'YGrid', 'on');
legend('CC');
xlabel('IMF');
ylabel('相关系数');

% 重构信号
cg_ev=imf(1,:)+imf(2,:);
figure(4);
plot(t,x);
hold on
plot(t,cg_ev,'r');
xlabel('t/s');ylabel('幅值');legend('原信号','重构信号');
%


function imf = emd(x)

x= transpose(x(:));%转置为行矩阵
imf = [];

while ~ismonotonic(x) %当x不是单调函数,分解终止条件
    x1 = x;
    sd = Inf;%均值

    %直到x1满足IMF条件,得c1
    % 两个条件:)
    %(1)在任意时刻,由局部极大值点形成的上包络线和由局部极小值点形成的下包络线的平均值为零,
    %     即上、下包络线相对于时间轴局部对称
    %(2)在整个时程内极值点个数与过零点个数相等或最多相差1

    while (sd > 0.1) || ~isimf(x1) %当标准偏差系数sd大于0.1或x1不是固有模态函数时,分量终止条件

        s1 = getspline(x1);%上包络线
        s2 = -getspline(-x1);%下包络线

        x2 = x1-(s1+s2)/2;%此处的x2为中间信号,即为链接中的C1,1,(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;% 余量
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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 = 1;
end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function u = isimf(x)
%u=0表示x不是固有模式函数,u=1表示x是固有模式函数
N  = length(x);
u1 = sum(x(1:N-1).*x(2:N) < 0);% 确定过零点个数:受限判断相邻两个点的成绩,小于0则表示一正一副,置1,反之置0.
u2 = length(findpeaks(x))+length(findpeaks(-x));% 确定极值点个数
if abs(u1-u2) > 1 % 在整个时程内极值点个数与过零点个数相等或最多相差1
    u = 0;
else
    u = 1;
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function s = getspline(x)
%三次样条函数拟合成元数据包络线
N = length(x);
p = findpeaks(x);
s = spline([0 p N+1],[0 x(p) 0],1:N);%spline(x,y,xq),x 坐标; y - x 坐标处的函数值; xq - 查询点
% y处前后取0的作用是使包络线在起始点与末尾处不发散。
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function n = findpeaks(x)
% 寻找极值点
n    = find(diff(diff(x) > 0) < 0);
u    = find(x(n+1) > x(n));
n(u) = n(u)+1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 

标签:EMD,end,imf,findpeaks,length,分解,Matlab,x1
From: https://www.cnblogs.com/Nigel-ff/p/17464375.html

相关文章

  • Matlab GUI编程入门
    先吐槽下,我看得是出版社的优秀教材,虽说是二等奖,《Matlab与控制系统仿真实践》第3版,北京航空航天大学出版社赵广元编著,想着入门,找本国内的书快速入门。谁知道一个上午就这么过去了。以一个控制系统典型环节演示为例,课本从P95开始:一、准备工作,设置下显示组件的名称与图标HOME-......
  • MATLAB偏最小二乘回归(PLSR)和主成分回归(PCR)分析光谱数据|附代码数据
    全文链接:http://tecdat.cn/?p=2655最近我们被客户要求撰写关于偏最小二乘回归(PLSR)和主成分回归(PCR)的研究报告,包括一些图形和统计输出。此示例显示如何在matlab中应用偏最小二乘回归(PLSR)和主成分回归(PCR),并讨论这两种方法的有效性当存在大量预测变量时,PLSR和PCR都是对因变量建模......
  • 基于mfcc和DTW语音信息特征提取算法matlab仿真
    1.算法仿真效果matlab2022a仿真结果如下:2.算法涉及理论知识概要在语音识别(SpeechRecognition)和话者识别(SpeakerRecognition)方面,最常用到的语音特征就是梅尔倒谱系数(Mel-scaleFrequencyCepstralCoefficients,简称MFCC)。根据人耳听觉机理的研究发现,人耳对不同频率的声波有不......
  • 基于mfcc和DTW语音信息特征提取算法matlab仿真
    1.算法仿真效果matlab2022a仿真结果如下:   2.算法涉及理论知识概要       在语音识别(SpeechRecognition)和话者识别(SpeakerRecognition)方面,最常用到的语音特征就是梅尔倒谱系数(Mel-scaleFrequencyCepstralCoefficients,简称MFCC)。根据人耳听觉机理的研究发......
  • 基于Matlab的毫米波雷达多路径核心代码
    ✅作者简介:热爱科研的Matlab仿真开发者,修心和技术同步精进,matlab项目合作可私信。......
  • 【反演】基于遗传算法实现均匀地层模型随钻电磁波测井反演附matlab代码
    ✅作者简介:热爱科研的Matlab仿真开发者,修心和技术同步精进,matlab项目合作可私信。......
  • matlab 更改中文语言,但是常规设置没有桌面语言的问题
    没有如下图红框的选项的原因是,操作系统的语言设置没有全部设置成中文造成的。尤其是时间和语言中的区域设置中的区域格式。解决方法:设置->区域 ......
  • matlab r2022b 在vbox下显示plot上下镜像的问题
    现象:matlabr2022b在vbox环境下plot显示会出现上下镜像的现象。原因:默认开启的OpenGL渲染模式造成的。解决:关闭OpenGL渲染模式或者代码改用"painters"渲染模式。set(gcf,'renderer','painters'); ......
  • m基于BBV网络的节点强度分布算法matlab仿真
    1.算法仿真效果matlab2022a仿真结果如下:   2.算法涉及理论知识概要       随着互联网的发展和数据规模的不断增大,网络科学在各个领域中得到了广泛应用。在网络科学中,节点强度是一个重要的指标,它用于描述一个节点在网络中的重要性或中心性。本文提出了一种基于......
  • m基于BBV网络的节点强度分布算法matlab仿真
    1.算法仿真效果matlab2022a仿真结果如下:2.算法涉及理论知识概要随着互联网的发展和数据规模的不断增大,网络科学在各个领域中得到了广泛应用。在网络科学中,节点强度是一个重要的指标,它用于描述一个节点在网络中的重要性或中心性。本文提出了一种基于BBV网络的节点强度分布算法......