AutoCAD 3DMAX C语言 Pro/E UG JAVA编程 PHP编程 Maya动画 Matlab应用 Android
Photoshop Word Excel flash VB编程 VC编程 Coreldraw SolidWorks A Designer Unity3D
 首页 > VC编程

VC++编程实现对波形数据的频谱分析

51自学网 2015-08-30 http://www.wanshiok.com
  摘要: 本文介绍了采用离散傅立叶变换(DFT)实现对采样得到的波形数据文件进行频谱分析的一般方法,并且为了提高运算效率、节省中间存储单元,最终采用了"时间抽选奇偶分解"的"库利-图基算法"实现快速离散傅立叶变换,对采样数据进行了高效的频谱分析,并用Microsoft Visual C++ 6.0编写实现。

  关键字:Microsoft Visual C ++ 6.0、离散傅立叶变换、快速傅立叶变换、采样

  一、 引言

  频谱分析是电子工程上一个非常重要的手段,许多计算机辅助电路分析(CAA)类软件都具备这种分析能力,以便电子工程师能清楚的看到某波形的频谱分布情况。而要对一个输入信号源作频谱分析,将其由时域信号转变为频域信号,就必然要用到傅立叶分析,而无论是在时域还是在频域,都要对连续函数进行积分运算。很显然,要通过计算机实现此变换必须预先通过抽样将原始的连续数据转变为离散数据,并将计算范围收缩到一个有限区间。因此在允许一定程度近似的条件下,可以使用"离散傅立叶变换(DFT)"对波形数据进行频谱分析。

  二、 快速傅立叶变换(FFT)算法构成原理

  要计算一个N点的离散傅立叶变换需要同一个N*N点的W矩阵(关于W矩阵请参阅信号与系统方面的书籍)相运算,随着N值的增大,运算次数显著上升,当点数达到1024时,需要进行复数乘法运算1,048,576次,显然这种算法在实际运用中无法保证当点数较大时的运算速度,无法满足对信号的实时处理。

  根据W矩阵中W元素的周期性和对称性我们可以将一个N点的DFT运算分解为两组N/2点的DFT运算,然后取和即可,为进一步提高效率,将上述两个矩阵按奇偶顺序逐级分解下去。当采样点数为2的指数次方M时,可分解为M级子矩阵运算,全部工作量仅为:

  复数乘法:M*N/2次

  复数加法:N*M次

  而直接DFT需要的运算量为:

  复数乘法:N*N次

  复数加法:N*(N-1)次

  当点数N为几十个点时FFT的优势还不明显,而一旦达到几千、几百个点时优势是十分明显的:

  N=1024时:DFT需1048576次运算,FFT仅需5120次运算,改善比204.8。

  N=2048时:DFT需4194304次运算,FFT仅需11264次运算,改善比达到372.4。

  三、 "时间抽选奇偶分解快速离散傅立叶变换"的程序实现

  当采样点数较多时,如变换前和变换后的序列都按自然顺序排列,则中间运算过程会占用大量的中间存储单元,造成效率的低下和存储单元的浪费。根据FFT的实现原理我们可以对采样序列进行逐次奇偶抽选,打乱以前的次序重新排序,然后按此顺序参加运算,可以实现"即位运算"提高存储单元的利用率。

  (一) 复数的描述方法

  进行傅立叶变换时不可避免的要用到复数,而在VC中并没有现成的可用于表示复数的数据类型,可以自己定义一个含有两个成员变量的数据结构来表示复数,这两个成员变量可分别用于表示复数的实部与虚部:

typedef struct tagComplex{
 float Re; //复数的实部
 float Im; //复数的虚部
}Complex;

  (二) 倒序的实现

  在进行快速傅立叶变换时,可以将输入的时域序列和输出的频域序列都按照自然顺序排列;也可以按照"蝴蝶图"所描述的计算方法对输入的时域序列按奇偶分解后的序列排序而输出的频域序列仍是按自然顺序排列的;还有一中方式是输入的时域序列是不进行抽选的自然序列,而输出的频域序列则是按奇偶分解后的顺序排列的。这三种方式各有优点,第一种对输入、输出不需要进一步排序,但由于自然排序不符合"蝴蝶图"运算规律,会占用大量中间存储单元。而后两种则无须中间存储单元,但需要倒一次序。权衡利弊,当采样点较多时还是采用后两种方式好,多一次倒序运算对现在的高性能计算机而言并不是什么负担。下面代码用于对原始采样序列的时间抽选奇偶分解工作,其中A、N分别表示指向采样序列复数数组的指针和序列的长度。

int NV2=N/2;
int NM1=N-1;
int I,J,K=0;
Complex T;//用于中介的复数变量T
I=J=1;
while(I<=NM1)
{
 if(I {
  T=A[J-1];//将A[J-1]的内容和A[I-1]的内容互换,借助于中间变量T
  A[J-1]=A[I-1];
  A[I-1]=T;
 }
 K=NV2;
 while(K {
  J-=K;
  K/=2;
 }
 J+=K;
 I++;
}


 

<

 

 

 
说明
:本教程来源互联网或网友上传或出版商,仅为学习研究或媒体推广,wanshiok.com不保证资料的完整性。
上一篇:VC++实现动画弹出/弹入式窗口  下一篇:利用消息机制实现VC与Delphi之间的通讯