一般來說，對均勻間隔采樣的數據的功率譜估計，可以采用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);]](http://thewe.net/tex/%28h_i,%20t_i%29 "(h_i, t_i)")
（所有數學公式采用latex語法書寫），
定義其數學期望和方差分別為
定義歸一化周期圖為
![[;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)} \} ;]](http://thewe.net/tex/P_N%28%5Comega%29%20=%20%5Cfrac%7B1%7D%7B2%20%5Cdelta%20%5E2%7D%20%20%5C%7B%20%5Cfrac%7B%5B%5Csum_i%20%28h_i%20-%20h%29%20%5Ccos%20%5Comega%28t_i%20-%20%5Ctau%29%5D%5E2%20%7B%5Csum_i%20%5Ccos%5E2%20%5Comega%20%28t_i%20-%20%5Ctau%29%7D%20%20+%20%5Cfrac%7B%5B%5Csum_i%20%28h_i%20-%20h%29%20%5Csin%20%5Comega%28t_i%20-%20%5Ctau%29%5D%5E2%20%7B%5Csum_i%20%5Csin%5E2%20%5Comega%20%28t_i%20-%20%5Ctau%29%7D%20%20%5C%7D "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;]](http://thewe.net/tex/%5Ctau "\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 閱讀(2597)
評論(3) 編輯 收藏 引用 所屬分類:
Numerical C++