非均匀取样数据的功率谱估计方法 --c++实现

一般来说，对均匀间隔采样的数据的功率谱估计，可以采用DFT/AR/MA/ARMA/MUSIC等方法估计；
如果由于实验条件限制，或者偶尔的数据丢失，采样的数据间隔时间／空间不是均匀间隔的，此时，大致可以有两种处理方式：
1）用插值的方法将数据转换为等间隔的，然后按等间隔采样方法处理；这种方法最大的缺陷是，在低频成分处会出现假的凸起。
2）采用Lomb方法处理；
本文将给出Lomb方法公式推导的结果，和C++代码实现；至于具体推导过程，请参阅

Jerry D. Scargle, Studies in Astronomical Time Series Analysis. II. Statistical Aspect of Spectral Analysis of Unevenly Spaced Data, the Astrophysical Journal, vol. 263, pp. 835--853

这里不再多说。

Lomb方法

假设数据点集为$(h_i, t_i)$（所有数学公式采用latex语法书写），

定义其数学期望和方差分别为

$$

h=/frac{1}{N} /sum_{i} h_i //

/delta^2=/frac{1}{N-1} /sum_{i}(h_i - h)^{2}

$$

定义归一化周期图为

$$

P_N(/omega) = /frac{1}{2 /delta ^2}

/{

/frac{[/sum_i (h_i - h) /cos /omega(t_i - /tau)]^2}{/sum_i /cos^2 /omega (t_i - /tau)}

+

/frac{[/sum_i (h_i - h) /sin /omega(t_i - /tau)]^2}{/sum_i /sin^2 /omega (t_i - /tau)}

 /}

$$

其中$/tau$由下式定义

$$

/tan 2 /omega /tau = /frac{/sum_i /sin 2 /omega t_i}{/sum_i /cos 2 /omega t_i}

$$

Lomb方法的C++实现：

 

 1 //lomb.hpp
 2 
 3 #ifndef LOMB_HPP_INCLUDED__
 4 #define LOMB_HPP_INCLUDED__
 5 
 6 #include <vector>
 7 using std::vector;
 8 
 9 //Warning:
10 //  If you would rather a self-defined type than
11 //  a builtin type, these methods MUST be offered:
12 //
13 //  operator=( const Type& other )
14 //  operator=( const long double& val )
15 //  operator+( const Type& lhs, const Type& rhs)
16 //  operator-( const Type& lhs, const Type& rhs)
17 //  operator*( const Type& lhs, const Type& rhs)
18 //  operator/( const Type& lhs, const Type& rhs)
19 //  sin( const Type& val )
20 //  cos( const Type& val )
21 //
22 template
23 <
24     typename Type = long double
25 >
26 class Lomb
27 {
28     public:
29         vector<Type> spectrum() const;
30 
31     public:
32         Lomb(   const vector<Type>& _h,
33                 const vector<Type>& _t );
34         ~Lomb();
35 
36     private:
37         struct Lomb_Data;
38         Lomb_Data* lomb_data;
39 };
40 
41 #include "lomb.tcc"
42 
43 #endif // LOMB_HPP_INCLUDED__
44

 

 

 

//lomb.tcc

#include "iterator.hpp"

#include <cassert>
#include <cmath>

#include <numeric>
#include <iterator>
#include <algorithm>
using namespace std;

template
<
    typename Type
>
struct Lomb<Type>::Lomb_Data
{
    vector<Type> h;
    vector<Type> t;

    Type f() const;
    Type tau( const unsigned int ) const;
    Type e() const;
    Type d() const;

};

template
<
    typename Type
>
Type Lomb<Type>::Lomb_Data::f()const
{
    const Type max = *( max_element(
            t.begin(), t.end() ) );
    const Type min = *( min_element(
            t.begin(), t.end() ) );
    const unsigned int N = h.size();
    const Type ans = N / ( max - min ) ;

    return ans;
}

template
<
    typename Type
>
Type Lomb<Type>::Lomb_Data::tau( const unsigned int n )const
{

    const Type omega =
        2.0L * 3.14159265358979L * ( this -> f() ) * n ;
    const Type sin_ =
        sin2_accumulate( t.begin(),
                         t.end(),
                         omega
                );
    const Type cos_ =
        cos2_accumulate( t.begin(),
                         t.end(),
                         omega
                );
    const Type ans = atan( sin_ / cos_ ) / ( 2.0L * omega );

    return ans;
}

