当前位置:文档之家› 频谱与功率谱的概念-FFT与相关系数的C++代码

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

频谱和功率谱有什么区别与联系谱是个很不严格的东西,常常指信号的Fourier变换,是一个时间平均(time average)概念功率谱的概念是针对功率有限信号的(能量有限信号可用能量谱分析),所表现的是单位频带内信号功率随频率的变换情况。

保留频谱的幅度信息,但是丢掉了相位信息,所以频谱不同的信号其功率谱是可能相同的。

有两个重要区别:1.功率谱是随机过程的统计平均概念,平稳随机过程的功率谱是一个确定函数;而频谱是随机过程样本的Fourier 变换,对于一个随机过程而言,频谱也是一个“随机过程”。

(随机的频域序列)2.功率概念和幅度概念的差别。

此外,只能对宽平稳的各态历经的二阶矩过程谈功率谱,其存在性取决于二阶局是否存在并且二阶矩的Fourier变换收敛;而频谱的存在性仅仅取决于该随机过程的该样本的Fourier变换是否收敛。

频谱分析(也称频率分析),是对动态信号在频率域内进行分析,分析的结果是以频率为坐标的各种物理量的谱线和曲线,可得到各种幅值以频率为变量的频谱函数F(ω)。

频谱分析中可求得幅值谱、相位谱、功率谱和各种谱密度等等。

频谱分析过程较为复杂,它是以傅里叶级数和傅里叶积分为基础的。

功率谱功率谱是个什么概念?它有单位吗?随机信号是时域无限信号,不具备可积分条件,因此不能直接进行傅氏变换。

一般用具有统计特性的功率谱来作为谱分析的依据。

功率谱与自相关函数是一个傅氏变换对。

功率谱具有单位频率的平均功率量纲。

所以标准叫法是功率谱密度。

通过功率谱密度函数,可以看出随机信号的能量随着频率的分布情况。

像白噪声就是平行于w轴,在w 轴上方的一条直线。

功率谱密度,从名字分解来看就是说,观察对象是功率,观察域是谱域,通常指频域,密度,就是指观察对象在观察域上的分布情况。

一般我们讲的功率谱密度都是针对平稳随机过程的,由于平稳随机过程的样本函数一般不是绝对可积的,因此不能直接对它进行傅立叶分析。

可以有三种办法来重新定义谱密度,来克服上述困难。

一是用相关函数的傅立叶变换来定义谱密度;二是用随机过程的有限时间傅立叶变换来定义谱密度;三是用平稳随机过程的谱分解来定义谱密度。

三种定义方式对应于不同的用处,首先第一种方式前提是平稳随机过程不包含周期分量并且均值为零,这样才能保证相关函数在时差趋向于无穷时衰减,所以lonelystar说的不全对,光靠相关函数解决不了许多问题,要求太严格了;对于第二种方式,虽然一个平稳随机过程在无限时间上不能进行傅立叶变换,但是对于有限区间,傅立叶变换总是存在的,可以先架构有限时间区间上的变换,在对时间区间取极限,这个定义方式就是当前快速傅立叶变换(FFT)估计谱密度的依据;第三种方式是根据维纳的广义谐和分析理论:Generalized harmonic analysis, Acta Math, 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<iostream.h>#include "math.h"class complex//复数类声明{private:double real;double image;public:complex(double r=0.0,double i=0.0)//构造函数{real=r;image=i;}complex operator+(complex c2);//+重载为成员函数complex operator-(complex c2);//-重载为成员函数complex operator*(const complex& other);//重载乘法//complex operator/(const complex& other);//重载除法void display();};complex complex::operator +(complex c2)//重载的实现{complex c;c.real=c2.real+real;c.image=c2.image+image;return complex(c.real,c.image);}complex complex::operator -(complex c2)//重载的实现{complex c;c.real=real-c2.real;c.image=image-c2.image;return complex(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 temp;}void complex::display(){cout<<"("<<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;//n<<l表示n左移l位即系列的长度n<<=l;nv2=n>>1;nm1=n-1;i=0;j=0;for(i=0;i<nm1;i++) //.....{if(i<j) //'''变址运算{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<lei;j++){//第m级中的相同系数的蝶形int aa=1;for(i=j;i<n;i+=le){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<nm1+1;i++){cout<<i<<endl;a[i].display();}return 0;}实例:void main(){complex A[64];for (int i=1;i<=64;i++){A[i-1]=complex(cos(i*0.1),cos(i*0.1+1));}fft(A,6);for(int ii=0;ii<64;ii++){cout<<ii<<endl;A[ii].display();}}本例与matlab中的FFT对比过,运算结果无误。

=======================================================================1.#include <math.h>2.#define DOUBLE_PI 6.2831853071795864769252867665593.4.// 快速傅里叶变换5.// data 长度为 (2 * 2^n), data 的偶位为实数部分, data 的奇位为虚数部分6.// isInverse表示是否为逆变换7.void FFT(double * data, int n, bool isInverse = false)8.{9.int mmax, m, j, step, i;10.double temp;11.double theta, sin_htheta, sin_theta, pwr, wr, wi, tempr, tempi;12. n = 2 * (1 < n);13.int nn = n >> 1;14.// 长度为1的傅里叶变换, 位置交换过程15. j = 1;16.for(i = 1; i n; 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; m mmax; 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是一组复数,偶数位存储的是复数的实数部分,奇数位存储的是复数的虚数部分。

相关主题