rbotiqulei.m
上传用户:cqjiali1
上传日期:2022-07-09
资源大小:78k
文件大小:3k
源码类别:

波变换

开发平台:

Matlab

  1. %--------------------------------------------------------------------------
  2. % 利用小波变换检测R波,并给出模极大曲线
  3. %--------------------------------------------------------------------------
  4. clear all;
  5. web -browser http://www.ilovematlab.cn/thread-11204-1-1.html
  6. load M.mat
  7. ecgdata=M(:,2)';
  8. plot(ecgdata)
  9. points=4000;%截取的点数
  10. range=2:1:32;%尺度范围
  11. sr=360;%抽样率
  12. signal=ecgdata(1:points)';
  13. figure,subplot(311)
  14. plot(signal),axis tight,grid on;
  15. % 对心电信号作小波变换,选'db1'为分析用小波
  16. wf='db1';
  17. subplot(312)
  18. ccfs = cwt(signal,range,wf); 
  19. imagesc(ccfs);
  20. %title('Continuous Transform, coefficients.') 
  21. colormap(gray(128)); 
  22. %ylabel('Scale')
  23. swd=ccfs;
  24. ddw=zeros(size(swd));
  25. % 对ECG信号进行分段
  26. pddw=ddw;
  27. nddw=ddw;
  28. posw=swd.*(swd>0);
  29. pdw=((posw(:,1:points-1)-posw(:,2:points))<0);
  30. pddw(:,2:points-1)=((pdw(:,1:points-2)-pdw(:,2:points-1))>0);
  31. negw=swd.*(swd<0);
  32. ndw=((negw(:,1:points-1)-negw(:,2:points))>0);
  33. nddw(:,2:points-1)=((ndw(:,1:points-2)-ndw(:,2:points-1))>0);
  34. ddw=pddw|nddw;
  35. ddw(:,1)=1;
  36. ddw(:,points)=1;
  37. wpeak=ddw.*swd;
  38. wpeak(:,1)=wpeak(:,1)+1e-10;
  39. wpeak(:,points)=wpeak(:,points)+1e-10;
  40. %figure,imagesc(wpeak);
  41. sam=sum(ccfs,1);
  42. posi=sam.*(sam>0);
  43. posi=(posi>max(posi)/2);
  44. nega=sam.*(sam<0);
  45. nega=-1*(nega<min(nega)/2);
  46. interva=posi+nega;
  47. loca=find(interva);
  48. diff=interva(loca(1:length(loca)-1))-interva(loca(2:length(loca)));
  49. seg=[1,round((loca(find((diff==2)))+loca(find((diff==2))+1))/2),points];%分段的结果保存在seg中
  50. subplot(313)
  51. plot(interva,'r--');
  52. hold on;
  53. plot(signal/100);
  54. hold on
  55. stem(seg,ones(length(seg))),grid on;
  56. % 分别对每一段进行R波检测
  57. for pha=1:length(seg)-1
  58.     block=wpeak(:,[seg(pha):seg(pha+1)]);
  59.     [maxvalue,maxpos]=max(block,[],2);
  60.     coefposi=polyfit(range,maxpos',1);
  61.     [minvalue,minpos]=min(block,[],2);% 找出模极大的位置和值
  62.     coefnega=polyfit(range,minpos',1);% 对模极大点进行直线拟合找到最大线
  63.     Rpeak(pha)=round((coefposi(2)+coefnega(2))/2+seg(pha)-1)/sr;% 每一段R波时刻存在Rpeak中
  64. end  
  65. % 绘制模极大随尺度变化的曲线
  66. figure
  67. subplot(223)
  68. plot(log2(range),log2(abs(minvalue./sqrt(range'))),'b--');hold on;
  69. plot(log2(range),log2(abs(maxvalue./sqrt(range'))),'r');grid on;
  70. %xlabel('log2(s)'),ylabel('log2(maxima)'), 
  71. %title('estimation of the singularity at the R peak'),, hold off;
  72. R_R=Rpeak(2:length(Rpeak))-Rpeak(1:length(Rpeak)-1);
  73. subplot(221)
  74. imagesc(ccfs(:,[seg(pha):seg(pha+1)]));colormap(gray);
  75. block=zeros(size(block));
  76. for i=1:length(range)
  77.     block(i,maxpos(i))=1;
  78.     block(i,minpos(i))=-1;
  79. end
  80. subplot(222)
  81. imagesc(block);colormap(gray);
  82. subplot(224)
  83. plot(maxpos,range,'r+'),hold on,plot(minpos,range,'k*'),grid on;
  84. xr=[0, 1, range];
  85. hold on,plot(xr*(coefnega(1))+coefnega(2),xr,'-');
  86. hold on,plot(xr*(coefposi(1))+coefposi(2),xr,'--');