一种有效的包络线算法

happy

1. function [up,down] = envelope(x,y,interpMethod)
   
2. 
   
3. %ENVELOPE gets the data of upper and down envelope of the known input (x,y).
   
4. %
   
5. %   Input parameters:
   
6. %    x               the abscissa of the given data
   
7. %    y               the ordinate of the given data
   
8. %    interpMethod    the interpolation method
   
9. %
   
10. %   Output parameters:
    
11. %    up      the upper envelope, which has the same length as x.
    
12. %    down    the down envelope, which has the same length as x.
    
13. %
    
14. %   See also DIFF INTERP1
    
15. 
    
16. %   Designed by: Lei Wang, <[email]WangLeiBox@hotmail.com[/email]>, 11-Mar-2003.
    
17. %   Last Revision: 21-Mar-2003.
    
18. %   Dept. Mechanical & Aerospace Engineering, NC State University.
    
19. % $Revision: 1.1 $  $Date: 3/21/2003 10:33 AM $
    
20. 
    
21. if length(x) ~= length(y)
    
22.     error('Two input data should have the same length.');
    
23. end
    
24. 
    
25. if (nargin < 2)|(nargin > 3), error('Please see help for INPUT DATA.');
    
26. elseif (nargin == 2),  interpMethod = 'linear'; end
    
27.    
    
28. % Find the extreme maxim values and the corresponding indexes
    
29. %----------------------------------------------------
    
30. extrMaxValue = y(find(diff(sign(diff(y)))==-2)+1);
    
31. extrMaxIndex =   find(diff(sign(diff(y)))==-2)+1;
    
32. 
    
33. % Find the extreme minim values and the corresponding indexes
    
34. %----------------------------------------------------
    
35. extrMinValue = y(find(diff(sign(diff(y)))==+2)+1);
    
36. extrMinIndex =   find(diff(sign(diff(y)))==+2)+1;
    
37. 
    
38. up = extrMaxValue; up_x = x(extrMaxIndex);
    
39. down = extrMinValue; down_x = x(extrMinIndex);
    
40. 
    
41. % Interpolation of the upper/down envelope data
    
42. %----------------------------------------------------
    
43. up = interp1(up_x,up,x,interpMethod);
    
44. down = interp1(down_x,down,x,interpMethod);

复制代码

happy

使用例子：

1. %DEMOENVELOPE shows how to use function envelope to obtain the
   
2. %   upper/down envelope of a given data and plot the envelope.
   
3. %
   
4. %   See also EVELOPE
   
5. 
   
6. %   Designed by: Lei Wang, <[email]WangLeiBox@hotmail.com[/email]>, 11-Mar-2003.
   
7. %   Last Revision: 21-Mar-2003.
   
8. %   Dept. Mechanical & Aerospace Engineering, NC State University.
   
9. % $Revision: 1.1 $  $Date: 3/21/2003 10:38 AM $
   
10. 
    
11. clc;
    
12. 
    
13. % Load a signal waveform
    
14. %--------------------------------------------
    
15. load data.txt data;
    
16. t = data(:,1); % time series
    
17. y = data(:,2); % signal data
    
18. figure(1); plot(t,y,'b-');  title('The original signal waveform','FontSize',18);
    
19. 
    
20. % Call function envelope to obtain the envelope data
    
21. %--------------------------------------------
    
22. [up,down] = envelope(t,y,'linear');
    
23. 
    
24. % Show the envelope alone
    
25. %--------------------------------------------
    
26. figure(2); plot(t,up); hold on;
    
27. plot(t,down); title('The envelope of the given signal data','FontSize',18); hold off;
    
28. 
    
29. % Show the original signal and its envelope
    
30. %--------------------------------------------
    
31. figure(3); plot(t,y,'g-'); hold on; plot(t,up,'r-.'); plot(t,down,'r-.');
    
32. title('The envelope vs the given signal data','FontSize',18); hold off;

复制代码

qiuqia17

interpMethod 包含几种?
  除了 linear

eight

“help interp1”, and you will find it

[ 本帖最后由 ChaChing 于 2010-5-3 00:33 编辑 ]

qiuqia17

我贴
'nearest'Nearest neighbor interpolation
'linear'Linear interpolation (default)
'spline'Cubic spline interpolation
'pchip'Piecewise cubic Hermite interpolation
'cubic'(Same as 'pchip')
'v5cubic'Cubic interpolation used in MATLAB 5. This method does not extrapolate. Also, if x is not equally spaced, 'spline' is used

junzifei

采用hilbert变换也可以
abs(hilbert(x));就可以
不过问题的重点好像不在这里

2365215

如果x不是单调的怎么办？
我用extrMaxValue = y(find(diff(sign(diff(y)))==-2)+1);
extrMaxIndex = find(diff(sign(diff(y)))==-2)+1;
提取了一下y值，然后再和x一一对应画出的图不是外轮廓线啊

我把图传了一下，请帮我看看。我试着用下轮廓那个做了一下，得到的是一样的图。

上轮廓
extrMaxValue = U2(find(diff(sign(diff(U2(1,:)')))==+2)+1);
extrMaxIndex = find(diff(sign(diff(U2(1,:)')))==+2)+1;

下轮廓
extrMaxValue = U2(find(diff(sign(diff(U2(2,:)')))==-2)+1);
extrMaxIndex = find(diff(sign(diff(U2(2,:)')))==-2)+1;

U2是一个二行1840列的数据

我是这样画图的
for i=1:1840
for j=1:245
if i == extrMaxIndex(j)
plot(Stick(i),extrMaxIndex(j),'.')
end
end
end

sogooda

那就先排序，再做嘛。
doc sort

2365215

谢谢，是我没看懂才会出现问题。 已经搞定。

好像用smooth函数也可以达到类似的效果。(5/12)

ChaChing

首先感谢 Happy教授 !
使用find(diff(sign(diff(y)))==-2)求取极大值的确绝妙
但使用上若有重复点即发生一些小问题如附图
各位大大对如何解决是否有好意见?

zsqlp

你的包络线太精确了，呵呵，有没有一种模糊的啊，画包络线一般就是为了好看，比如区分波包啊之类的，你的有点过头了啊。呵呵

sspitty

matlab好像有个自带的啊 效果看上去还要好~~有没有那代码？？？？

impulse

matlab很容易实现，但一般在做包络时先要对信号进行滤波，最近研究LabVIEW声学振动包时才知道的。

ChaChing

回复 sspitty 的帖子
"matlab好像有个自带的", 那个函数?:@)

回复 junzifei 的帖子
"abs(hilbert(x));就可以", 建议说清楚些! :@)   

禾木目

这程序如果是曲线上升的，得到的包络线就不对了呀