template
<
    typename Type
>
Type Lomb<Type>::Lomb_Data::e()const
{
    const Type ans =
        mean(
            h.begin(),
            h.end(),
            Type_Keep<Type>()
            );

    return ans;
}

template
<
    typename Type
>
Type Lomb<Type>::Lomb_Data::d()const
{
    const Type ans =
        variance(
            h.begin(),
            h.end(),
            Type_Keep<Type>()
            );

    return ans;
}

template
<
    typename Type
>
Lomb<Type>::Lomb( const vector<Type>& _h,
            const vector<Type>& _t
          )
{
    assert( _t.size() == _h.size() );
    lomb_data = new Lomb_Data;
    lomb_data -> h = _h;
    lomb_data -> t = _t;
}

template
<
    typename Type
>
Lomb<Type>::~Lomb()
{
    delete lomb_data;
}

template
<
    typename Type
>
vector<Type> Lomb<Type>::spectrum() const
{
    const unsigned int N = (*lomb_data).h.size();

    const Type f = lomb_data -> f();
    const Type e = lomb_data -> e();

    const Type d = lomb_data -> d();

    vector<Type> ans;

    for ( unsigned int i = 0; i < N; ++i )
    {
        const Type tau = lomb_data -> tau( i );
        const Type omega = 6.283185307179586476L *
                                i * f;

        const Type h_cos_square =
            h_cos_square_accumulate(
                                        (*lomb_data).h.begin(), (*lomb_data).h.end(),
                                        (*lomb_data).t.begin(), (*lomb_data).t.end(),
                                        e,
                                        tau,
                                        omega
                                    );

        const Type cos_square =
            cos_square_accumulate(
                                    (*lomb_data).t.begin(),
                                    (*lomb_data).t.end(),
                                    tau,
                                    omega
                                 );

        const Type h_sin_square =
            h_sin_square_accumulate(
                                        (*lomb_data).h.begin(), (*lomb_data).h.end(),
                                        (*lomb_data).t.begin(), (*lomb_data).t.end(),
                                        e,
                                        tau,
                                        omega
                                    );

        const Type sin_square =
            sin_square_accumulate(
                                    (*lomb_data).t.begin(),
                                    (*lomb_data).t.end(),
                                    tau,
                                    omega
                                 );

        ans.push_back( sqrt(
                                (
                                    h_cos_square/cos_square +
                                    h_sin_square/sin_square
                                ) / (2.0L*d)
                            ) / N
                     );
    }

    return ans;
}

 

 

 

  1 //iterator.hpp
  2 
  3 #ifndef ITERATOR_HPP_INCLUDED__
  4 #define ITERATOR_HPP_INCLUDED__
  5 
  6 //this struct is just employed
  7 //to keep the Type information
  8 template
  9 <
 10     typename Type
 11 >
 12 struct Type_Keep;
 13 
 14 
 15 //calculate
 16 //      /sum_{i} /sin 2 /omega t_i
 17 template
 18 <
 19     typename InputItor,
 20     typename Type
 21 >
 22 Type sin2_accumulate(
 23                         InputItor begin,
 24                         InputItor end,
 25                         Type omega
 26                     );
 27 
 28 //calculate
 29 //      /sum_{i} /cos 2 /omega t_i
 30 template
 31 <
 32     typename InputItor,
 33     typename Type
 34 >
 35 Type cos2_accumulate(
 36                         InputItor begin,
 37                         InputItor end,
 38                         Type omega
 39                     );
 40 
 41 //calculate
 42 //      /frac{1}{N} /sum_{i} h_i
 43 template
 44 <
 45     typename InputItor,
 46     typename Type
 47 >
 48 Type mean(
 49             InputItor begin,
 50             InputItor end,
 51             Type_Keep<Type>
 52          );
 53 
 54 //calculate
 55 //      /frac{1}{N-1} /sum_{i}(h_i - h)^{2}
 56 template
 57 <
 58     typename InputItor,
 59     typename Type
 60 >
 61 Type variance(
 62                 InputItor begin,
 63                 InputItor end,
 64                 Type_Keep<Type>
 65              );
 66 
 67 //calcuate
 68 //       [/sum_i (h_i - h) /sin /omega (t_i - /tau)]^{2}
 69 template
 70 <
 71     typename Type,
 72     typename InputItor
 73 >
 74 Type h_sin_square_accumulate(     InputItor h_start, InputItor h_end,
 75                                 InputItor t_start, InputItor t_end,
 76                                 Type h,
 77                                 Type tau,
 78                                 Type omega
 79                             );
 80 
 81 //calcuate
 82 //       [/sum_i (h_i - h) /cos /omega (t_i - /tau)]^{2}
 83 template
 84 <
 85     typename Type,
 86     typename InputItor
 87 >
 88 Type h_cos_square_accumulate(     InputItor h_start, InputItor h_end,
 89                                 InputItor t_start, InputItor t_end,
 90                                 Type h,
 91                                 Type tau,
 92                                 Type omega
 93                             );
 94 
 95 //calculate
 96 //      /sum_{i} /cos ^{2} /omega (t_i - /tau)
 97 template
 98 <
 99     typename Type,
