频谱与功率谱的概念FFT与相关系数的C代码.docx
《频谱与功率谱的概念FFT与相关系数的C代码.docx》由会员分享，可在线阅读，更多相关《频谱与功率谱的概念FFT与相关系数的C代码.docx（15页珍藏版）》请在冰豆网上搜索。
频谱与功率谱的概念FFT与相关系数的C代码
频谱和功率谱有什么区别与联系
谱是个很不严格的东西,常常指信号的Fourier变换,是一个时间平均（timeaverage）概念功率谱的概念是针对功率有限信号的（能量有限信号可用能量谱分析）,所表现的是单位频带内信号功率随频率的变换情况。
保留频谱的幅度信息,但是丢掉了相位信息,所以频谱不同的信号其功率谱是可能相同的。
有两个重要区别:
1.功率谱是随机过程的统计平均概念,平稳随机过程的功率谱是一个确定函数;而频谱是随机过程样本的Fourier变换,对于一个随机过程而言,频谱也是一个“随机过程”。
（随机的频域序列）
2.功率概念和幅度概念的差别。
此外,只能对宽平稳的各态历经的二阶矩过程谈功率谱,其存在性取决于二阶局是否存在并且二阶矩的Fourier变换收敛;而频谱的存在性仅仅取决于该随机过程的该样本的Fourier变换是否收敛。
频谱分析（也称频率分析）,是对动态信号在频率域内进行分析,分析的结果是以频率为坐标的各种物理量的谱线和曲线,可得到各种幅值以频率为变量的频谱函数F（ω）。
频谱分析中可求得幅值谱、相位谱、功率谱和各种谱密度等等。
频谱分析过程较为复杂,它是以傅里叶级数和傅里叶积分为基础的。
功率谱
功率谱是个什么概念?
它有单位吗?
随机信号是时域无限信号,不具备可积分条件,因此不能直接进行傅氏变换。
一般用具有统计特性的功率谱来作为谱分析的依据。
功率谱与自相关函数是一个傅氏变换对。
功率谱具有单位频率的平均功率量纲。
所以标准叫法是功率谱密度。
通过功率谱密度函数,可以看出随机信号的能量随着频率的分布情况。
像白噪声就是平行于w轴,在w轴上方的一条直线。
功率谱密度,从名字分解来看就是说,观察对象是功率,观察域是谱域,通常指频域,密度,就是指观察对象在观察域上的分布情况。
一般我们讲的功率谱密度都是针对平稳随机过程的,由于平稳随机过程的样本函数一般不是绝对可积的,因此不能直接对它进行傅立叶分析。
可以有三种办法来重新定义谱密度,来克服上述困难。
一是用相关函数的傅立叶变换来定义谱密度;二是用随机过程的有限时间傅立叶变换来定义谱密度;三是用平稳随机过程的谱分解来定义谱密度。
三种定义方式对应于不同的用处,首先第一种方式前提是平稳随机过程不包含周期分量并且均值为零,这样才能保证相关函数在时差趋向于无穷时衰减,所以lonelystar说的不全对,光靠相关函数解决不了许多问题,要求太严格了;对于第二种方式,虽然一个平稳随机过程在无限时间上不能进行傅立叶变换,但是对于有限区间,傅立叶变换总是存在的,可以先架构有限时间区间上的变换,在对时间区间取极限,这个定义方式就是当前快速傅立叶变换（FFT）估计谱密度的依据;第三种方式是根据维纳的广义谐和分析理论:
Generalizedharmonicanalysis,ActaMath,55（1930）,117-258,利用傅立叶-斯蒂吉斯积分,对均方连续的零均值平稳随机过程进行重构,在依靠正交性来建立的。
另外,对于非平稳随机过程,也有三种谱密度建立方法,由于字数限制,功率谱密度的单位是G的平方/频率。
就是就是函数幅值的均方根值与频率之比。
是对随机振动进行分析的重要参数。
功率谱密度的国际单位是什么?
如果是加速度功率谱密度,加速度的单位是m/s^2,
那么,加速度功率谱密度的单位就是（m/s^2）^2/Hz,
而Hz的单位是1/s,经过换算得到加速度功率谱密度的单位是m^2/s^3.
同理,如果是位移功率谱密度,它的单位就是m^2*s,
如果是弯矩功率谱密度,单位就是（N*m）^2*s
位移功率谱——m^2*s
速度功率谱——m^2/s
加速度功率谱——m^2/s^3
求FFT实现的C++算法#include
#include"math.h"
classcomplex//复数类声明
{
private:
doublereal;
doubleimage;
public:
complex（doubler=0.0,doublei=0.0）//构造函数
{
real=r;
image=i;
}
complexoperator+（complexc2）;//+重载为成员函数
complexoperator-（complexc2）;//-重载为成员函数
complexoperator*（constcomplex&other）;//重载乘法
//complexoperator/（constcomplex&other）;//重载除法
voiddisplay（）;
};
complexcomplex:
:
operator+（complexc2）//重载的实现
{
complexc;
c.real=c2.real+real;
c.image=c2.image+image;
returncomplex（c.real,c.image）;
}
complexcomplex:
:
operator-（complexc2）//重载的实现
{
complexc;
c.real=real-c2.real;
c.image=image-c2.image;
returncomplex（c.real,c.image）;
}
complexcomplex:
:
operator*（constcomplex&other）
{
complextemp;
temp.real=real*other.real-image*other.image;
temp.image=image*other.real+real*other.image;
returntemp;
}
voidcomplex:
:
display（）
{
cout<<"（"<}
intfft（complex*a,intl）
//此处的l是级数数
{
constdoublepai=3.;
complexu,w,t;
intn=1,nv2,nm1,k,le,lei,ip;
inti,j,m;
doubletmp;
//n<n<<=l;
nv2=n>>1;
nm1=n-1;
i=0;
j=0;
for（i=0;i{
if（i{
t=a[j];
a[j]=a[i];
a[i]=t;
}
k=nv2;//求下一个倒位序数
while（k<=j）
{//原理是从高位加1,向低位进位
j-=k;
k>>=1;
}
j+=k;
}
//基2的输入倒位序,输出自然顺序的时间抽取FFT运算
le=1;
for（m=1;m<=l;m++）//当前运算的是第m级
{
lei=le;
le<<=1;
u=complex（1,0）;
tmp=pai/lei;
w=complex（cos（tmp）,-sin（tmp））;
//第m级中的不同系数的蝶形运算
for（j=0;j{
//第m级中的相同系数的蝶形
intaa=1;
for（i=j;i{
ip=i+lei;
t=a[ip]*u;
a[ip]=a[i]-t;
a[i]=a[i]+t;
}
//每级的旋转因子,由j和w来控制
u=u*w;
}
}
for（i=0;i{
cout<a[i].display（）;
}
return0;
}
实例:
voidmain（）
{
complexA[64];
for（inti=1;i<=64;i++）
{
A[i-1]=complex（cos（i*0.1）,cos（i*0.1+1））;
}
fft（A,6）;
for（intii=0;ii<64;ii++）
{
cout<A[ii].display（）;
}
}
本例与matlab中的FFT对比过,运算结果无误。
=======================================================================
1.#include
2.#defineDOUBLE_PI6.925286766559
3.
4.//快速傅里叶变换
5.//data长度为（2*2^n）,data的偶位为实数部分,data的奇位为虚数部分
6.//isInverse表示是否为逆变换
7.voidFFT（double*data,intn,boolisInverse=false）
8.{
9.intmmax,m,j,step,i;
10.doubletemp;
11.doubletheta,sin_htheta,sin_theta,pwr,wr,wi,tempr,tempi;
12.n=2*（113.intnn=n>>1;
14.//长度为1的傅里叶变换,位置交换过程
15.j=1;
16.for（i=1;in;i+=2）
17.{
18.if（j>i）
19.{
20.temp=data[j-1];
21.data[j-1]=data[i-1];
22.data[i-1]=temp;
23.data[j]=temp;
24.data[j]=data[i];
25.data[i]=temp;
26.}
27.//相反的二进制加法
28.m=nn;
29.while（m>=2&&j>m）
30.{
31.j-=m;
32.m>>=1;
33.}
34.j+=m;
35.}
36.//Danielson-Lanczos引理应用
37.mmax=2;
38.while（n>mmax）
39.{
40.step=mmax<1;
41.theta=DOUBLE_PI/mmax;
42.if（isInverse）
43.{
44.theta=-theta;
45.}
46.sin_htheta=sin（0.5*theta）;
47.sin_theta=sin（theta）;
48.pwr=-2.0*sin_htheta*sin_htheta;
49.wr=1.0;
50.wi=0.0;
51.for（m=1;mmmax;m+=2）
52.{
53.for（i=m;i=n;i+=step）
54.{
55.j=i+mmax;
56.tempr=wr*data[j-1]-wi*data[j];
57.tempi=wr*data[j]+wi*data[j-1];
58.data[j-1]=data[i-1]-tempr;
59.data[j]=data[i]-tempi;
60.data[i-1]+=tempr;
61.data[i]+=tempi;
62.}
63.sin_htheta=wr;
64.wr=sin_htheta*pwr-wi*sin_theta+wr;
65.wi=wi*pwr+sin_htheta*sin_theta+wi;
66.}
67.mmax=step;
68.}
69.}
70.
71./*输入数据为data,data是一组复数,偶数位存储的是复数的实数部分,奇数位存储的是复数的虚数部分。
data的长度与n相匹配。
注意:
这里的n并非是data的长度,data的实际长度为（2*2^n）,存储了N=2^n个复数。
72.
73.输出也存放在data中。
74.
75.以正向傅里叶变换为例,作为输入data中存储的是以delta为时间间隔时域函数的振幅抽样值。
经过函数
计算后data中存放输出,存储的是以1/（N*delta）为频率间隔频域像函数值。
频率范围为
0Hz,1/（N*delta）,2/（N*delta）...（N/2-1）/N*delta,+/-1/delta,-（N/2-1）/N*delta...-2/（N*delta）,-1/（N*delta）。
注意这是一个中间大两边小的排列。
76.
77.如果将isInverse设置为true则计算逆傅里叶变换。
*/
为了看明白那堆积分变换,不得不把复变函数扫了一遍,可看完了,才发现原来这堆变换说白了只是一些数字游戏,Examda提示:
也没用到啥复变函数的知识。
最后,用C++程序实现了下FFT,也算告一段落,代码如下:
#include
#include
#include
usingnamespacestd;
constdoublePI=3.46;
intn;//数据个数=2的logn次方
intlogn;
///复数结构体
structstCompNum
{
doublere;
doubleim;
};
stCompNum*pData1=NULL;
stCompNum*pData2=NULL;
///Examda提示:
正整数位逆序后输出
intreverseBits（intvalue,intbitCnt）
{
inti;
intret=0;
for（i=0;i{
ret|=（value&0x1）<<（bitCnt-1-i）;
value>>=1;
}
returnret;
}
voidmain（）
{
ifstreamfin（"data.txt"）;
inti,j,k;
//inputlogn
fin>>logn;
//calculaten
for（i=0,n=1;i//mallocmemoryspace
pData1=newstCompNum[n];
pData2=newstCompNum[n];
//inputrawdata
for（i=0;i>pData1[i].re;
for（i=0;i>pData1[i].im;
//FFTtransform
intcnt=1;
for（k=0;k{
for（j=0;j{
intlen=n/cnt;
doublec=-2*PI/len;
for（i=0;i{
intidx=len*j+i;
pData2[idx].re=pData1[idx].re+pData1[idx+len/2].re;pData2[idx].im=pData1[idx].im+pData1[idx+len/2].im;}
for（i=len/2;i{
doublewcos=cos（c*（i-len/2））;
doublewsin=sin（c*（i-len/2））;
intidx=len*j+i;
stCompNumtmp;
tmp.re=pData1[idx-len/2].re-pData1[idx].re;
tmp.im=pData1[idx-len/2].im-pData1[idx].im;
pData2[idx].re=tmp.re*wcos-tmp.im*wsin;
pData2[idx].im=tmp.re*wsin+tmp.im*wcos;
}
}
cnt<<=1;
stCompNum*pTmp=NULL;
pTmp=pData1;
pData1=pData2;
pData2=pTmp;
}
//resort
for（i=0;i{
intrev=reverseBits（i,logn）;
stCompNumtmp;
if（rev>i）
{
tmp=pData1[i];
pData1[i]=pData1[rev];
pData1[rev]=tmp;
}
}
//outputresultdata
for（i=0;icout<for（i=0;icout<//freememoryspace
delete[]pData1;
delete[]pData2;
fin.close（）;
system（"pause"）;
}
输入文件data.txt的内容如下:
4
2.24.56.78.510.212.314.516.219.321.225.229.436.439.245.250.1
0000000000000000
相关系数C++程序
该程序按照信号处理书中公式编写;
doublecorrcoef（double*d1,double*d2,intdataL）
{
inti;
floatxy=0,x=0,y=0,xsum=0,ysum=0;
doublecorrc;
for（i=0;i{
xsum+=d1[i];
ysum+=d2[i];
}
xsum/=dataL;
ysum/=dataL;
for（i=0;i{
d1[i]-=xsum;
d2[i]-=ysum;
x+=d1[i]*d1[i];
y+=d2[i]*d2[i];
xy+=d1[i]*d2[i];
}
corrc=abs（xy）/（sqrt（x）*sqrt（y））;
returncorrc;
}
升级会员 