频谱与功率谱的概念FFT与相关系数的C代码.docx
- 文档编号:27553944
- 上传时间:2023-07-02
- 格式:DOCX
- 页数:15
- 大小:20.08KB
频谱与功率谱的概念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*(1 13.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 for(i=0;i //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;i cout< for(i=0;i cout< //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; }
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 频谱 功率 概念 FFT 相关系数 代码
![提示](https://static.bdocx.com/images/bang_tan.gif)