频谱与功率谱的概念-FFT与相关系数的C++代码

1、频谱和功率谱有什么区分与联系Fourier 变换， 是一个时间平均time average概念 功率谱的概念是针对功率有限信号的(能量有限信号可用能量谱分析)，所表现的是单位频带内信号功率随频率的变换状况。保存频谱的幅度信息，但是丢掉了相位信息，所以频谱不同的信号其功率谱是可能一样的。有两个重要区分：功率谱是随机过程的统计平均概念，平稳随机过程的功率谱是一个确定函数；而频谱是随机过程样本的Fourier变换，对于一个随机过程而言，频谱也是一个“随机过程”。随机的频域序列功率概念和幅度概念的差异。此外，只能对宽平稳的各态历经的二阶矩过程谈功率谱，其存在性取决于二阶局 Fourier Fourie

2、r 变换是否收敛。频谱分析也称频率分析，是对动态信号在频率域内进展分析，分析的结果是以频率为坐标的各种物理量的谱线和曲线，可得到各种幅值以频率为变量的频谱函数 F()。频谱分析中可求得幅值谱、相位谱、功率谱和各种谱密度等等。频谱分析过程较为简单，它是以傅里叶级数和傅里叶积分为根底的。功率谱功率谱是个什么概念？它有单位吗?随机信号是时域无限信号，不具备可积分条件，因此不能直接进展傅氏变换。一般用具有统计特性的功率谱来作为谱分析的依据。功率谱与自相关函数是一个傅氏变换对。功率谱具有单位频率的平均功率量纲。所以标准叫法是功率谱密度w 轴，在 w 轴上方的一条直线。功率谱密度，从名字分解来看就是说，观

3、看对象是功率，观看域是谱域，通常指频域，密度，就是指观看对象在观看域上的分布状况。一般我们讲的功率谱密度都是针对平稳随机过程的，由于平稳随机过程的样本函数一般不是确定可积的，因此不能直接对它进展傅立叶分析。可以有三种方法来重定义谱密度，来抑制上述困难。 一是用相关函数的傅立叶变换来定义谱密度；二是用随机过程的有限时间傅立叶变换来定义谱密度；三是用平稳随机过程的谱分解来定义谱密度。三种定义方式对应于不同的用处，首先第一种方式前提是平稳随机过程不包含周期重量并且均值为零，这样才能保证相关函数在时差趋向于无lonelystar 说的不全对，光靠相关函数解决不了很多问题，要求太严格了；对于其次种方式，

4、虽然一个平稳随机过程在无限时间上不能进展傅立叶变换，但是对于有限区间，傅立叶变换总是存在的，可以先架构有限时间区间上的变换，在对时间区间取极限，这个定义方式就是当前快速傅立叶变换FFT估量谱密度的依据；第三种 方式是依据维纳的广义谐和分析理论：Generalized harmonic analysis, Acta Math, 55(1930), 117-258,利用傅立叶-斯蒂吉斯积分，对均方连续的零均值平稳随机过程进展重构，在依靠正交性来建立的。另外，对于非平稳随机过程，也有三种谱密度建立方法，由于字数限制，功率谱密度的单位是 G 的平方/频率。就是就是函数幅值的均方根值与频率之比。是对随机

5、振动进展分析的重要参数。功率谱密度的国际单位是什么?假设是加速度功率谱密度，加速度的单位是 m/s2,那么，加速度功率谱密度的单位就是(m/s2)2/Hz,Hz 1/s,m2/s3.同理，假设是位移功率谱密度，它的单位就是 m2*s,假设是弯矩功率谱密度，单位就是(N*m)2*s 位移功率谱m2*s速度功率谱m2/s加速度功率谱m2/s3#include#include “math.h“classcomplex/复数类声明private:double real; double image;public:complex(double r=0.0,double i=0.0)/构造函数real=r;

6、 image=i;complex operator+(complex c2);/+重载为成员函数complexoperator-(complexc2);/-重载为成员函数complexoperator*(constcomplex&other);/重载乘法/complexoperator/(constcomplex&other);/重载除法void display;complex complex:operator +(complex c2)/重载的实现complex c; c.real=c2.real+real; c.image=c2.image+image;return complex(c.r

7、eal,c.image);complexcomplex:operator -(complex c2)/重载的实现complex c; c.real=real-c2.real;c.image=image-c2.image; returncomplex(c.real,c.image);complex complex:operator*(const complex& other)complex temp;temp.real=real*other.real-image*other.image; temp.image=image*other.real+real*other.image;return te

