请选择 进入手机版 | 继续访问电脑版

Hello Mat

 找回密码
 立即注册

QQ登录

只需一步,快速开始

查看: 138|回复: 0

6-EMD端点效应优化(相似特征尺度波形延拓法)(MATLAB视频)

[复制链接]

645

主题

760

帖子

4123

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4123
发表于 2018-9-1 22:19:38 | 显示全部楼层 |阅读模式
6-EMD端点效应优化(相似特征尺度波形延拓法)(MATLAB视频)
百度网盘链接:https://pan.baidu.com/s/11owgsJ42vJ4t-DQRSp9HXQ
具体链接在halcom.cn论坛,联系人QQ:3283892722
该论坛是一个学习交流平台,我会逐一的和大家分享学习。
欢迎大家录制视频,你可在论坛进行打赏分享。
视频专用播放器:http://halcom.cn/forum.php?mod=viewthread&tid=258&extra=page%3D1\
主函数如下:
  1. clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
  2. warning off         % 消除警告
  3. format short        % 数据类型
  4. tic                 % 运算计时
  5. %% 初始化信号
  6. [t, signal] = gen_signal(5e3, 1e3);
  7. figure(1), plot(t,signal,'linewidth',1)
  8. xlabel('t');  ylabel('实际信号');  grid on;
  9. %% 相似特征尺度沿拓
  10. %% 极值点和过零点位置提取
  11. [indmin, indmax, indzer] = extr(signal, t); % 极值点和过零点位置提取
  12. figure(2), plot(t,signal,'linewidth',1);
  13. hold on
  14. plot(t(indmin), signal(indmin), 'm.-')
  15. plot(t(indmax), signal(indmax), 'r.-')
  16. plot(t(indzer), signal(indzer), 'g.-')
  17. xlabel('t');  ylabel('实际信号');  grid on;
  18. legend('原始信号','极小值','极大值','过零点')
  19. axis tight;
  20. hold off;
  21. %% 相似特征尺度沿拓
  22. % 左端点沿拓
  23. [ S_index_z, indmin_index_z, percentage_z ] = find_Si_Z( signal, indmin, indmax, 0 );
  24. [ S_index_f, indmin_index_f, percentage_f ] = find_Si_F( signal, indmin, indmax, 0 );
  25. percentage = [percentage_z; percentage_f];
  26. [indexa, indexb] = find( abs( percentage(:,1) )<0.2*5000 );
  27. [minP, index] = min( abs( percentage(indexa, 7) ) );
  28. indmin_index = percentage(indexa(index), 2);
  29. S_index = percentage(indexa(index), 3);
  30. if(indmax(indmin_index-1)<=S_index)
  31.     newdata = signal( indmax(indmin_index-1):1:S_index );
  32. else
  33.     newdata = signal( indmax(indmin_index-1):-1:S_index );
  34. end
  35. welddata_left = weld( signal, newdata, 0 );

  36. % 右端点沿拓
  37. [ S_index_z, indmin_index_z, percentage_z ] = find_Si_Z( signal, indmin, indmax, 1 );
  38. [ S_index_f, indmin_index_f, percentage_f ] = find_Si_F( signal, indmin, indmax, 1 );
  39. percentage = [percentage_z; percentage_f];
  40. [indexa, indexb] = find( abs( percentage(:,1) )<0.2*5000 );
  41. [minP, index] = min( abs( percentage(indexa, 7) ) );
  42. indmin_index = percentage(indexa(index), 2);
  43. S_index = percentage(indexa(index), 3);
  44. if(indmin(indmin_index+1)<=S_index)
  45.     newdata = signal( S_index : indmin(end) );
  46. else
  47.     newdata = signal( S_index:indmin(indmin_index+1) );
  48. end
  49. welddata_right = weld( signal, newdata, 1 );

  50. welddata = [welddata_left, signal, welddata_right];

  51. figure(3),hold on
  52. plot(1:length(welddata_left), welddata_left,'r.-');
  53. plot(length(welddata_left): length(welddata_left)+length(signal)-1 ,signal,'b-');
  54. plot(length(welddata)-length(welddata_right)-1:length(welddata)-2,welddata_right,'r.-');
  55. xlabel('t');  ylabel('实际信号');  grid on;
  56. axis tight;
  57. hold off;
