郎锐

频谱分析是电子工程上一个非常重要的分析手段，许多计算机辅助电路分析（CAA）类软件都具备这种分析能力，以便电子工程师能清楚地看到某波形的频谱分布情况。要对一个输入信号源作频谱分析，将其由时域信号转变为频域信号，就必然要用到傅立叶变换。这样，无论是在时域还是在频域，都要对连续函数进行积分运算。很显然，要通过计算机实现这种变换就需要预先通过抽样将原始的连续数据转变为离散数据，并将计算范围收缩到一个有限区间。因此，在允许一定程度近似的条件下，可以使用“离散傅立叶变换（DFT）”对波形数据进行频谱分析。  
算法构成原理
要计算一个N点的离散傅立叶变换需要同一个N×N点的W矩阵（关于W矩阵请参阅信号与系统方面或数学方面的书籍）相运算，随着N值的增大，运算次数显著上升，当点数达到1024时，需要进行复数乘法运算1048576次。显然这种算法在实际运用中无法保证当点数较大时的运算速度，无法满足对信号的实时处理要求。
根据W矩阵中W元素的周期性和对称性我们可以将一个N点的DFT运算分解为两组N/2点的DFT运算，然后取和即可。为进一步提高效率，将上述两个矩阵按奇偶顺序逐级分解下去。当采样点数为2的指数次方M时，可分解为M级子矩阵运算，全部工作量如下：
复数乘法：M×N/2次
复数加法：N×M次
直接采用DFT算法需要的运算量为：
复数乘法：N×N次
复数加法：N×(N－1)次
当点数N为几十个点时快速傅立叶交换(FFT)的优势还不明显，而一旦N达到几千时优势是十分明显的：
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;
//用于中介的复数变量T
Complex T;
I=J=1;
while(I<=NM1)
{
if(I< J)
{
//借助于中间变量T，将A[J－1]的内容和A[I－1]的内容互换
T=A[J－1];
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(……)等成员函数将其读取到缓