8、mp;void complex:displaycout“(“real“,“image“)“endl;int fft(complex *a,int l)/l 是级数数const double pai=3.141592653589793;complex u,w,t;int n=1,nv2,nm1,k,le,lei,ip;int i,j,m;double tmp;/nl 表示n l 位即系列的长度n1;nm1=n-1; i=0;j=0;for(i=0;inm1;i+)/.if(ij)/”变址运算t=aj; aj=ai;ai=t;k=nv2;while(k=1;j+=k;/2 的输入倒位序，输出自然挨

9、次的时间抽取FFT运算le=1; /m级lei=le; le=1;u=complex (1,0);tmp=pai/lei;w=complex (cos(tmp),-sin(tmp);for(j=0;jlei;j+)int aa=1;for(i=j;in;i+=le)ip=i+lei; t=aip*u; aip=ai-t;ai=ai+t; u=u*w;for(i=0;inm1+1;i+)coutiendl;ai.display;return 0;实例：void maincomplex A64;for (int i=1;i=64;i+)Ai-1=complex(cos(i*0.1),cos(i*0

10、.1+1);fft(A,6);for(int ii=0;ii64;ii+)coutiiendl;Aii.display;.3.=#include#defineDOUBLE_PI/快速傅里叶变换/data长度为 (2*2n)data 的偶位为实数局部data的奇位为虚数局部表示是否为逆变换voidFFT(double *data, intn,boolisInverse= false)intmmax,m,j,step,i;doubletemp;doubletheta,sin_htheta,sin_theta,pwr,wr,wi,tempr,tempi;n=2*(11;14. 1的傅里叶变换位置交换

11、过程15.j=1;16.for(i=1;i n;i+=2)17.18.if(ji)19.temp=dataj-1;dataj-1=datai-1;datai-1=temp;dataj=temp;dataj=datai;datai=temp;26./ 相反的二进制加法m=nn;while(m=2&jm)30.31.j-=m;32.m=1;33.34.j+=m;35./DanielsonLanczos 引理应用mmax=2;while(nmmax)39.step=mmax1;theta=DOUBLE_PI/mmax;if(isInverse)43.44.theta=-theta;45.sin_ht

12、heta=sin(0.5*theta);sin_theta=sin(theta);pwr=-2.0*sin_htheta*sin_htheta;49.wr=1.0;50.wi=0.0;51.for(m=1;m mmax;m+=2)52.53.for(i=m;i=n;i+=step)54.j=i+mmax;tempr=wr*dataj-1-wi*dataj;tempi=wr*dataj+wi*dataj-1;dataj-1=datai-1-tempr;dataj=datai-tempi;datai-1+=tempr;datai+=tempi;62.sin_htheta=wr;wr=sin_hth

13、eta*pwr-wi*sin_theta+wr;wi=wi*pwr+sin_htheta*sin_theta+wi;66.67.mmax=step;68.69. 70.71. ndata的长度，data的实际长度为(2 * 2n),存储了 N = 2n 个复数。72.73.data中。74.75.delta为时间间隔时域函数的振幅抽样值。经过函数1/(N*delta)为频率间隔频域像函数值。频率范围为0Hz,1/(N*delta),2/(N*deltaN/21/N*delta+/1/deltaN/21/N*delta-2/(N*delta)-1/(N*delta)。留意这是一个中间大两边小的排

14、列。76.77.true则计算逆傅里叶变换。*/提示: 也没用到啥复变函数的学问。最终，用C+FFT，也算告一段落，代码如下：#include #include #include usingnamespacestd;const double PI = 3.14159265358979323846; 次方int logn;/ 复数构造体struct stCompNumdouble re; doubleim;stCompNum* pData1 = NULL;stCompNum*pData2=NULL;/ Examda提示: 正整数位逆序后输出int reverseBits(int value, i

15、nt bitCnt)int i;int ret = 0;for(i=0; ibitCnt; i+)ret |= (value & 0 x1) = 1;return ret;void mainifstream fin(“data.txt“);int i,j,k;/ input lognfinlogn;/ calculate nfor(i=0, n=1; ilogn; i+) n *= 2;/ malloc memory space pData1 = new stCompNumn;pData2=newstCompNumn;/ input raw datafor(i=0; ipData1i.re;

16、for(i=0;ipData1i.im;/ FFT transformint cnt = 1;for(k=0; klogn; k+)for(j=0; jcnt; j+)int len = n / cnt;double c = - 2 * PI / len;for(i=0; ilen/2; i+)int idx = len * j + i;pData2idx.re = pData1idx.re + pData1idx + len/2.re; pData2idx.im=pData1idx.im+pData1idx+len/2.im;for(i=len/2; ilen; i+)double wcos

17、 = cos(c * (i - len/2);double wsin = sin(c * (i - len/2); int idx = len * j + i;stCompNum tmp;tmp.re = pData1idx - len/2.re - pData1idx.re; tmp.im = pData1idx - len/2.im - pData1idx.im;pData2idx.re = tmp.re * wcos - tmp.im * wsin; pData2idx.im = tmp.re * wsin + tmp.im * wcos;cnt = 1;stCompNum* pTmp

18、= NULL;pTmp = pData1;pData1 = pData2;pData2 = pTmp;/ resortfor(i=0; i i)tmp = pData1i; pData1i = pData1rev;pData1rev = tmp;/ output result datafor(i=0; in; i+) coutpData1i.re“t“;coutendl;for(i=0; in; i+) coutpData1i.im“t“;coutendl;/ free memory spacedelete pData1; delete pData2; fin.close;system(“pause“);data.txt 的内容如下：42.2 4.5 6.7 8.5 10.2 12.3 14.5 16.2 19.3 21.2 25.2 29.4 36.4 39.2 45.2 50.10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0相关系数C+程序该程序依据信号处理书中公式编写；double corrcoef(double *d1,