存中，并将其转化为复变量存放于复变量数组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算子的值
W.Re=(float)cos(PI/(1.0＊LE1));
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;
//复数运算A×U
T.Re=(float)A[IP－1].Re＊U.Re－A
[IP－1].Im＊U.Im;
T.Im=(float)A[IP－1].Re＊U.Im＋A
[IP－1].Im＊U.Re;
//复数运算A－T
A[IP－1].Re=(float)A[I－1].Re－
T.Re;
A[IP－1].Im=(float)A[I－1].Im－
T.Im;
//复数运算A＋T
A[I－1].Re＋=T.Re;
A[I－1].Im＋=T.Im;
I＋=LE;
}
float temp=U.Re;
//复数运算U×W
U.Re=(float)U.Re＊W.Re－U.Im＊
W.Im;
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
……
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方法的需要，先对波形进行了抽样和截断，这样再用
程序分析采样数据必然会引起误差。从分析结果还可以看出，频率越高，误差波动也越大，此分析结果产生的误差在允许范围之内，是一个可以允许的近似。
本
程序在Windows 98、Microsoft Visual C＋＋ 6.0下编译通过。
