%
一个求包络线和包络谱的程序
%
现代希尔伯特变换解调分析：
%
带通滤波；希尔伯特变换获得信号时域的包络线；用
fft
变换获得包络谱
%
如何获得包络线？
%
信号经希尔伯特变换不能直接得到包络，设信号
x
的希尔伯特变换为
y
，则平方和
%x.^2+y.^2(
或者再开根号，直接取平方和的效果为好
)
才是信号
x
的包络。
%
构造实验数据
clear all;close all;
t=0:0.005:1*pi;
fs=10000;
s=4*sin(2*200*pi*t).*(sin(2*4500*pi*t))+25*(sin(2*4500*pi*t));
figure(1);
subplot(211);plot(t,s);title('
原始信号
');
%
运用小波方法滤波
[c,l]=wavedec(s,1,'db10');
d1=wrcoef('d',c,l,'db10');
a1=0;
subplot(212);plot(d1);title('
滤波后重构的高频信号
');
%
希尔伯特变换求包络线
y=hilbert(d1);
y1=abs(d1+y*j);               %
这是取得包络线的三种方程。看一看哪种效果好。
%y1=abs(y);                     %
或者
z=x.^2+y.^2;
有的取得是
abs(y),
但是不推荐用。
%y1=d1.^2+y.^2;                 %
通过分析，
该方程在包络谱中的效果最好，
即取二者平
方和。
figure(2);
subplot(211);
hold on
plot(t,s);
plot(t,y1,'r');title('
包络线
');
hold off
%FFT
求包络谱
N=1024;
p=abs(fft(y1,N));
subplot(212);
plot((0:N/2-1)/N*fs,p(1:N/2));%
只需取到半频，即
fs/2
%f=(0:N-1)*fs/N; plot(f,p);   %
横坐标是在
fs
上，其中以
fs/2
为轴中心对称。
title('
包络谱
');xlabel('
频率
');ylabel('
功率谱
');
%
对比信号直接的傅里叶变换功率谱与包络谱

如图，对一个边长为1的单位正三角形，在初始时刻，三个动点各自位于三角形的三个顶点处。随后同时出发，以相同的速度沿着三角形的边运动。三个动点的连线在运动的过程中始终形成一个新的，小一点的正三角形。效果如下：
求中间空白区域的面积。
分析
直观上来看，中间空白区域很像一个莱洛三角形(如下图深色部分)，但稍微分析一下，就知道这一定不是莱洛三角形：
莱洛三角形的三条边是三条对应60°圆周角的圆弧，而围出空白区域的三条曲边，我们不妨取其中的一条看看。如下，
这明显不是圆弧，看上去更像一条抛物线。
事实上，如果学过微积分或微分几何，很容易通过解析几何的方法计算曲线族的包络:
建立上图所示的坐标系，AB=BC=AC=1。点E和点F是两个动点。满足AF+AE=1。设AF=t，则
于是直线EF的点斜式方程为
令
我们就得到了一组含有参数t的直线族
它的包络曲线的方程为
消去参数t,可得
这是一个开口向下，以y轴为对称轴，(0，√3/4)为顶点，且经过点(-1/2，0)和(1/2，0)的 抛物线！
根据对称性，另外两条包络曲线也是形状相同的抛物线。
既然要计算面积，就需要知道他们交点的坐标。
我们不必再去计算另外两个抛物线的方程，根据如下对称性，
蓝色抛物线和红色抛物线的交点，也是蓝色抛物线和正三角形一条高的交点。
这条高的方程很容易求，从而和抛物线的交点为(-1/6,2√3/9)。过程略。
连接三个交点，这样待求面积转化为一个边长为1/3的正三角形和三个小曲边三角形的面积。
即
大功告成。
过程稍显复杂。难点主要在于包络曲线的计算。
如果已经知道是抛物线，那么直接根据抛物线经过的三个特殊点(也就是正三角形底边的两个顶点B和C，底边上中位线的中点D)得到其方程，就不必求偏导联立方程组计算了。
怎样确定这是一条抛物线呢？
注意到求出的抛物线恰好以三角形的中心为焦点，再联系抛物线的光学几何性质，不难得到其切线恰好将三角形的周长按1:2分成两部分。根据同一法立得，包络曲线正是这条抛物线。
最后，将题目稍作改编，难度就比较小，且适合高中生做：
对上图所示正三角形，一条抛物线经过其两个顶点。抛物线以该三角形的中心为焦点，以该三角形的高为对称轴。取抛物线上位于在三角形内部的任意一点，做切线DE。
证明：线段DE将正三角形的周长按固定比1:2分成两部分。
这个问题留给你们。

包络线去除是光谱识别中一个重要的步骤，算法如下：
    （1）通过求导得到光谱曲线上所有极大值点，即“峰”值点，然后比较大小，得到极大值点中的最大值点；
    （2）以最大值点作为包络线的一个端点，计算该点与长波方向（波长增长的方向）各个极大值连线的斜率，以斜率最大点作为包络线下一个端点，再以此点为起点循环，直至最后一点；
    （3）以最大值点作为包选线的一个端点，向短波（波长减少的方向）进行类似计算，以斜率最小点为下个端点，再以此点为起点循环，直至曲线上的开始点；
    （4）沿长波方向连接所有端点，可形成曲线的包络线，用实际光谱反射率去除包络线上相应波段的反射率值，可得到包络线消除法归一化后的值。
/**
    求取包络线，去除包络线(为了保存所有去除包络线数据)--->为了提高速度
    r:读取的光谱数据
    w:读取的波长
    h:输出的去包络线数据
*/

template <typename T> void RsEnvelopingLine(T* r, float* w, float* h, int bands)
{
    int i = 0, j = 1, m;
    int k = 0;
    float tempHH = 0;
    while (i < bands - 1)
    {
        if (j == bands - 1)
        {
            //从i-j循环
            float temp1 = ((float)r[j] - (float)r[i]) / (w[j] - w[i]);
            float temp2 = (float)r[i] - temp1 * w[i];
            for (k = i; k <= j; k++)
            {
                tempHH = temp1 * w[k] + temp2;
                if (tempHH == 0)
                {
                    h[k] = 0.0f;
                }
                else
                {
                    h[k] = (float)r[k] / tempHH;
                }
            }
            break;
        }
        else
        {
            for (m = j + 1; m <bands; m++)
            {
                float temp1 = ((float)r[j] - (float)r[i]) / (w[j] - w[i]);
                float temp2 = (float)r[i] - temp1*w[i];
                float y = temp1 * w[m] + temp2;
                if (y >= r[m])
                {
                    if (m == bands - 1)
                    {
                        //从i-j循环
                        for (k = i; k <= j; k++)
                        {
                            tempHH = temp1*w[k] + temp2;
                            if (tempHH == 0)
                            {
                                h[k] = 0.0f;
                            }
                            else
                            {
                                h[k] = (float)r[k] / tempHH;
                            }
                        }
                        i = j;
                        break;
                    }
                }
                else
                {
                    break;
                }
            }
        }
        j = j + 1;
    }
}

函数：
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 返回输入序列x的上下包络线，作为其解析信号的大小。              %
% x的解析信号是利用希尔伯特实现的离散傅里叶变换得到的。           %
% 这个函数首先移除x的均值，然后在计算包络线之后再把它加回来。      %
% 如果x是一个矩阵，那么包络线对x的每一列都是独立的。             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[yupper,ylower] = envelope(x)




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 返回用解析信号的大小确定的x的包络线。                        %
% 解析信号通过使用长度为fl的Hilbert FIR滤波器对x进行滤波来计算。 %
% 如果只指定两个参数，则使用此语法。                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[yupper,ylower] = envelope(x,fl,'analytic')



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 返回x的上、下均方根包络。包络是使用长度wl样本的滑动窗口确定的。  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[yupper,ylower] = envelope(x,wl,'rms')



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 返回x的上峰和下峰包络线。                                  %
% 包络线由至少np个样本分离的局部最大值上的样条插值确定。         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[yupper,ylower] = envelope(x,np,'peak')



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 在没有输出参数的情况下，绘制信号及其上下包络线。               %
% 该语法接受以前语法中的任何输入参数。                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
envelope(___)

例1： [up,lo] = envelope(q)
% 信号生成
t = 0:1/2000:2-1/2000;
q = chirp(t-2,4,1/2,6,'quadratic',100,'convex').*exp(-4*(t-1).^2);
plot(t,q)


% 包络线绘制
[up,lo] = envelope(q);
hold on
plot(t,up,t,lo,'linewidth',1.5)
legend('q','up','lo')
hold off



例2： [up,lo] = envelope(q,100,'analytic');
% 信号生成
t = 0:1/1000:3;
q1 = sin(2*pi*7*t).*exp(-t/2);
q2 = chirp(t,30,2,5).*exp(-(2*t-3).^2)+2;
q = [q1;q2]';
plot(t,q)


% 包络线绘制
[up,lo] = envelope(q,100,'analytic');
hold on
plot(t,up,'-',t,lo,'--')
hold off


例3：不同参数的包络线区别
% 信号生成
t = 0.5:-1/100:-2.49;
z = airy(t*10).*exp(-t.^2);
plot(z)

% 无参数
figure
envelope(z)

% 参数：analytic,使用50-tap Hilbert filter计算信号的包络线
figure
envelope(z,50,'analytic')

% 参数：rms，使用40个样本的移动窗口来计算信号的RMS包络
figure
envelope(z,40,'rms')

% 参数：peak，确定峰包络，在被至少10个样本分隔的局部极大值上使用非结点条件的样条插值。
figure
envelope(z,10,'peak')










例4：不同采样点的影响





参考文献：

https://www.mathworks.com/help/signal/ref/envelope.html