100     typename InputItor
101 >
102 Type cos_square_accumulate(     InputItor t_start,
103                                 InputItor t_end,
104                                 Type tau,
105                                 Type omega
106                             );
107 
108 //calculate
109 //      /sum_{i} /cos ^{2} /omega (t_i - /tau)
110 template
111 <
112     typename Type,
113     typename InputItor
114 >
115 Type sin_square_accumulate(     InputItor t_start,
116                                 InputItor t_end,
117                                 Type tau,
118                                 Type omega
119                             );
120 
121 #include "iterator.tcc"
122 
123 #endif // ITERATOR_HPP_INCLUDED__
124 

 

 

 

//iterator.tcc

#include <cmath>

template
<
    typename Type
>
struct Type_Keep
{
    typedef Type T;
};

template
<
    typename Type
>
struct Sin
{
    Type operator()( const Type val )const
    {
        return sin( val );
    }
};

template
<
    typename Type
>
struct Cos
{
    Type operator()( const Type val )const
    {
        return cos( val );
    }
};

template
<
    typename InputItor,
    typename Type,
    typename Function
>
static Type f2_accumulate(
                            InputItor begin,
                            InputItor end,
                            Type omega,
                            Function f
                        )
{
    Type ans = Type();
    while( begin != end )
    {
        ans += f( 2.0L * omega * (*begin++) );
    }
    return ans;
}

template
<
    typename InputItor,
    typename Type
>
Type sin2_accumulate(
        InputItor begin,
        InputItor end,
        Type omega
        )
{
    return f2_accumulate(
                            begin,
                            end,
                            omega,
                            Sin<Type>()
                         );
}

template
<
    typename InputItor,
    typename Type
>
Type cos2_accumulate(
        InputItor begin,
        InputItor end,
        Type omega
        )
{
    return f2_accumulate(
                            begin,
                            end,
                            omega,
                            Sin<Type>()
                        );
}

template
<
    typename InputItor,
    typename Type
>
Type mean(
        InputItor begin,
        InputItor end,
        Type_Keep<Type>
        )
{
    Type sum = Type();
    unsigned int cnt = 0;
    while ( begin != end )
    {
        sum += *begin++;
        ++cnt;
    }

    return sum / cnt;
}

template
<
    typename InputItor,
    typename Type
>
static Type diff_square_accumulate(
                    InputItor begin,
                    InputItor end,
                    Type mean
                    )
{
    Type ans = Type();
    while( begin != end )
    {
        ans += ( mean - *begin ) * ( mean - *begin ) ;
        ++begin;
    }
    return ans;
}

template
<
    typename InputItor
>
static unsigned int difference(
                InputItor begin,
                InputItor end
                )
{
    unsigned int cnt = 0;
    while( begin++ != end )
        ++cnt;

    return cnt;
}

template
<
    typename InputItor,
    typename Type
>
Type variance(
            InputItor begin,
            InputItor end,
            Type_Keep<Type>
            )
{
    const Type average =
        mean( begin, end, Type_Keep<long double>() );
    const Type power =
        diff_square_accumulate( begin, end, average );
    const unsigned int size =
        difference( begin, end );
    const Type ans = power / (size-1);

    return ans;
}

