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

　   摘要： 本文介绍了采用离散傅立叶变换（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<J)
　 {
　　 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)
　 {
　　 J-=K;
　　 K/=2;
　 }
　 J+=K;
　 I++;
}
　　（三） 时域信号的频谱分析

　　首先要将从外设输入或采集的时域波形数据经抽样量化后，通过CFile类的Open（……）、Read(……)等成员函数将其读取到缓存中，并将其转化为复变量存放于复变量数组A中，同时需要验证以下数据量的长度是否为2的整数次幂，如若不是则必须用0来补齐，否则无法用"蝴蝶图"进行分解运算。下面代码用于完成对原始采样时域序列的快速傅立叶变换，A、M分别表示指向原始采样数据数组的指针和序列长度的2的整数次幂：

……
Complex U,W,T;
int LE,LE1,I,J,IP;
int N=(int)pow(2,M);
//在此采用的是时间抽选奇偶分解方式，所以在参加运算前首先要对时间序列进行倒序
ReverseOrder(A,N);
int L=1;
while(L<=M)
{
　 LE=(int)pow(2,L);
　 LE1=LE/2;
　 U.Re=1.0f;
　 U.Im=0.0f;
　 W.Re=(float)cos(PI/(1.0*LE1));//计算W算子的值
　 W.Im=(float)-1.0*sin(PI/(1.0*LE1));
　 if(abs(W.Re)<1.0e-12)
　　 W.Re=0.0f;
　 if(abs(W.Im)<1.0e-12)
　　 W.Im=0.0f;
　　 J=1;
　　 while(J<=LE1)
　　 {
　　　 I=J;
　　　 while(I<=N)
　　　 {
　　　　 IP=I+LE1;
　　　　 T.Re=(float)A[IP-1].Re*U.Re-A[IP-1].Im*U.Im;//计算复数运算A*U
　　　　 T.Im=(float)A[IP-1].Re*U.Im+A[IP-1].Im*U.Re;
　　　　 A[IP-1].Re=(float)A[I-1].Re-T.Re;//计算复数运算A-T
　　　　 A[IP-1].Im=(float)A[I-1].Im-T.Im;
　　　　 A[I-1].Re+=T.Re;//计算复数运算A+T
　　　　 A[I-1].Im+=T.Im;
　　　　 I+=LE;
　　　 }
　　　 float temp=U.Re;
　　　 U.Re=(float)U.Re*W.Re-U.Im*W.Im;//计算复数运算U*W
　　　 U.Im=(float)temp*W.Im+U.Im*W.Re;
　　　 J++;
　　 }
　　 L++;
　 }
……
　　上述代码执行完毕时，原先存放着时域数值的复变量数组内存放的就是经过分析后的频域值了，对此数据可以通过绘图将频域波形直观的显示出来，也可以将其存成数据文件，以备进一步使用。

　　四、 测试及运算结果分析

　　编译运行程序，打开一三角脉冲的数据文件，并将分析结果保存，该三角脉冲幅度为1，持续时间2毫秒，采样时抽样时间间隔是20微秒，延拓周期（数据记录长度）为10毫秒，采样点数目500点，取2的整数次幂512个样点。下附该三角脉冲频谱的计算结果及误差分析：

　　频率（Hz） FFT求得　　　　 X(f)　　　　 误差
　　0.00 　　1.00006E-03 　1.00000E-03　 6.10352E-08
　　100.00 　9.67593E-04 　9.67531E-04　 6.14332E-08
　　200.00 　8.75203E-04 　8.75150E-04　 6.25092E-08
　　300.00 　7.36904E-04 　7.36849E-04　 6.39413E-08
　　400.00 　5.72852E-04 　5.72787E-04　 6.52926E-08
　　500.00 　4.05351E-04 　4.05285E-04　 6.61362E-08
　　600.00 　2.54638E-04 　2.54572E-04　 6.61847E-08
　　700.00 　1.35403E-04 　1.35338E-04　 6.53870E-08
　　800.00 　5.47602E-05 　5.46963E-05　 6.39612E-08
　　900.00 　1.20072E-05 　1.19448E-05　 6.23453E-08
　　1000.00　 6.10719E-08　 1.17757E-32　 6.53870E-08
　　1100.00　 8.05672E-06　 7.99613E-06　 6.05985E-08
　　1200.00　 2.43706E-05　 2.43095E-05　 6.11450E-08
　　1300.00　 3.93026E-05　 3.92400E-05　 6.25965E-08
　　1400.00　 4.68226E-05　 4.67581E-05　 6.45128E-08
　　1500.00　 4.50979E-05　 4.50316E-05　 6.62543E-08
　　1600.00　 3.58664E-05　 3.57992E-05　 6.71930E-08
　　1700.00　 2.30135E-05　 2.29466E-05　 6.69399E-08
　　1800.00　 1.08697E-05　 1.08042E-05　 6.55073E-08
　　1900.00　 2.74348E-06　 2.68014E-05　 6.33390E-08
　　2000.00　 6.11826E-08　 1.17757E-32　 6.11826E-08
　　2100.00　 2.25379E-06　 2.19395E-06　 5.98376E-08
　　2200.00　 7.29243E-06　 7.23256E-06　 5.98625E-08
　　2300.00　 1.25974E-05　 1.25360E-05　 6.13467E-08
　　2400.00　 1.59746E-05　 1.59107E-05　 6.38421E-08
　　2500.00　 1.62779E-05　 1.62114E-05　 6.64915E-08
　　2600.00　 1.36254E-05　 1.35571E-05　 6.83226E-08
　　2700.00　 9.16539E-06　 9.09679E-06　 6.86075E-08
　　2800.00　 4.53216E-06　 4.46500E-06　 6.71550E-08
　　2900.00　 1.21487E-06　 1.15945E-06　 6.44190E-08
　
　　注：在此，FFT运算结果都倍乘了系数10毫秒（0.01秒）。

　　在分析结果中产生了误差，是由于待分析的连续时间信号不具备离散性或周期性，也可能有无限长度。为了适应FFT方法的需要，对波形进行了抽样和截断，这样再用程序分析采样数据必然会引入误差，从分析结果可以看出，频率越高，误差波动也越大，此分析结果产生的误差在允许范围之内，是一个可以满意的近似。实践证明，本程序的算法是正确可靠的。

　　小结：

　　DFT尤其是FFT的应用已遍及各个科学领域，"DFT的应用"与 "FFT的应用"几乎成为同义语。通过本文介绍和程序示例可以清楚的看到FFT方法在直接处理离散信号数据的作用，而且也可以很好的用于对连续时间信号分析的逼近。本程序在Windows 2000 Professional下、由Microsoft Visual C++ 6.0编译通过。