对波形数据的频谱分析

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

2、换需要同一个 NXN点的W矩阵(关于W矩阵 请参阅信号与系统方面或数学方面的书籍) 相运算，随着N值的增大，运算次数 显著上 升，当点数达到1024时，需要进行复数乘法运算1048576次。显然这种 算法在实际运用中无法保证当点数较大时的运算速度， 无法满足对信号的实时处 理要求。根据W矩阵中顿素的周期性和对称性我们可以将一个 N点的DFT!算分解为 两组N/2点的DFT运算，然后取和即可。为进一步提高效率，将上述两个矩阵按 奇偶顺序逐级分解下去。当采样点数为2的指数次方M时，可分解为M级子矩阵 运算，全部工作量如下：复数乘法：M< N/2次复数加法：NX M次直接采用DFT算法需要的运

3、算量为:复数乘法：NX N次复数加法：NX (N- 1)次当点数N为几十个点时快速傅立叶交换(FFT)的优势还不明显，而一旦 N达到 几千时优势是十分明显的：N=1024时：DFT需1048576次运算，FFT仅需5120次运算，改善比为204.8。N=2048时：DFT需4194304次运算，FFT仅需11264次运算，改善比达到 372.4。当采样点数较多时，如变换前和变换后的序列都按自然顺序排列，则中间运算 过程会占用大量的中间存储单元，造成效率的低下和存储单元的浪费。根据FFT的实现原理我们可以对采样序列进行逐次奇偶抽选，打乱以前的次序重新排序， 然后按此顺序参加运算，以“即位运算”提

