波形包络提取与峰值提取

在工作中经常会对信号极值点进行提取，这里就需要用到信号波形包络提取和峰值提取方法。翻阅了一些资料，下面简单总结。

1、包络检测又称幅度解调，当调制信号已知时，常用同步解调的方式，又称“相敏解调”。当载波信号未知时，可以采用包络检波法，又称包络解调，适用于普通调幅信号的解调，指产生的输出信号与已调信号包络线成正比的幅度解调。包络检波原理如图

当 Ui(t) < Uo(t-) 时，RC * dUo/dt = Uo

化成差分方程为：

代码如下：

/**
 * 包络检波，模拟了硬件全波检波的过程
 * rc = 0 时初始化

 *参数 rct，对应的是硬件电路中的RC时间常数，要根据待检测的包络信号的频带来确定

 **/
double env_3(double x, double rct)
{
    static double old_y = 0.0;
    if(rct == 0.0)
    {
        old_y = 0.0;
    }
    else
    {
        x = fabs(x);
        if(x > old_y)
        {
            old_y = x;
        }
        else
        {
            old_y *= rct / ( rct + 1 );
        }
    }
    return old_y;
}
void env_4(double x[], double y[], int N, double rct)
{
    double xx = 0.0;
    int i;
    y[0] = fabs(x[0]);
    for(i = 1; i < N; i++)
    {
        xx = fabs(x[i]);
        if( xx > y[i-1])
        {
            y[i] = xx;
        }
        else
        {
            y[i] = y[i-1] * rct / ( rct + 1 );
        }
    }
}

代码效果如图

  2、峰值提取

得到包络曲线后，常需要对峰值进行提取。一个是逐级扫描，扫描的依据就是斜率的变化，好处是绝对值比较，可以避免极值在平滑处 的影响；还有一个就是抓住问题精髓，将斜率变化，转化为求导 处理，缺点是在峰不明锐的时候，同一点会有多个值，需要对数据的幅度放大，滤波。

逐级斜率扫描matlab代码如下

function [t,px pk vx valley]=exseek(curv,wd,th)
%%输入一维数组和滤波宽度wd(推荐值50)，输出峰值pk和位置px
a=zeros(size(curv));
a(2:end)=diff(curv);
n=length(a)-1;
k=1;j=1;
px=[];pk=[];vx=[];valley=[];t=[];
b=zeros(3*n+1,1);b(1:n)=a(1:end-1);b(n+1:2*n+1)=a;b(2*n+2:end)=a(2:end);

for i=(wd+1):(3*n+1-wd)
    temp=b((i-wd):(i+wd));
    if b(i)==max(temp)&&b(i)>th
        pk(k)=b(i);
        px(k)=i;
        k=k+1;
    end
    if b(i)==min(temp)&&b(i)<-th
        valley(j)=b(i);
        vx(j)=i-1;
        j=j+1;
    end
end

if ~isempty(px)
    if length(px)==length(vx)
        if px(1)<vx(1)
            t=px/2+vx/2;
        elseif px(1)>vx(1)
            t=px(1:end-1)/2+vx(2:end)/2;
        else
            msgbox('px(1)与vx(1)重合，请检查')
        end
    else
        msgbox('px与vx数目不等，请检查')
    end
end
t=t-n-1;
t(t<n)=[];t(t>2*n+1)=[];
plot(a);hold on
if ~isempty(px)
    for j=1:length(pk)
        plot(px(j),pk(j),'or')
    end
end
if ~isempty(vx)
    for j=1:length(vx)
        plot(vx(j),valley(j),'*g')
        plot(t(j),0,'.r')
    end
end
hold off

第二种方法matlab代码如下

x = 0:.1:4*pi;
y = @(x) sin(x)./cos(x);

y0 = y(x);
yy1 = diff(y0);
yy1 = sign(yy1);
yy1 = diff(yy1);
f = find(yy1<0)+1; % 峰
g = find(yy1>0)+1; % 谷

hold on;
plot(x,y0);
plot(x(f),y(x(f)),'ro');
plot(x(g),y(x(g)),'go');
hold off;
grid on;