复制代码
正向匹配:
  1. function [ S_index, indmin_index, percentage ] = find_Si_Z( signal, indmin, indmax, flag )
  2. % 正向沿拓
  3. % 相似三角波沿拓
  4. if(flag==0)   % 左端点沿拓
  5.    % 左端点
  6.    x1 = 1;          y1 = signal(x1);  % s0
  7.    x2 = indmin(1);  y2 = signal(x2);  % n0
  8.    x3 = indmax(1);  y3 = signal(x3);  % m0
  9. %    S1U1 = abs( (x2-x1)^2+(y2-y1)^2 );
  10. %    S1V1 = abs( (x3-x1)^2+(y3-y1)^2 );
  11. %    U1V1 = abs( (x3-x2)^2+(y3-y2)^2 );
  12.    % y=k(x-x1)+y1
  13.    k0 = (y2-y1)./(x2-x1+eps);  % 直线s0-n0
  14.    Q0_x=x3;
  15.    Q0_y = k0.*(Q0_x-x1)+y1;
  16.    h0 = abs( y3-Q0_y );
  17.    r01 = abs( sqrt( (x1-Q0_x)^2+(y1-Q0_y)^2 ) );
  18.    r02 = abs( sqrt( (x2-Q0_x)^2+(y2-Q0_y)^2 ) );
  19.    % 匹配相似波
  20.    t = min( [indmin(1), indmax(1)] );
  21.    percentage = [];
  22.    k = 1;
  23.    for i=2:min( [length(indmin), length(indmax)] )
  24.        x2 = indmin(i);  y2 = signal(x2);  % ni
  25.        x3 = indmax(i);  y3 = signal(x3);  % mi
  26.        for j=indmax(i-1):indmax(i)
  27.            x1 = j;       y1 = signal(x1);  % si
  28. %            SiUi = abs( (x2-x1)^2+(y2-y1)^2 );
  29. %            SiVi = abs( (x3-x1)^2+(y3-y1)^2 );
  30. %            UiVi = abs( (x3-x2)^2+(y3-y2)^2 );
  31.            % y=k(x-x1)+y1
  32.            ki = (y2-y1)./(x2-x1+eps);  % 直线si-ni
  33.            Qi_x=x3;
  34.            Qi_y = ki.*(Qi_x-x1)+y1;
  35.            hi = abs( y3-Qi_y );
  36.            ri1 = abs( sqrt( (x1-Qi_x)^2+(y1-Qi_y)^2 ) );
  37.            ri2 = abs( sqrt( (x2-Qi_x)^2+(y2-Qi_y)^2 ) );
  38.            
  39.            percentage(k, 1) = hi*(r01+r02)-h0*(ri1+ri2);
  40.            percentage(k, 2) = i;
  41.            percentage(k, 3) = j;
  42.            percentage(k, 4) = hi;
  43.            percentage(k, 5) = ri1;
  44.            percentage(k, 6) = ri2;
  45.            percentage(k, 7) = abs( ri1/(ri2+eps) - r01/(r02+eps) );
  46.            percentage(k, 8) = abs( abs(k0)-abs(ki) );
  47.            k=k+1;
  48.        end
  49.    end
  50.    [indexa, indexb] = find( abs( percentage(:,1) )<0.2*5000 );
  51.    [minP, index] = min( abs( percentage(indexa, 7) ) );
  52.    indmin_index = percentage(indexa(index), 2);
  53.    S_index = percentage(indexa(index), 3);
  54. elseif(flag==1)   % 右端点沿拓
  55.    % 右端点
  56.    x1 = length(signal);   y1 = signal(x1);  % s0
  57.    x2 = indmin(end);      y2 = signal(x2);  % n0
  58.    x3 = indmax(end);      y3 = signal(x3);  % m0
  59. %    S1U1 = abs( (x2-x1)^2+(y2-y1)^2 );
  60. %    S1V1 = abs( (x3-x1)^2+(y3-y1)^2 );
  61. %    U1V1 = abs( (x3-x2)^2+(y3-y2)^2 );
  62. %    k1 = (y2-y1)./(x2-x1+eps);
  63. %    k2 = (y3-y1)./(x3-x1+eps);
  64.    % y=k(x-x1)+y1
  65.    k0 = (y2-y1)./(x2-x1+eps);  % 直线s0-n0
  66.    Q0_x=x3;
  67.    Q0_y = k0.*(Q0_x-x1)+y1;
  68.    h0 = abs( y3-Q0_y );
  69.    r01 = abs( sqrt( (x1-Q0_x)^2+(y1-Q0_y)^2 ) );
  70.    r02 = abs( sqrt( (x2-Q0_x)^2+(y2-Q0_y)^2 ) );
  71.    % 匹配相似波
  72.    percentage = [];
  73.    k = 1;
  74.    for i=min( [length(indmin), length(indmax)] )-1 : -1 : 1
  75.        x2 = indmin(i);  y2 = signal(x2);  % ni
  76.        x3 = indmax(i);  y3 = signal(x3);  % mi
  77.        for j=indmax(i):indmin(i+1)
  78.            x1 = j;     y1 = signal(x1);   % si
  79. %            SiUi = abs( (x2-x1)^2+(y2-y1)^2 );
  80. %            SiVi = abs( (x3-x1)^2+(y3-y1)^2 );
  81. %            UiVi = abs( (x3-x2)^2+(y3-y2)^2 );
  82. %            k11 = (y2-y1)./(x2-x1+eps);
  83. %            k22 = (y3-y1)./(x3-x1+eps);
  84.            % y=k(x-x1)+y1
  85.            ki = (y2-y1)./(x2-x1+eps);  % 直线si-ni
  86.            Qi_x=x3;
  87.            Qi_y = ki.*(Qi_x-x1)+y1;
  88.            hi = abs( y3-Qi_y );
  89.            ri1 = abs( sqrt( (x1-Qi_x)^2+(y1-Qi_y)^2 ) );
  90.            ri2 = abs( sqrt( (x2-Qi_x)^2+(y2-Qi_y)^2 ) );
  91.            
  92.            percentage(k, 1) = hi*(r01+r02)-h0*(ri1+ri2);
  93.            percentage(k, 2) = i;
  94.            percentage(k, 3) = j;
  95.            percentage(k, 4) = hi;
  96.            percentage(k, 5) = ri1;
  97.            percentage(k, 6) = ri2;
  98.            percentage(k, 7) = abs( ri1/(ri2+eps) - r01/(r02+eps) );
  99.            percentage(k, 8) = abs( abs(k0)-abs(ki) );
  100.            k=k+1;
  101.        end
  102.    end
  103.    [indexa, indexb] = find( abs( percentage(:,1) )<0.2*5000 );
  104.    [minP, index] = min( abs( percentage(indexa, 7) ) );
  105.    indmin_index = percentage(indexa(index), 2);
  106.    S_index = percentage(indexa(index), 3);
  107. end