4、高存储单元的利用率。复数的描述方法进行傅立叶变换时不可避免地要用到复数，而在 VC中并没有现成的可用于表 示复数的数据类型，因此需要自己定义一个含有两个成员变量的数据结构来表示 复数，这两个成员变量可分别用于表示复数的实部与虚部：typedef struct tagComplex(/复数的实部float Re;/复数的虚部float Im;Complex;倒序的实现在进行快速傅立叶变换时，可以将输入的时域序列和输出的频域序列都按照自 然顺序排列；也可以按照“蝴蝶图”所描述的计算方法对输入的时域序列按奇偶 分解后 的序列排序，而输出的频域序列仍是按自然顺序排列；还有一种方式是 输入的时域序列是自

5、然序列，而输出的频域序列则是按奇偶分解后的顺序排列。 这三种方式各 有优缺点：第一种对输入、输出不需要进一步排序，但由于自然 排序不符合“蝴蝶图”运算规律，会占用大量中间存储单元；而后两种则无需中 间存储单元，但需要 倒序。权衡利弊，当采样点较多时还是采用后两种方式好， 多一次倒序运算对现在的高性能计算机而言并不是什么负担。下面代码用于对原 始采样序列的时间抽选奇 偶分解工作，其中A、N分别表示指向采样序列复数数 组的指针和序列的长度。int NV2=N/2;int NM1=N 1;int I,J,K=0;/用于中介的复数变量TComplex T;I=J=1;while(I<=NM1)(

6、if(I< J)(/借助于中间变量T,将AJ-1的内容和AI -1的内容互换T=AJ 1;AJ 1=AI 1;AI 1=T;K=NV2;while(K< J)(J =K;K/=2;J + =K;I + ;时域信号的频谱分析首先要将从外设输入或采集的时域波形数据经抽样量化后,通过CFile类的Open (, )、Read(,)等成员函数将其读取到缓存中，并将其转化为复变量存放于复变量数组A中。同时需要验证数据量的长度是否为 2的整数次籍，如 不是则用0来补齐，否则无法用“蝴蝶图”进行分解运算。下面代码用于完成对原始采样时域序列的快速傅立叶变换，A、M分别表示指向原始采样数据数组的指针

7、和序列长度的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;i

8、f(abs(W.Im)< 1.0e12)W.Im=0.0f;J=1;while(J<=LE1)(I=J;while(I<=N)(IP=I + LE1;/复数运算AX UT.Re=(float)AIP 1.Re 大 U.Re AIP 1.Im 大 U.Im;T.Im=(float)AIP 1.Re 大 U.Im + AIP 1.Im 大 U.Re;/复数运算A TAIP 1.Re=(float)AI 1.Re -T.Re;AIP 1.Im=(float)AI 1.Im 一T.Im;/复数运算A+ TAI 1.Re + =T.Re;AI 1.Im + =T.Im;I + =LE;

9、float temp=U.Re;/复数运算UX WU.Re=(float)U.Re 大 W.R U.Im 大W.Im;U.Im=(float)temp大 W.Im+ U.Im *W.Re;J Hk ;L + ;,上述代码执行完毕时，原先存放着时域数值的复变量数组内存放的就是经过分析后的频域值，利用此数据可以通过绘图将频域波形直观地显示出来，也可以将其存成数据文件，以备进一步使用。测试及运算结果分析编译运行程序，分析一个三角脉冲的数据文件，并保存分析结果。该三角脉冲 幅度为1,持续时间2毫秒，抽样时间间隔是20微秒，延拓周期(数据记录长 度)为10毫秒，采样点数为500,取2的整数次籍512个采

10、样点。下附该三角 脉冲频谱的计算结果及误差分析：频率(Hz) FFT求得X(f) 误差0.00 1.00006E 03 1.00000E 03 6.10352E 08100.00 9.67593E 04 9.67531E 04 6.14332E 08200.00 8.75203E 04 8.75150E 04 6.25092E 08300.007.36904E一047.36849E一046.39413E一08400.005.72852E一045.72787E一046.52926E一08500.004.05351E一044.05285E一046.61362E一08600.002.54638E一0

11、42.54572E一046.61847E一082700.009.16539E一069.09679E一066.86075E一082800.004.53216E一064.46500E一066.71550E一082900.001.21487E一061.15945E一066.44190E一08注：此处FFT运算结果都乘以了系数10毫秒（0.01秒）。从上述数据中可以看出，在分析结果中产生了误差。这是由丁待分析的连续时 间信号不具备离散性或周期性，也可能有无限长度。为了适应FFT方法的需要,先对 波形进行了抽样和截断，这样再用程序分析采样数据必然会引起误差。从 分析结果还可以看出，频率越高，误差波动也越

12、大，此分析结果产生的误差在允 许范围之 内，是一个可以允许的近似。本程序在 Windows 98、Microsoft Visual C + + 6.0 下编译通过Top2楼MPU（黑眼圈）回复于 2003-11-29 21:51:44 得分0pDoc->myterms是迭代次数.差不多1000就成为正方形.IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII/II CFourierView drawingvoid CFourierView:OnDraw(CDC* pDC)(CFourierDoc* pDoc

13、 = GetDocument();ASSERT_VALID(pDoc);/ all remaining code for Fourier Series int i,j,ang; double y, yp; CBrush newbrush; CBrush* oldbrush; CPen newpen; CPen* oldpen;/ common color dialog box structure information / allow initial color value to be setdlg1.m_cc.Flags |= CC_FULLOPEN | CC_RGBINIT; dlg1.m

14、_cc.rgbResult = pDoc->mycolor;pDC->SetMapMode(MM_ISOTROPIC); pDC->SetWindowExt(500,500);pDC->SetViewportExt(m_cxClient,-m_cyClient); pDC->SetViewportOrg(m_cxClient/20,m_cyClient/2);ang=0;yp=0.0;newpen.CreatePen(BS_SOLID,1,RGB(0,0,0); oldpen=pDC->SelectObject(&newpen);/ draw x &

15、amp; y coordinate axes pDC->MoveTo(0,240);pDC->LineTo(0,-240);pDC->MoveTo(0,0);pDC->LineTo(400,0);pDC->MoveTo(0,0);/ draw actual Fourier waveform for (i=0; i<=400; i+) for (j=1; j<=pDoc->myterms; j+) y=(250.0/(2.0*j)-1.0)*sin(j*2.0)-1.0)*(ang*2.0*3400.0);yp=yp+y; pDC->LineTo(i,(int) yp); yp-=yp; ang+;/ create brush from common color dialog box select ion/ for waveform fillnewbrush.CreateSolidBrush(pDoc->myc