template
<
    typename Type,
    typename InputItor,
    typename Function
>
static Type h_f_square_accumulate(     InputItor h_start, InputItor h_end,
                                    InputItor t_start, InputItor t_end,
                                    Type h,
                                    Type tau,
                                    Type omega,
                                    Function f
                            )
{
    Type ans = Type();
    while( ( h_start != h_end ) && ( t_start != t_end ) )
    {
        ans += ( *h_start++ - h ) * f( omega * ( *t_start++ - tau) );
    }
    return ans*ans;
}


template
<
    typename Type,
    typename InputItor
>
Type h_sin_square_accumulate(     InputItor h_start, InputItor h_end,
                                InputItor t_start, InputItor t_end,
                                Type h,
                                Type tau,
                                Type omega
                            )
{
    return h_f_square_accumulate(
                                    h_start, h_end,
                                    t_start, t_end,
                                    h,
                                    tau,
                                    omega,
                                    Sin<Type>()
                                );
}

template
<
    typename Type,
    typename InputItor
>
Type h_cos_square_accumulate(     InputItor h_start, InputItor h_end,
                                InputItor t_start, InputItor t_end,
                                Type h,
                                Type tau,
                                Type omega
                            )
{
    return h_f_square_accumulate(
                                    h_start, h_end,
                                    t_start, t_end,
                                    h,
                                    tau,
                                    omega,
                                    Cos<Type>()
                                );
}

template
<
    typename Type,
    typename InputItor,
    typename Function
>
static Type f_square_accumulate(     InputItor t_start,
                                    InputItor t_end,
                                    Type tau,
                                    Type omega,
                                    Function f
                            )
{
    Type ans = Type();
    while( t_start != t_end )
    {
        const Type tmp = f( omega * ( *t_start++ - tau ) );
        ans += tmp * tmp;
    }
    return ans;
}

template
<
    typename Type,
    typename InputItor
>
Type cos_square_accumulate(     InputItor t_start,
                                InputItor t_end,
                                Type tau,
                                Type omega
                            )
{
    return f_square_accumulate(
                                    t_start,
                                    t_end,
                                    tau,
                                    omega,
                                    Cos<Type>()
                                );
}

template
<
    typename Type,
    typename InputItor
>
Type sin_square_accumulate(     InputItor t_start,
                                InputItor t_end,
                                Type tau,
                                Type omega
                            )
{
    return f_square_accumulate(
                                    t_start,
                                    t_end,
                                    tau,
                                    omega,
                                    Sin<Type>()
                                );
}

 

 

实例

下面这个例子中，将从文件phi.dat中输入采样数据，从文件rem.dat中输入采样空间间隔，使用Lomb算法处理后，得到的空间频谱将输出到文件spectrum.dat中去

 

 1 #include "lomb.hpp"
 2 
 3 #include <iostream>
 4 #include <fstream>
 5 #include <vector>
 6 #include <iterator>
 7 #include <algorithm>
 8 
 9 using namespace std;
10 
11 int main()
12 {
13     ifstream phi_src( "phi.dat" );
14     ifstream rem_src( "rem.dat" );
15     ofstream spectrum_dst( "spectrum.dat" );
16 
17     vector<long double> h;
18     vector<long double> t;
19 
20     copy( istream_iterator<long double>(phi_src), istream_iterator<long double>(), back_inserter(t) );
21     copy( istream_iterator<long double>(rem_src), istream_iterator<long double>(), back_inserter(h) );
22 
23     Lomb<long double>* lomb = new Lomb<long double>( h, t );
24 
25     vector<long double> v_spectrum = lomb -> spectrum();
26     copy( v_spectrum.begin(), v_spectrum.end(), ostream_iterator<long double>(spectrum_dst, "/n") );
27 
28     delete lomb;
29     phi_src.close();
30     rem_src.close();
31     spectrum_dst.close();
32 
33     return 0;
34 }
35 

 

 

备注

Lomb算法的复杂度是O(n^2)，若采样数据比较长，可以采用fft方法简化复杂度到O(nlogn)，与大数乘法中的fft用法一致，此处不再多说。