复制代码

反向匹配:
  1. function [ S_index, indmin_index, percentage ] = find_Si_F( signal, indmin, indmax, flag )
  2. % 反向沿拓
  3. % 相似三角波沿拓
  4. if(flag==0)   % 左端点沿拓
  5.    % 左端点
  6.    x1 = 1;          y1 = signal(x1);  % s0
  7.    x2 = indmin(1);  y2 = signal(x2);  % n0
  8.    x3 = indmax(1);  y3 = signal(x3);  % m0
  9. %    S1U1 = abs( (x2-x1)^2+(y2-y1)^2 );
  10. %    S1V1 = abs( (x3-x1)^2+(y3-y1)^2 );
  11. %    U1V1 = abs( (x3-x2)^2+(y3-y2)^2 );
  12.    % y=k(x-x1)+y1
  13.    k0 = (y2-y1)./(x2-x1+eps);  % 直线s0-n0
  14.    Q0_x=x3;
  15.    Q0_y = k0.*(Q0_x-x1)+y1;
  16.    h0 = abs( y3-Q0_y );
  17.    r01 = abs( sqrt( (x1-Q0_x)^2+(y1-Q0_y)^2 ) );
  18.    r02 = abs( sqrt( (x2-Q0_x)^2+(y2-Q0_y)^2 ) );
  19.    % 匹配相似波
  20.    t = min( [indmin(1), indmax(1)] );
  21.    percentage = [];
  22.    k = 1;
  23.    for i=2:min( [length(indmin), length(indmax)] )
  24.        x2 = indmin(i-1);  y2 = signal(x2);  % ni
  25.        x3 = indmax(i);  y3 = signal(x3);  % mi
  26.        for j=indmax(i-1):indmax(i)
  27.            x1 = j;       y1 = signal(x1);  % si
  28. %            SiUi = abs( (x2-x1)^2+(y2-y1)^2 );
  29. %            SiVi = abs( (x3-x1)^2+(y3-y1)^2 );
  30. %            UiVi = abs( (x3-x2)^2+(y3-y2)^2 );
  31.            % y=k(x-x1)+y1
  32.            ki = (y2-y1)./(x2-x1+eps);  % 直线si-ni
  33.            Qi_x=x3;
  34.            Qi_y = ki.*(Qi_x-x1)+y1;
  35.            hi = abs( y3-Qi_y );
  36.            ri1 = abs( sqrt( (x1-Qi_x)^2+(y1-Qi_y)^2 ) );
  37.            ri2 = abs( sqrt( (x2-Qi_x)^2+(y2-Qi_y)^2 ) );
  38.            
  39.            percentage(k, 1) = hi*(r01+r02)-h0*(ri1+ri2);
  40.            percentage(k, 2) = i;
  41.            percentage(k, 3) = j;
  42.            percentage(k, 4) = hi;
  43.            percentage(k, 5) = ri1;
  44.            percentage(k, 6) = ri2;
  45.            percentage(k, 7) = abs( ri1/(ri2+eps) - r01/(r02+eps) );
  46.            percentage(k, 8) = abs( abs(k0)-abs(ki) );
  47.            k=k+1;
  48.        end
  49.    end
  50.    [indexa, indexb] = find( abs( percentage(:,1) )<0.2*5000 );
  51.    [minP, index] = min( abs( percentage(indexa, 7) ) );
  52.    indmin_index = percentage(indexa(index), 2);
  53.    S_index = percentage(indexa(index), 3);
  54. elseif(flag==1)   % 右端点沿拓
  55.    % 右端点
  56.    x1 = length(signal);   y1 = signal(x1);  % s0
  57.    x2 = indmin(end);      y2 = signal(x2);  % n0
  58.    x3 = indmax(end);      y3 = signal(x3);  % m0
  59. %    S1U1 = abs( (x2-x1)^2+(y2-y1)^2 );
  60. %    S1V1 = abs( (x3-x1)^2+(y3-y1)^2 );
  61. %    U1V1 = abs( (x3-x2)^2+(y3-y2)^2 );
  62. %    k1 = (y2-y1)./(x2-x1+eps);
  63. %    k2 = (y3-y1)./(x3-x1+eps);
  64.    % y=k(x-x1)+y1
  65.    k0 = (y2-y1)./(x2-x1+eps);  % 直线s0-n0
  66.    Q0_x=x3;
  67.    Q0_y = k0.*(Q0_x-x1)+y1;
  68.    h0 = abs( y3-Q0_y );
  69.    r01 = abs( sqrt( (x1-Q0_x)^2+(y1-Q0_y)^2 ) );
  70.    r02 = abs( sqrt( (x2-Q0_x)^2+(y2-Q0_y)^2 ) );
  71.    % 匹配相似波
  72.    percentage = [];
  73.    k = 1;
  74.    for i=min( [length(indmin), length(indmax)] )-1 : -1 : 1
  75.        x2 = indmin(i+1);  y2 = signal(x2);  % ni
  76.        x3 = indmax(i);  y3 = signal(x3);  % mi
  77.        for j=indmax(i):indmin(i+1)
  78.            x1 = j;     y1 = signal(x1);   % si
  79. %            SiUi = abs( (x2-x1)^2+(y2-y1)^2 );
  80. %            SiVi = abs( (x3-x1)^2+(y3-y1)^2 );
  81. %            UiVi = abs( (x3-x2)^2+(y3-y2)^2 );
  82. %            k11 = (y2-y1)./(x2-x1+eps);
  83. %            k22 = (y3-y1)./(x3-x1+eps);
  84.            % y=k(x-x1)+y1
  85.            ki = (y2-y1)./(x2-x1+eps);  % 直线si-ni
  86.            Qi_x=x3;
  87.            Qi_y = ki.*(Qi_x-x1)+y1;
  88.            hi = abs( y3-Qi_y );
  89.            ri1 = abs( sqrt( (x1-Qi_x)^2+(y1-Qi_y)^2 ) );
  90.            ri2 = abs( sqrt( (x2-Qi_x)^2+(y2-Qi_y)^2 ) );
  91.            
  92.            percentage(k, 1) = hi*(r01+r02)-h0*(ri1+ri2);
  93.            percentage(k, 2) = i;
  94.            percentage(k, 3) = j;
  95.            percentage(k, 4) = hi;
  96.            percentage(k, 5) = ri1;
  97.            percentage(k, 6) = ri2;
  98.            percentage(k, 7) = abs( ri1/(ri2+eps) - r01/(r02+eps) );
  99.            percentage(k, 8) = abs( abs(k0)-abs(ki) );
  100.            k=k+1;
  101.        end
  102.    end
  103.    [indexa, indexb] = find( abs( percentage(:,1) )<0.2*5000 );
  104.    [minP, index] = min( abs( percentage(indexa, 7) ) );
  105.    indmin_index = percentage(indexa(index), 2);
  106.    S_index = percentage(indexa(index), 3);
  107. end
复制代码







本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
程序:算法QQ  3283892722
群智能算法视频售价:38元,链接http://halcom.cn/forum.php?mod=forumdisplay&fid=73
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则


Python|Opencv|MATLAB|Halcom.cn  

GMT+8, 2018-12-15 17:00 , Processed in 0.139478 second(s), 36 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.

快速回复 返回顶部 返回列表