非均匀取样数据的功率谱估计方法

一般来说，对均匀间隔采样的数据的功率谱估计，可以采用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=$

定义归一化周期图为

$[;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 =$

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用法一致，此处不再多说。

posted on 2009-01-02 21:20 Wang Feng 阅读(2503) 评论(3) 编辑收藏引用所属分类: Numerical C++

常用链接

留言簿(4)

随笔分类

随笔档案

Link List

搜索

最新评论

阅读排行榜

评论排行榜