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就可以算出来了