风控建模招聘:用VC＋＋实现对波形数据的频谱分析

● 郎锐
01-7-4 上午 10:47:00
频谱分析是电子工程上一个非常重要的分析手段，许多计算机辅助电路分析（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(……)等成员函数将其读取到缓存中，并将其转化为复变量存放于复变量数组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下编译通过。