detrend去趋势函数的Matlab、Python与C实现

趋势分量对频域分析的影响

在对信号做频域分析时，如果有趋势项的存在，会对分析形成干扰。仿真一个有斜率的线性信号，仿真信号及其经过傅里叶转换后的频谱如下图（左图为时域图、右图为频域图）。由此可见斜坡信号中低频信号幅值较大，随着频率增加，幅值减小。

我们再来仿真一个线性信号与正弦信号的叠加信号。正弦信号周期为100s，幅值为0.2。该信号的时域图和频谱图如下图所示（左图为时域图、右图为频域图）。由于我们已知正弦信号的周期，因此这个信号应该在频率为0.01Hz处，如红圈标注所示。但是我们会发现红圈处的幅值并不是0.2，对比上图可知，这其实是由于线性信号在0.01Hz处也有一个分量，所以两个信号叠加之后，幅值比0.2大。

由上述示例可知，趋势项的存在会影响频谱分析，因此需要对信号或者数据去趋势。

detrend去趋势函数（Matlab、Python）

Maltab对该函数的官方介绍detrend, Remove polynomial trend。detrend函数通过设置参数，不仅可以移除线性趋势，也能移除多阶趋势。

下面是一个官方示例：

t = 0:20;
x = 3*sin(t) + t;
y = detrend(x);
plot(t,x,t,y,t,x-y,':k')
legend('Input Data','Detrended Data','Trend','Location','northwest') 

Python中也有这个函数，官方介绍scipy.signal.detrend，不过只能移除线性趋势。如果要移除多阶趋势，需要使用obspy.signal.detrend.polynomial。

detrend的C语言实现

在博客 预处理丨去趋势(Matlab和C++)中提到了来自lppier 的C代码实现，不过这里的detrend函数目前只能移除线性趋势，代码如下。

/************************************************************************************
    Function    : void detrend_IP(T *y, T *x, int m)
    Description : Remove the linear trend of the input floating point data. Note that this
                  will initialize a work buffer inside the function. So if you are calling
                  this many, many times, create your work buffer in the calling scope and call
                  detrend(T *y, T*x, int m) instead to avoid initializing memory over and over
                  again.
    Inputs      : y - Floating point input data
                 m - Input data length
    Outputs     : y - Data with linear trend removed
    Copyright   : DSO National Laboratories
    History     : 01/02/2008, TCK, Adapted from HYC code
                  01/12/2008, TCK, Added in return value
                  25/01/2016, Pier, Changed into template type, removed need for work buffer
    *************************************************************************************/
template<typename T>
void Detrend::detrend_IP(T *y, int m)
{
    T xmean, ymean;
    int i;
    T temp;
    T Sxy;
    T Sxx;

    T grad;
    T yint;

    std::unique_ptr<T[]> x(new T[m]);

    /********************************
    Set the X axis Liner Values
    *********************************/
    for (i = 0; i < m; i++)
        x[i] = i;

    /********************************
    Calculate the mean of x and y
    *********************************/
    xmean = 0;
    ymean = 0;
    for (i = 0; i < m; i++)
    {
        xmean += x[i];
        ymean += y[i];
    }
    xmean /= m;
    ymean /= m;

    /********************************
    Calculate Covariance
    *********************************/
    temp = 0;
    for (i = 0; i < m; i++)
        temp += x[i] * y[i];
    Sxy = temp / m - xmean * ymean;

    temp = 0;
    for (i = 0; i < m; i++)
        temp += x[i] * x[i];
    Sxx = temp / m - xmean * xmean;

    /********************************
    Calculate Gradient and Y intercept
    *********************************/
    grad = Sxy / Sxx;
    yint = -grad * xmean + ymean;

    /********************************
    Removing Linear Trend
    *********************************/
    for (i = 0; i < m; i++)
        y[i] = y[i] - (grad * i + yint);
}

该段代码的实现原理实际使用了简单线性回归模型。
趋势项模型如下：

优化的指标函数设置为残差的平方，残差越小，则找了效果最好的线性回归：

为了使残差最小，通过对b求偏导，使其导数为0，则：

对a求偏导，使其导数为0，则：

将公式对照代码，就可以知道斜率、截距参数是如何计算得到的。