頻譜與功率譜的概念-FFT與相關(guān)系數(shù)的C++代碼

1、頻譜和功率譜有什么區(qū)分與聯(lián)系Fourier 變換， 是一個時間平均time average概念 功率譜的概念是針對功率有限信號的(能量有限信號可用能量譜分析)，所表現(xiàn)的是單位頻帶內(nèi)信號功率隨頻率的變換狀況。保存頻譜的幅度信息，但是丟掉了相位信息，所以頻譜不同的信號其功率譜是可能一樣的。有兩個重要區(qū)分：功率譜是隨機過程的統(tǒng)計平均概念，平穩(wěn)隨機過程的功率譜是一個確定函數(shù)；而頻譜是隨機過程樣本的Fourier變換，對于一個隨機過程而言，頻譜也是一個“隨機過程”。隨機的頻域序列功率概念和幅度概念的差異。此外，只能對寬平穩(wěn)的各態(tài)歷經(jīng)的二階矩過程談功率譜，其存在性取決于二階局 Fourier Fourie

2、r 變換是否收斂。頻譜分析也稱頻率分析，是對動態(tài)信號在頻率域內(nèi)進展分析，分析的結(jié)果是以頻率為坐標的各種物理量的譜線和曲線，可得到各種幅值以頻率為變量的頻譜函數(shù) F()。頻譜分析中可求得幅值譜、相位譜、功率譜和各種譜密度等等。頻譜分析過程較為簡單，它是以傅里葉級數(shù)和傅里葉積分為根底的。功率譜功率譜是個什么概念？它有單位嗎?隨機信號是時域無限信號，不具備可積分條件，因此不能直接進展傅氏變換。一般用具有統(tǒng)計特性的功率譜來作為譜分析的依據(jù)。功率譜與自相關(guān)函數(shù)是一個傅氏變換對。功率譜具有單位頻率的平均功率量綱。所以標準叫法是功率譜密度w 軸，在 w 軸上方的一條直線。功率譜密度，從名字分解來看就是說，觀

3、看對象是功率，觀看域是譜域，通常指頻域，密度，就是指觀看對象在觀看域上的分布狀況。一般我們講的功率譜密度都是針對平穩(wěn)隨機過程的，由于平穩(wěn)隨機過程的樣本函數(shù)一般不是確定可積的，因此不能直接對它進展傅立葉分析?？梢杂腥N方法來重定義譜密度，來抑制上述困難。 一是用相關(guān)函數(shù)的傅立葉變換來定義譜密度；二是用隨機過程的有限時間傅立葉變換來定義譜密度；三是用平穩(wěn)隨機過程的譜分解來定義譜密度。三種定義方式對應于不同的用處，首先第一種方式前提是平穩(wěn)隨機過程不包含周期重量并且均值為零，這樣才能保證相關(guān)函數(shù)在時差趨向于無lonelystar 說的不全對，光靠相關(guān)函數(shù)解決不了很多問題，要求太嚴格了；對于其次種方式，

4、雖然一個平穩(wěn)隨機過程在無限時間上不能進展傅立葉變換，但是對于有限區(qū)間，傅立葉變換總是存在的，可以先架構(gòu)有限時間區(qū)間上的變換，在對時間區(qū)間取極限，這個定義方式就是當前快速傅立葉變換FFT估量譜密度的依據(jù)；第三種 方式是依據(jù)維納的廣義諧和分析理論：Generalized harmonic analysis, Acta Math, 55(1930), 117-258,利用傅立葉-斯蒂吉斯積分，對均方連續(xù)的零均值平穩(wěn)隨機過程進展重構(gòu)，在依靠正交性來建立的。另外，對于非平穩(wěn)隨機過程，也有三種譜密度建立方法，由于字數(shù)限制，功率譜密度的單位是 G 的平方/頻率。就是就是函數(shù)幅值的均方根值與頻率之比。是對隨機

5、振動進展分析的重要參數(shù)。功率譜密度的國際單位是什么?假設(shè)是加速度功率譜密度，加速度的單位是 m/s2,那么，加速度功率譜密度的單位就是(m/s2)2/Hz,Hz 1/s,m2/s3.同理，假設(shè)是位移功率譜密度，它的單位就是 m2*s,假設(shè)是彎矩功率譜密度，單位就是(N*m)2*s 位移功率譜m2*s速度功率譜m2/s加速度功率譜m2/s3#include#include “math.h“classcomplex/復數(shù)類聲明private:double real; double image;public:complex(double r=0.0,double i=0.0)/構(gòu)造函數(shù)real=r;

6、 image=i;complex operator+(complex c2);/+重載為成員函數(shù)complexoperator-(complexc2);/-重載為成員函數(shù)complexoperator*(constcomplex&other);/重載乘法/complexoperator/(constcomplex&other);/重載除法void display;complex complex:operator +(complex c2)/重載的實現(xiàn)complex c; c.real=c2.real+real; c.image=c2.image+image;return complex(c.r

7、eal,c.image);complexcomplex:operator -(complex c2)/重載的實現(xiàn)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 是級數(shù)數(shù)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 的偶位為實數(shù)局部data的奇位為虛數(shù)局部表示是否為逆變換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 個復數(shù)。72.73.data中。74.75.delta為時間間隔時域函數(shù)的振幅抽樣值。經(jīng)過函數(shù)1/(N*delta)為頻率間隔頻域像函數(shù)值。頻率范圍為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則計算逆傅里葉變換。*/提示: 也沒用到啥復變函數(shù)的學問。最終，用C+FFT，也算告一段落，代碼如下：#include #include #include usingnamespacestd;const double PI = 3.14159265358979323846; 次方int logn;/ 復數(shù)構(gòu)造體struct stCompNumdouble re; doubleim;stCompNum* pData1 = NULL;stCompNum*pData2=NULL;/ Examda提示: 正整數(shù)位逆序后輸出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 的內(nèi)容如下：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相關(guān)系數(shù)C+程序該程序依據(jù)信號處理書中公式編寫；double corrcoef(double *d1,