A. 功率譜是什麼
功率譜
周期運動在功率譜中對應尖鋒，混沌的特徵是譜中出現"雜訊背景"和寬鋒。它是研究系統從分岔走向混沌的重要方法。 在很多實際問題中（尤其是對非線性電路的研究）常常只給出觀測到的離散的時間序列X1, X2, X3,...Xn,那麼如何從這些時間序列中提取前述的四種吸引子（零維不動點、一維極限環、二維環面、奇怪吸引子）的不同狀態的信息呢？ 我們可以運用數學上已經嚴格證明的結論，即擬合。我們將N個采樣值加上周期條件Xn+i=Xi，則自關聯函數（即離散卷積）為 然後對Cj完成離散傅氏變換，計算傅氏系數。 Pk說明第k個頻率分量對Xi的貢獻，這就是功率譜的定義。當採用快速傅氏變換演算法後，可直接由Xi作快速傅氏變換，得到系數 然後計算 ，由許多組{Xi}得一批{Pk'}，求平均後即趨近前面定義的功率譜Pk。 從功率譜上，四種吸引子是容易區分的，如圖12 (a)，(b)對應的是周期函數，功率譜是分離的離散譜 (c)對應的是准周期函數，各頻率中間的間隔分布不像(b)那樣有規律。 (d)圖是混沌的功率譜，表現為"雜訊背景"及寬鋒。 考慮到實際計算中，數據只能取有限個，譜也總以有限分辨度表示出來，從物理實驗和數值計算的角度看，一個周期十分長的解和一個混沌解是難於區分的，這也正是功率譜研究的主要弊端。
B. 一個關於128點的快速傅立葉的c語言程序
這是我寫的1024點的快速傅里葉變換程序，下面有驗證，你把數組
double
A[2049]={0};
double
B[1100]={0};
double
powerA[1025]={0};
改成
A[256]={0};
B[130]={0};
power[129]={0};就行了，
void
FFT(double
data[],
int
nn,
int
isign)
的程序可以針對任何點數，只要是2的n次方
具體程序如下：
#include
<iostream.h>
#include
"math.h"
#include<stdio.h>
#include<string.h>
#include
<stdlib.h>
#include
<fstream.h>
#include
<afx.h>
void
FFT(double
data[],
int
nn,
int
isign)
{
//復數的快速傅里葉變換
int
n,j,i,m,mmax,istep;
double
tempr,tempi,theta,wpr,wpi,wr,wi,wtemp;
n
=
2
*
nn;
j
=
1;
for
(i
=
1;
i<=n
;
i=i+2)
//這個循環進行的是碼位倒置。
{
if(
j
>
i)
{
tempr
=
data[j];
tempi
=
data[j
+
1];
data[j]
=
data[i];
data[j
+
1]
=
data[i
+
1];
data[i]
=
tempr;
data[i
+
1]
=
tempi;
}
m
=
n
/
2;
while
(m
>=
2
&&
j
>
m)
{
j
=
j
-
m;
m
=
m
/
2;
}
j
=
j
+
m;
}
mmax
=
2;
while(
n
>
mmax
)
{
istep
=
2
*
mmax;
//這里表示一次的數字的變化。也體現了級數，若第一級時，也就是書是的第0級，其為兩個虛數，所以對應數組應該增加4，這樣就可以進入下一組運算
theta
=
-6.28318530717959
/
(isign
*
mmax);
wpr
=
-2.0
*
sin(0.5
*
theta)*sin(0.5
*
theta);
wpi
=
sin(theta);
wr
=
1.0;
wi
=
0.0;
for(
m
=
1;
m<=mmax;
m=m+2)
{
for
(i
=
m;
i<=n;
i=i+istep)
{
j
=
i
+
mmax;
tempr=double(wr)*data[j]-double(wi)*data[j+1];//這兩句表示蝶形因子的下一個數乘以W因子所得的實部和虛部。
tempi=double(wr)*data[j+1]+double(wi)*data[j];
data[j]
=
data[i]
-
tempr;
//蝶形單元計算後下面單元的實部，下面為虛部，注意其變換之後的數組序號與書上蝶形單元是一致的
data[j
+
1]
=
data[i
+
1]
-
tempi;
data[i]
=
data[i]
+
tempr;
data[i
+
1]
=
data[i
+
1]
+
tempi;
}
wtemp
=
wr;
wr
=
wr
*
wpr
-
wi
*
wpi
+
wr;
wi
=
wi
*
wpr
+
wtemp
*
wpi
+
wi;
}
mmax
=
istep;
}
}
void
main()
{
//本程序已經和MATLAB運算結果對比，准確無誤，需要注意的的是，計算中數組都是從1開始取得，丟棄了A[0]等數據
double
A[2049]={0};
double
B[1100]={0};
double
powerA[1025]={0};
char
line[50];
char
dataA[20],
dataB[20];
int
ij;
char
ch1[3]="\t";
char
ch2[3]="\n";
int
strl1,strl2;
CString
str1,str2;
ij=1;
//********************************讀入文件data1024.txt中的數據,
其中的數據格式見該文件
FILE
*fp
=
fopen("data1024.txt","r");
if(!fp)
{
cout<<"Open
file
is
failing!"<<endl;
return;
}
while(!feof(fp))
//feof(fp)有兩個返回值:如果遇到文件結束，函數feof（fp）的值為1，否則為0。
{
memset(line,0,50);
//清空為0
memset(dataA,0,20);
memset(dataB,0,20);
fgets(line,50,fp);
//函數的功能是從fp所指文件中讀入n-1個字元放入line為起始地址的空間內
sscanf(line,
"%s%s",
dataA,
dataB);
//我同時讀入了兩列值，但你要求1024個，那麼我就只用了第一列的1024個值
//dataA讀入第一列，dataB讀入第二列
B[ij]=atof(dataA);
//將字元型的dataA值轉化為float型
ij++;
}
for
(int
mm=1;mm<1025;mm++)//A[2*mm-1]是實部，A[2*mm]是虛部，當只要輸入實數時，那麼保證虛部A[mm*2]為零即可
{
A[2*mm-1]=B[mm];
A[2*mm]=0;
}
//*******************************************正式計算FFT
FFT(A,1024,1);
//********************************************寫入數據到workout.txt文件中
for
(int
k=1;k<2049;k=k+2)
{
powerA[(k+1)/2]=sqrt(pow(A[k],2.0)+pow(A[k+1],2.0));//求功率譜
FILE
*pFile=fopen("workout.txt","a+");
//?a+只能在文件最後補充，游標在結尾。沒有則創建
memset(ch1,0,15);
str1.Format("%.4f",powerA[(k+1)/2]);
if
(A[k+1]>=0)
str2.Format("%d\t%6.4f%s%6.4f
%s",(k+1)/2,A[k],"+",A[k+1],"i");//保存fft計算的頻譜，是復數頻譜
else
str2.Format("%d\t%6.4f%6.4f
%s",(k+1)/2,A[k],A[k+1],"i");
strl1=strlen(str1);
strl2=strlen(str2);
//
用
法:fwrite(buffer,size,count,fp);
//
buffer：是一個指針，對fwrite來說，是要輸出數據的地址。
//
size：要寫入的位元組數；
//
count:要進行寫入size位元組的數據項的個數；
//
fp:目標文件指針。
fwrite(str2,1,strl2,pFile);
fwrite(ch1,1,3,pFile);
fwrite(ch1,1,3,pFile);
fwrite(str1,1,strl1,pFile);
fwrite(ch2,1,3,pFile);
fclose(pFile);
}
cout<<"計算完畢，到fft_test\workout.txt查看結果"<<endl;
}
C. 用matlab實現功率譜
n=0:0.1:200;%設定信號時間長度為0到200秒，采樣間隔0.1，則采樣頻率為10HZ，點數2001
y=sin(2*pi*0.2*n)+sin(2*0.213*n);
Y=fft(y);%FFTPyy=Y.*conj(Y)/2000;%信號功率譜f=10*(0:1000)/2000;%計算橫軸頻率值figure(1)subplot(2,1,1),plot(n,y),title('信號'),xlabel('時間(S)')subplot(2,1,2),plot(f,Pyy(1:1001)),title('信號功率譜'),xlabel('頻率(Hz)')
請採納答案，支持我一下。
D. 請問有人知道 功率譜估計C++的源程序嗎到那裡可以下載
//lomb.hpp
#ifndef LOMB_HPP_INCLUDED__
#define LOMB_HPP_INCLUDED__
#include <vector>
using std::vector;
//Warning:
// If you would rather a self-defined type than
// a builtin type, these methods MUST be offered:
//
// operator=( const Type& other )
// operator=( const long double& val )
// operator+( const Type& lhs, const Type& rhs)
// operator-( const Type& lhs, const Type& rhs)
// operator*( const Type& lhs, const Type& rhs)
// operator/( const Type& lhs, const Type& rhs)
// sin( const Type& val )
// cos( const Type& val )
//
template
<
typename Type = long double
>
class Lomb
{
public:
vector<Type> spectrum() const;
public:
Lomb( const vector<Type>& _h,
const vector<Type>& _t );
~Lomb();
private:
struct Lomb_Data;
Lomb_Data* lomb_data;
};
#include "lomb.tcc"
#endif // LOMB_HPP_INCLUDED__
//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;
}
//iterator.hpp
#ifndef ITERATOR_HPP_INCLUDED__
#define ITERATOR_HPP_INCLUDED__
//this struct is just employed
//to keep the Type information
template
<
typename Type
>
struct Type_Keep;
//calculate
// \sum_{i} \sin 2 \omega t_i
template
<
typename InputItor,
typename Type
>
Type sin2_accumulate(
InputItor begin,
InputItor end,
Type omega
);
//calculate
// \sum_{i} \cos 2 \omega t_i
template
<
typename InputItor,
typename Type
>
Type cos2_accumulate(
InputItor begin,
InputItor end,
Type omega
);
//calculate
// \frac{1}{N} \sum_{i} h_i
template
<
typename InputItor,
typename Type
>
Type mean(
InputItor begin,
InputItor end,
Type_Keep<Type>
);
//calculate
// \frac{1}{N-1} \sum_{i}(h_i - h)^{2}
template
<
typename InputItor,
typename Type
>
Type variance(
InputItor begin,
InputItor end,
Type_Keep<Type>
);
//calcuate
// [\sum_i (h_i - h) \sin \omega (t_i - \tau)]^{2}
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
);
//calcuate
// [\sum_i (h_i - h) \cos \omega (t_i - \tau)]^{2}
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
);
//calculate
// \sum_{i} \cos ^{2} \omega (t_i - \tau)
template
<
typename Type,
typename InputItor
>
Type cos_square_accumulate( InputItor t_start,
InputItor t_end,
Type tau,
Type omega
);
//calculate
// \sum_{i} \cos ^{2} \omega (t_i - \tau)
template
<
typename Type,
typename InputItor
>
Type sin_square_accumulate( InputItor t_start,
InputItor t_end,
Type tau,
Type omega
);
#include "iterator.tcc"
#endif // ITERATOR_HPP_INCLUDED__
//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>()
);
}
#include "lomb.hpp"
#include <iostream>
#include <fstream>
#include <vector>
#include <iterator>
#include <algorithm>
using namespace std;
int main()
{
ifstream phi_src( "phi.dat" );
ifstream rem_src( "rem.dat" );
ofstream spectrum_dst( "spectrum.dat" );
vector<long double> h;
vector<long double> t;
( istream_iterator<long double>(phi_src), istream_iterator<long double>(), back_inserter(t) );
( istream_iterator<long double>(rem_src), istream_iterator<long double>(), back_inserter(h) );
Lomb<long double>* lomb = new Lomb<long double>( h, t );
vector<long double> v_spectrum = lomb -> spectrum();
( v_spectrum.begin(), v_spectrum.end(), ostream_iterator<long double>(spectrum_dst, "\n") );
delete lomb;
phi_src.close();
rem_src.close();
spectrum_dst.close();
return 0;
}
Lomb演算法的復雜度是O(n^2)，若采樣數據比較長，可以採用fft方法簡化復雜度到O(nlogn)，與大數乘法中的fft用法一致，此處不再多說。
從網上來的資料！
E. 請教懂信號，通信的大佬：頻譜-功率譜-功率譜密度到底是什麼
功率譜和頻譜的區別
1、計算
功率譜的計算需要信號先做自相關，然後再進行FFT運算。
頻譜的計算則是將信號直接進行FFT就行了。
2、方式
功率譜是對信號研究，不過它是從能量的方面來對信號研究的。
而頻譜也是用來形容信號的,只是的表示方式變了，從時域轉變成了頻域表示,也就是說一種信號的表示方式不同而已。
F. 計算功率譜密度的兩種方法,直接計算和periodogram法求解
是通過積分求得的！求解方法：1、直接法（又稱周期圖法），它是把隨機序列x(n)的N個觀測數據視為一能量有限的序列，直接計算x(n)的離散傅立葉變換，得X(k)，然後再取其幅值的平方，並除以N，作為序列x(n)真實功率譜的估計。2、間接法先由序列
G. 信號功率譜怎麼計算
用FFT求取信號頻譜的實部和虛部，實部的平方價虛部的平方就是功率譜。周期性連續信號x(t)的頻譜可表示為離散的非周期序列Xn,它的幅度頻譜的功率譜平方│Xn│2所排成的序列，就被稱之為該周期信號的「功率譜」。
對信號進行傅里葉變換，取sin部分為實部，cos部分為虛部，直接算實部和虛部的平方和，得到的就是頻域功率譜的分布 推薦使用matlab計算，因為一個函數FFT就可以算出來。
信號x(t)的功率譜密度計算方法：
1、先計算x(t)的傅立葉變換:X(jw),
2、取模:|X(jw)|,再平方:|X(jw)|^2,再除以樣本長度: |X(jw)|^2/T
3、就得到: x(t)的功率譜密度函數: Gxx(w)= |X(jw)|^2/T
(7)c語言功率譜擴展閱讀：
周期運動在功率譜中對應尖鋒，混沌的特徵是譜中出現"雜訊背景"和寬鋒。它是研究系統從分岔走向混沌的重要方法。 在很多實際問題中（尤其是對非線性電路的研究）常常只給出觀測到的離散的時間序列X1, X2, X3,Xn,那麼如何從這些時間序列中提取前述的四種吸引子（零維不動點、一維極限環、二維環面、奇怪吸引子）的不同狀態的信息。
可以運用數學上已經嚴格證明的結論，即擬合。我們將N個采樣值加上周期條件Xn+i=Xi，則自關聯函數（即離散卷積）為 然後對Cj完成離散傅氏變換，計算傅氏系數。 Pk說明第k個頻率分量對Xi的貢獻，這就是功率譜的定義。
H. 功率譜的定義
功率譜
周期運動在功率譜中對應尖鋒，混沌的特徵是譜中出現"雜訊背景"和寬鋒。它是研究系統從分岔走向混沌的重要方法。
在很多實際問題中（尤其是對非線性電路的研究）常常只給出觀測到的離散的時間序列X1,
X2,
X3,...Xn,那麼如何從這些時間序列中提取前述的四種吸引子（零維不動點、一維極限環、二維環面、奇怪吸引子）的不同狀態的信息呢？
我們可以運用數學上已經嚴格證明的結論，即擬合。我們將N個采樣值加上周期條件Xn+i=Xi，則自關聯函數（即離散卷積）為
然後對Cj完成離散傅氏變換，計算傅氏系數。
Pk說明第k個頻率分量對Xi的貢獻，這就是功率譜的定義。當採用快速傅氏變換演算法後，可直接由Xi作快速傅氏變換，得到系數
然後計算
，由許多組{Xi}得一批{Pk'}，求平均後即趨近前面定義的功率譜Pk。
從功率譜上，四種吸引子是容易區分的，如圖12
(a)，(b)對應的是周期函數，功率譜是分離的離散譜
(c)對應的是准周期函數，各頻率中間的間隔分布不像(b)那樣有規律。
(d)圖是混沌的功率譜，表現為"雜訊背景"及寬鋒。
考慮到實際計算中，數據只能取有限個，譜也總以有限分辨度表示出來，從物理實驗和數值計算的角度看，一個周期十分長的解和一個混沌解是難於區分的，這也正是功率譜研究的主要弊端。
I. FFT和功率譜的結果是什麼含義
FFT：將時域波形轉化為頻譜
功率譜：信號功率隨頻率的變化。
http://wenku..com/view/2c57730c6c85ec3a87c2c56d.html?from=rec&pos=0&weight=11&lastweight=4&count=5
J. 功率譜計算
對信號進行傅里葉變換，取sin部分為實部，cos部分為虛部，直接算實部和虛部的平方和，得到的就是頻域功率譜的分布 推薦使用matlab計算，因為一個函數FFT就可以算出來了