posts - 12,  comments - 54,  trackbacks - 0
一般来说,对均匀间隔采样的数据的功率谱估计,可以采用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 阅读(1889) 评论(3)  编辑 收藏 引用 所属分类: Numerical C++

FeedBack:
# re: 非均匀取样数据的功率谱估计方法
2010-07-09 14:53 | JannieBallard
One knows that men's life seems to be expensive, however some people require money for different things and not every man earns big sums cash. Hence to get fast <a href="http://bestfinance-blog.com/topics/business-loans">business loans</a> and just consolidation loans would be a correct solution.   回复  更多评论
  
# re: 非均匀取样数据的功率谱估计方法
2010-07-18 01:00 | term paper help
There is nothing unnatural If you buy custom essays from the custom research paper service, however this thing would aid you to complete academic results.   回复  更多评论
  
# re: 非均匀取样数据的功率谱估计方法
2010-07-21 11:49 | dissertation writing service
Direct on what you are going to do at the time being and tomorrow, and have loyalty that the extent will advance when it is damn well on tap. Or maybe buy thesis is what you need . Thanks.   回复  更多评论
  

<2009年1月>
28293031123
45678910
11121314151617
18192021222324
25262728293031
1234567

常用链接

留言簿(4)

随笔分类

随笔档案

Link List

搜索

  •  

最新评论

阅读排行榜

评论排行榜