用于确定非相干抽样数据的功率频谱的低泄露技术

摘要

本发明提供了一种技术，用于对从自动测试系统中抽样的波形的频率分量的幅度进行确定，该技术包括：对预计在抽样波形中发现的N个频率的列表进行汇编。在测试仪上运行的测试程序通常提供频率列表。该技术假定抽样波形与理想化波形模型一致，而该理想化波形模型在数学上与N个正弦波之和对应。构成该模型的N个正弦波中的每个均具有未知幅度，并具有一频率，该频率与频率列表中的N个频率中的一个相等。该技术试图在数学上通过线性最小二乘方算法使模型与实际抽样波形之差最小，从而求出N频率中的各方的未知幅度。

说明书

发明领域

本发明涉及电子装置测试设备和方法，具体涉及用于从电子测试信号中抽出测试信号频率分量的幅度的技术。

发明的背景

自动测试系统的测试程序一般要求测试仪对从在测装置(DUT)抽样的信号的功率频谱进行测定。在常规测试方案中，自动测试系统生成用于执行DUT输入的激励，并随着DUT响应该激励，对DUT的输出进行抽样。测试仪软件通过对所获得的抽样进行离散傅里叶变换(DFT)来计算抽样输出信号的功率频谱。

众所周知，每当抽样时钟与抽样信号不“相干”时，被称为“泄漏”的误差就会显现在由DFT生成的功率频谱中。如果抽样时钟的频率是抽样信号内存在的各频率的精确整倍数，则抽样时钟是“相干的”。泄漏是对截断频率，即：在抽样窗口内未完成整个循环的频率执行DFT的数学推论。随着谱线的误展宽、假波峰和假波谷(凸起部)的生成、以及功率频谱噪声最低限度的一般升高，会观察到泄漏。

已设计了多种方法来减少泄漏。一种方法是增加抽样速度。一般来说，抽样速度越高，感兴趣频率范围内的截断量就越小，泄漏误差就越小。该方法尽管有效，然而增加抽样速度来减少泄漏仅与这种增加的幅度成正比。而且该方法往往会大幅增加所用抽样设备的成本。

用于减少泄漏的另一常用技术是用开窗口函数来乘以抽样的数据序列。窗口函数具有使抽样数据序列围绕其端点递减的效果，从而消除会引起泄漏误差的不连续性。可使用不同的窗口函数，例如，Blackman，Hanning，或Hamming窗口函数，每种函数均具有其自身具体特性。窗口函数往往会减少与功率频谱中的波峰相隔一定距离的泄漏偏差，但也往往会生成较宽波峰。这样，这些窗口函数具有重新分配泄漏而不是完全消除泄漏的效果。而且，由于窗口函数实际改变执行DFT的数据，因而这些窗口函数往往会使频谱略微失真。

还有一种技术是按照与抽样信号频率相干的抽样速度对波形数据“重新抽样”。重新抽样工作是通过在按照某一速度抽样的实际点之间进行插值来实现的，以便在数学上构建一系列点，该一系列点看起来已按照某一不同速度进行抽样。尽管重新抽样对于减少泄漏会极其有效，然而它的计算工作量大，并且其精度会受到插值误差的影响。

还有一种减少泄漏的技术是改变抽样时钟速度，以使该速度与在抽样信号中发现的各频率的整倍数精确相等。该技术虽极其有效，但需要价格昂贵的硬件。当测试仪包括大量抽样时钟时，该方案尤其价格昂贵，而情况经常是这样。

自动测试设备(ATE或“测试仪”)的制造商一般通过提供价格不贵的常规测试问题解决方案来力图改善其产品。通过增加测试仪性能，同时降低测试仪成本，可获益匪浅。为此，目前强烈要求提供一种价格低廉的技术来减少由自动测试系统抽样的信号频谱中的泄漏。

发明概述

考虑到上述背景，本发明的目的是减少抽样信号中的泄漏，而无需大幅增加测试仪成本。

为了实现上述目的以及其他目的和优点，提供了一种用于对抽样波形频率含量进行分析的技术，该技术包括对预计在抽样波形中将要求发现的N个频率的列表进行汇编。假定抽样波形与一波形模型一致，而该波形模型在数学上与N个正弦波之和对应。N个正弦波中的各方均具有未知幅度和相位，并具有一频率，该频率与频率列表中的N个频率中的一不同频率相等。该技术可求出使模型与抽样数据最佳拟合的未知幅度和/或相位。

根据一个说明性实施例，当抽样的波形频率事先未知时，也可使用上述技术。根据该变形例，可对抽样波形进行傅里叶变换计算，以生成常规(rough)功率频谱。对该常规功率频谱中的波峰进行识别，并对其频率进行编译，以形成N个频率的列表。在对N个频率列表进行编译时，可考虑其他因素，例如，向从中获得抽样数据的装置施加的已知激励，以及其他伴随情况。然后，对生成的N个频率列表实施上述技术，以确定N个正弦波中的各方的精确幅度和/或相位。

附图说明

通过考虑以下说明和附图，将明白本发明的其他目的、优点和新颖性特征，在附图中：

图1是一种激励在测装置的输入并从在测装置的输出中进行信号抽样的自动测试系统的高度简化示意图；

图2是对根据本发明的一种处理进行说明的高级流程图，用于当频率分量的频率事先已知时，对从图1的自动测试系统抽样的波形的频率分量的幅度进行确定；

图3是对根据本发明的一种处理进行说明的流程图，该处理用于当频率分量事先未知时，对由图1的自动测试系统抽样的波形的功率频谱进行计算；以及

图4A～4C是根据不同模拟误差对用于减少泄漏的不同技术的性能进行比较的模拟功率频谱。

优选实施例的详细说明

技术

图1是示出一种使用自动测试系统110对在测装置(DUT)120进行测试的常规构成的高度简化方框图。该自动测试系统包括主计算机112。该主计算机112包括用于运行测试程序的测试仪软件。该测试程序对用于测试DUT的测试仪硬件资源进行控制。例如，测试程序可控制频率合成器114，以便向DUT的输入施加激励，并可控制数字化仪116，以便在DUT的输出生成的响应进行抽样。

主计算机一般把由数字化仪116获得的抽样数据存储在存储器内进行分析。测试程序，或者可用于测试程序的软件例行程序可对存储的抽样数据进行操作，以分析该抽样数据内容。常规上，测试程序将引导测试仪软件对抽样的数据进行离散傅里叶变换(DFT)。然后，测试程序将对DFT的结果进行测试。

图2示出了根据本发明的用于在ATE环境中对测试信号进行抽样和分析的一般处理。在步骤210，自动测试系统110向DUT 120的输入施加激励。

在步骤212，访问N个频率列表。该N个频率表示幅度和/或相位信息期望是已知的抽样波形的频率。数字N可以是任何正整数。并不是所有的N个频率都实际需要存在于抽样波形中。实际上，该技术可用于测定任何具体频率分量的有无。优选地是，该N个频率列表事先是已知的，并存储在测试程序内。

在步骤214，对抽样波形进行计算机建模。表示抽样波形的该计算机模型由总计的N个正弦波构成，以近似于实际抽样波形。该N个正弦波中的每个均采用以下公式：

                A_ksin(ω_ki)+B_kcos(ω_ki)        (公式1)

式中：

·  “A_k”和“B_k”是未知系数，

·  “k”是范围从1到N的指数，并表示N个频率分量中的一个，

·  ω_k表示第k频率分量(具体地说，ω_k＝2πF_k，式中，F_k是第k频率)，以及

·  “i”是识别具体抽样并表示时间的指数。

尽管公式1看起来是两个正弦波之和，然而公式1在数学上等于单个正弦波，该正弦波的频率等于ω_k/2π，幅度等于 $>>>A>2>>k>+>>B>2>>k>>>，相位等于B$_k和A_k的2自变量反正切。

假定公式1表示抽样波形中的N个正弦波中的每个，则整个抽样波形可使用以下公式来建模： $>>>\Sigma >>k>=>1>>N>>>(>>A>k>>cos>>\omega >k>>i>+>>B>k>>sin>>\omega >k>>i>)>>>>(公式2)$

在步骤216，对公式2的波形模型进行计算机处理，以获得在模型和实际抽样波形之间的最佳拟合。优选实施例采用线性最小二乘方技术，以使模型与数据拟合。特别是，步骤216试图使以下最小二乘方估计量最小： $>>>\Sigma >>i>=>0>>M>>[>>y>i>>->>\Sigma >>k>=>1>>N>>>(>>A>k>>cos>>\omega >k>>i>+>>B>k>>sin>>\omega >k>>i>)>>>]>2>>>>(公式3)式中，y$_i是抽样波形的第i抽样点，并且i的范围从0到M，其中，M表示抽样波形中的抽样总数。

为了最小化公式3，该技术认识到，当公式3的偏导数等于零时，即：当根据A_k和B_k的各值获得公式3的偏导数时，可实现最佳拟合。由于N个频率中的每个均存在A_k和B_k的值，因而根据这些A_k和B_k的各值获得公式3的偏导数，可生成2N方程组。 $>>>\Sigma >>i>=>0>>M>>>y>i>>cos>>\omega >j>>i>=>>\Sigma >>k>=>1>>N>>>(>>A>k>>>\Sigma >>i>=>0>>M>>cos>>\omega >k>>i>cos>>\omega >j>>i>+>>B>k>>>\Sigma >>i>=>0>>M>>>>sin>\omega >>k>>i>cos>>\omega >j>>i>)>>>>(公式4)$ >>>\Sigma >>i>=>0>>M>>>y>i>>sin>>\omega >j>>i>=>>\Sigma >>k>=>1>>N>>>(>>A>k>>>\Sigma >>i>=>0>>M>>cos>>\omega >k>>i>sin>>\omega >j>>i>+>>B>k>>>\Sigma >>i>=>0>>M>>sin>>\omega >k>>i>sin>>\omega >j>>i>)>>>>(公式5)式中，当“j”指数的范围从1到N时，公式4和公式5均重复N次。定义以下系数可简化说明：$$

令 $>>>cc>kj>>=>>\Sigma >>i>=>0>>M>>cos>>\omega >k>>i>cos>>\omega >j>>i>>>(公式6)$

令 $>>>sc>kj>>=>>\Sigma >>i>=>0>>M>>sin>>\omega >k>>i>cos>>\omega >j>>i>>>(公式7)$

令 $>>>cs>kj>>=>>\Sigma >>i>=>0>>M>>cos>>\omega >k>>i>sin>>\omega >j>>i>>>(公式8)$

令 $>>>ss>kj>>=>>\Sigma >>i>=>0>>M>>sin>>\omega >k>>i>sin>>\omega >j>>i>>>(公式9)$

这些系数可在数学上简化是通过认识到以下公式： $>>>cc>kj>>=>>(>>\Sigma >>i>=>0>>M>>cos>>(>>\omega >k>>+>>\omega >j>>)>>i>+>>\Sigma >>i>=>0>>M>>cos>>(>>\omega >k>>->>\omega >j>>)>>i>)>>/>2>>>$ >>>sc>kj>>=>>cs>kj>>=>>(>>\Sigma >>i>=>0>>M>>sin>>(>>\omega >k>>+>>\omega >j>>)>>i>+>>\Sigma >>i>=>0>>M>>sin>>(>>\omega >k>>->>\omega >j>>)>>i>)>>/>2>>>$ >>>ss>kj>>=>>(>>\Sigma >>i>=>0>>M>>cos>>(>>\omega >k>>+>>\omega >j>>)>>i>->>\Sigma >>i>=>0>>M>>cos>>(>>\omega >k>>->>\omega >j>>)>>i>)>>/>2>>>以及$ >>>\Sigma >>i>=>0>>M>>cos>\alpha i>=>>(>cos>\alpha >>(>N>->1>)>>->cos>\alpha N>->cos>\alpha >+>1>)>>/>2>>(>1>->cos>\alpha >)>>>>和$ >>>\Sigma >>i>=>0>>M>>sin>\alpha i>=>>(>sin>\alpha >>(>N>->1>)>>->sin>\alpha N>+>sin>\alpha >)>>/>2>>(>1>->cos>\alpha >)>>>>$$$$$

式中，α是ω的任何任意值。

使用在公式6～公式9中定义的系数来重写公式4和公式5，可生成以下矩阵：(公式10)通过确定矩阵X的逆矩阵，并把该逆矩阵乘以公式10左侧的矢量V，可解公式10求矢量ab中的A_k和B_k的每个值。

一旦对从1到N的k的每个值的A_k和B_k是已知的，就可通过计算 $>>>A>2>>k>+>>B>2>>k>>>来确定频率列表的每个第k个频率的幅度。通过计算B$_k和A_k的2自变量反正切，可确定各相位的值。

通过施加约束，即：2N＝M+1(频率数是抽样数的一半)，可减轻一些通过解公式10所施加的计算负担。通过把该约束放在合适位置，可将公式10的矢量V重写如下：(公式11)由于施加上述约束会迫使C和X^-1成为同一秩(rank)的方形矩阵，因而可把公式10和公式11进行组合以形成：

ab＝(X^-1C)y                      (公式12)

这样，可求出ab，而无须计算V。

建议使用以下方法求出ab：

·  首先，使用以下递归关系来快速构建C：

cos(a+1)ω_k＝2cosω_kcosaω_k-cos(a-1)ω_k

sin(a+1)ω_k＝2cosω_ksinaω_k-sin(a-1)ω_k

·  然后，通过L-U分解，计算X^-1；

·  把X^-1施加于C；

·  通过使X^-1C乘以y来计算ab。

一旦构建了X^-1C，计算ab就可大致要求N²倍增累积。

优选地是，本文所述技术是作为软件库中的一种函数来实施的。该函数优选地接收输入阵列，该输入阵列存储频率列表和所述抽样数据指针。该函数优选地返回包含A_k和B_k各值的阵列，从中可计算幅度和相位。或者，该函数直接返回幅度和相位。软件库优选地驻留在自动测试系统上，在该自动测试系统中，可访问在测试系统中运行的测试程序。

实施例

图4A～4C示出了根据本发明的最佳拟合技术的性能与用于减少泄漏的其他技术的关系曲线的模拟预测。在说明书的末尾提供了一份用于生成这些图表的数据的测试代码列表。图4A～4C中的三个图中的每个均对在四种不同条件下获得的单音调功率频谱作了比较：

1.对未校正数据(即：使用矩形窗口获得的，标为“未校正”)进行快速傅里叶变换(FFT)，

2.对由Hanning窗口成形的数据(“已开窗口”)进行FFT，

3.对已重新抽样数据(“已重新抽样”)进行FFT，以及

4.本文所述的最佳拟合技术(“最小二乘方”)。

各图的水平轴表示频率，具体地说，表示频率点(Bin)0～63。为了与采用FFT的方法进行直接比较，使用N＝64来操作最佳拟合技术，其中，N个频率中的每个均与FFT点(Bin)对应。垂直轴表示以db为单位的幅度。

各图均示出了从第9频率点(Bin)即：即未进行相干抽样的频率点(Bin)的中心略微偏斜的单音调。特别是，ω_k＝2π(k-1)(1+ε)，式中，对图4A中ε＝10^-6，对图4B为ε＝10^-9，对图4C为ε＝10^-12。这些图均表明了最佳拟合技术的强度。与其他技术相比，最佳拟合技术可保持极窄波峰，而在波峰周围没有升高区域(“裙部”)。

优点

所揭示技术与常规DFT比起来有许多优点，特别是在自动测试设备方面。使用所述技术，可利用价格较低廉的测试仪电子装置进行高精度频谱分析。抽样时钟无需与待测定频率相干，并可大大消除频谱泄漏。与把频率分类为宽度有限的各点(Bin)的常规DFT相比，所揭示技术不采用频率点(Bin)，而采用离散频率。这样，所揭示技术可按照使用常规DFT无法实现的方式来分辨间隔极密的频率。

该技术还具有可升级的，因为其计算时间随着N的函数，即：待分析频率数而变化。这样，如果仅分析少量频率，则可较快实施该技术。此外，由于该技术采用最佳拟合算法，因而可使用该技术来确定在抽样窗口内未充分完成一个循环的频率分量的幅度和相位。要求事先规定频率，一般在自动测试设备中不会起到不利作用，在该自动测试设备中，由DUT生成的频率事先在很大程度上是已知的，并且测试仪生成用于驱动DUT的激励。

尽管X^-1C需要大量时间来计算，然而每当分析波形时，无需对其重新计算。只要频率列表和抽样数M保持恒定，就可通过检索X^-1C的存储拷贝并将其乘以y，对新获得的抽样数据进行分析。提供的库可包括针对不同频率和抽样数的各种不同X^-1C组合。用户可从这些组合中进行选择，以便对波形进行快速分析。

替代例

在对一个实施例作了说明之后，可进行各种替代实施例或变化。

例如，如上所述，在自动测试设备方面可使用抽样波形分析技术。然而，该技术可更普遍地适用于期望进行频率信息分析的任何抽样数据。尽管根据图1的具体测试方案对该技术作了说明，然而该技术不限于任何具体测试方案。

本文揭示的实施例对使用线性最小二乘方来获得在波形模型(公式2)和实际抽样数据之间的最佳拟合作了规定。然而，可使用其他最佳拟合技术，例如，柯西-洛伦兹(Cauchy-Lorentz)分布，以及试图使模型和抽样数据之差的绝对值最小的技术。因此，本发明不限于使用最小二乘方。

此外，本说明还提供了用于处理矩阵和执行必要计算的软件。或者，可提供专用硬件电路或处理器，以便更有效地执行这些功能。

如上所述，该技术要求N个频率列表事先是已知的。然而，可通过对抽样数据进行DFT并检查结果来避免该要求。这种变形例如图3所示。在步骤310，进行DFT。在步骤312，在从DFT得出的功率频谱中识别波峰。然后，把与波峰对应的频率附加给频率列表，用于进行更准确的分析。步骤316和318如上所述进行：假定抽样数据与模型一致，并且获得在抽样和模型之间的最佳拟合。该技术无需在空白中(vacuum)进行。也可将伴随情况加以考虑，例如，施加给DUT的激励频率，这些频率的谐波，以及DUT的已知特性。

上述技术假定按照匀速进行波形抽样。然而，根据替代实施例，也可按照非匀速进行波形抽样。特别是，用“t_i”项(即：实际抽样时间)替代上述方程和矩阵中的离散指数“i”，可实现任意非均匀抽样。对于非均匀抽样，可以不使用公式9之后的计算简化；然而，该变化对所述技术的其余部分是透明的。

这些替代例和变形例以及其他例等均由本发明人作了设想，并用于包含在本发明的范围内。因此，应理解的是，上述说明仅是作为例子，并且本发明应仅由所附权利要求的精神和范围来限制。

计算机列表

以下提供了用于生成这些图表的数据的测试代码的软件列表：

 /* standard header files */ #include ＜stdio.h＞ #include ＜math.h＞ #include ＜sys/time.h＞ #include ＜image.h＞ #define N 128/* Number of samples of original waveform*/ #define M 64 /* Number of frequency bins for least-squares fit routine*/ #define M2 (2*M) /* Total number of bins(real，and imaginary) */ #define R 32 /* Number of samples to expand by in resampling routine*/ #define F 1024 /* Number of samples infilter in resampling routine*/ doubleraw_capture[N]；/* The original samplas*//* These are used in the resampling routine only */double extended_capture[2*N]，/* basically raw_capture[]，but with extra samples*/ expanded_capture[R*N]，/* this is (hopefully) equivalent to raw_capture[]*/　　/*but with R times as many samples*/ filter[1024]， /* uses a cos^2 filter 1024 points*/ resampled_capture[N]； /* this is (hopefully) raw_capture[] but sampled at */　　/*the correct rate*/double Mx[M2][M2]，/* This corresponds to matrix M */ Mi[M2][M2]，/* inverse of matrix M*/ MiC[M2][M2]， /* This is the product of M^-1C */ x[M2]，xmags[M]，xargs[M]， /* uncorrected dft-mags and phases (arguments)*/ xw[M2]，xwmags[M]，xwargs[M]，/* uncorrected dft with hamming window*/ xs[M2]，xsmags[M]，xsargs[M]，/* resampled dft*/ xp[M2]，xpmags[M]，xpargs[M]；/* least squares dft*//* Computes sum(cos w*i)i＝0 to n */double sumcos(w，n)double w；{return w ？ 0.5*(cos(w*(n-1))-cos(w*n)-cos(w)+1)/(1-cos(w)：n；}/* Computes sum(sin w*i) i＝0 to n */double sumsin(w，n)double w；{return w ？ 0.5*(sin(w*(n-1))-sin(w*n)+sin(w))/(1-cos(w))：0；}/* Makes 1024 point cos^2 filter-used for resampling algorithm *//* This is not used by least squares */make_filter(){int i； for(i＝0；i＜1024；i++)<!-- SIPO <DP n="12"> --><dp n="d12"/> filter[i]＝(1-cos(2+M_PI*i/1024))/2； } /* 1d interpolation routine-used for resampling algorithm.*/ /* Basically Laqrange’s interpolation */ /* This is not used by least squares */ /* tries to return y[a] given y[0]...y[n-1] as input */ double interpolate_1d(a，y，n) double a，y[]； int n； {int k，j，xl，xr；double product，sum＝0；xl＝ftoi(a-2.0)；xr＝ftoi(a+2.0)；xl＝xl＞＝0？xl：0；xr＝xl＜＝n？xr：n；for(k＝xl；k＜＝xr；k++){product＝y[k]；　　for(j＝xl；j＜＝xr；j++)　　 if(j！＝k)　　product*＝(a-j)/(k-j)； sum+＝product；}return sum；｝main(){int i，j，k，m； struct timeval tp0，tp1，tp2，tp3； /* Used for timing */ double f0＝1.0+0.001e-9； /* Frequency error of 0.001ppb */ gettimeofday(&tp0，NULL)； for(i＝0；i＜N；i++)raw_capture[i]＝cos(2*M_PI*9*f0*i/N)； gettimeofday(&tp1，NULL)； printf(″Time to construct raw capture：％6uus\n″，diff(tp1，tp0))； gettimeofday(&tp0，NULL)； for(m＝0；m＜M；m++) {x[2*m]＝x[2*m+1]＝0；for(i＝0；i＜N；i++){x[2*m]+＝raw_capture[i]*cos(2*M_PI*m*i/N)； x[2*m+1]+＝raw_capture[i]*sin(2*M_PI*m*i/N)； }} gettimeofday(&tp1，NULL)； printf(″Time to construct ％i-point uncorrected dft：％6uus\n″，M，diff(tp1，tp0))； for(i＝0；i＜M；i++) {xmags[i]＝hypot(x[2*i+1]，x[2*i+0])/M；xargs[i]＝180/M_PI*atan2.(x[2*i+1]，x[2*i+0])； ｝/* This section is not least squares fit，but is a crude resampling algorithm */#ifdef R gettimeofday(&tp0，NULL)； for(i＝0；i＜2*N；i++)extended_capture[i]＝cos(2*M_PI* 9*f0*(i-N/2)/N)； gettimeofday(&tp1，NULL)； printf(″Time to construct raw capture(with N/2 padding)：％6uus\n″，diff(tp1，tp0))； gettimeofday(&tp0，NULL)； make_filter()；<!-- SIPO <DP n="13"> --><dp n="d13"/> gettimeofday(&tp1，NULL)； printf(″Time to construct cos^2 filter：％6uus\n″，diff(tp1，tp0))； gettimeofday(&tp0，NULL)； for(i＝0；i＜N*R；i++) {int e0＝N/2-F/R/2+(i+R-1)/R；expanded_capture[i]＝0；for(j＝0；j＜F/R；j++){int fi＝(N*R-i)％R+j*R； int ej＝e0+j； expanded_capture[i]+＝filter[fi]*extended_capture[ej]； }} gettimeofday(&tp1，NULL)； printf(″Time to expand padded capture：％6uus\n″，diff(tp1，tp0))； gettimeofday(&tp1，NULL)； for(i＝0；i＜N；i++)resampled_capture[i]＝interpolate_ld(i*R/f0，expanded_capture，N*R)； gettimeofday(&tp1，NULL)； printf(″Time to resample from expanded capture：％6uus\n″，diff(tp1，tp0))； gettimeofday(&tp0，NULL)； for(m＝0；m＜M；m++) {xs[2*m]＝xs[2*m+1]＝0；for(i＝0；i＜N；i++){xs[2*m]+＝resampled_capture[i]*cos(2*M_PI*m*i/N)； xs[2*m+1]+＝resampled_capture[i]*sin(2*M_PI*m*i/N)； }} gettimeofday(&tp1，NULL)； printf(″Time to construct ％i-point resampled dft：％6uus\n″，M，diff(tp1，tp0))； for(i＝0；i＜M；i++) {xsmags[i]＝hypot(xs[2*i+1]，xs[2*i+0])/M；xsargs[i]＝180/M_PI*atan2(xs[2*i+1]，xs[2*i+0])； ｝#endif/* This section is not least squares fit，but instead uses a hanning window */ gettimeofday(&tp0，NULL)； for(m＝0；m＜M；m++) {xw[2*m]＝xw[2*m+1]＝0；for(i＝0；i＜N；i++){xw[2*m]+＝raw_capture[i]*cos(2*M_PI*m*i/N)*(1-(1+cos(2*M_PI*i/N))/2)； xw[2*m+1]+＝raw_capture[i]*sin(2*M_PI*m*i/N)*(1-(1+cos(2*M_PI*i/N))/2)； }} gettimeofday(&tp1，NULL)； printf(″Time to construct ％i-point windowed dft：％6uus\n″，M，diff(tp1，tp0))； for(i＝0；i＜M；i++) {xwmags[i]＝hypot(xw[2*i+1]，xw[2*i+0])/M；xwargs[i]＝180/M_PI*atan2(xw[2*i+1]，xw[2*i+0])； ｝ /* The least-squares section */ gettimeofday(&tp0，NULL)； for(i＝0；i＜M；i++)for(j＝0；j＜M；j++){Mx[2*i][2*j]＝0.5*(sumcos(2*M_PI*(i+j)*f0/N，N)+sumcos(2*M_PI*(i-j)*f0/N，N))； Mx[2*i][2*j+1]＝0.5*(sumsin(2*M_PI*(i+j)*f0/N，N)-sumsin(2*M_PI*(i-j)*f0/N，N))； Mx[2*i+1][2*j]＝0.5*(sumsin(2*M_PI*(i+j)*f0/N，N)+sumsin(2*M_PI*(i-j)*f0/N，N))； Mx[2*i+1][2*j+1]＝0.5*(-sumcos(2*M_PI*(i+j)*f0/N，N)+sumcos(2*M_PI*(i-j)*f0/N，N))； ｝Mx[1][1]＝Mx[0][0]；/* fix that nasty problem with sin */gettimeofday(&tp1，NULL)；printf(″Time to construct matrix：％6uus/n″，diff(tp1，tp0))；<!-- SIPO <DP n="14"> --><dp n="d14"/> gettimeofday(&tp0，NULL)； invert2M( )； gettimeofday(&tp1，NULL)； printf(″Time to invert matrix：％6uus\n″，diff(tp1，tp0))； gettimeofday(&tp0，NULL)； for(i＝0；i＜M2；i++)for(j＝0；j＜M2；j++){MiC[i][j]＝0.0； for(k＝0；k＜M；k++)　　MiC[i][j]+＝Mi[i][2*k]*cos(2*M_PI*k*f0*j/N)+　　Mi[i][2*k+1]*sin(2*M_PI*k*f0*j/N)；} gettimeofday(&tp1，NULL)； printf(″Time to construct composite matrix：％6uus\n″，diff(tp1，tp0))； gettimeofday(&tp0，NULL)； for(i＝0；i＜M2；i++) {xp[i]＝0；for(j＝0；j＜M2；j++) xp[i]+＝MiC[i][j]*raw_capture[j]； } gettimeofday(&tp1，NULL)； printf(″Time to multiply composite matrix on result：％6uus\n″，diff(tp1，tp0))； for(i＝0；i＜M；i++) {xpmags[i]＝hypot(xp[2*i+1]，xp[2*i+0])；xpargs[i]＝180/M_PI*atan2(xp[2*i+1]，xp[2*i+0])；}}/* These matrix inversion routines should be replaced by L-U decomposition*//* These unfortunately use a crude Gauss-Jordon algorithm with no pivoting */initMi( ){int i，j； for(i＝0；i＜M2；i++) for(j＝0；j＜M2；j++)　　Mi[i][j]＝i＝＝j；}invert2M( ){int i，j； initMi( )； for(i＝0；i＜M2；i++) {divrow(i，Mx[i][i])；for(j＝i+1；j＜M2；j++) if(Mx[j][i]) {divrow(j，Mx[j][i])；　　decrow(j，i)； } } for(i＝M2-1；i；i--)for(j＝i-1；j＞＝0；j--) submultdrow(j，Mx[j][i]，i)；}divrow(i，denom)int i；float denom；{int k； for(k＝0；k＜M2；k++) {Mx[i][k]/＝denom；Mi[i][k]/＝denom；}}<!-- SIPO <DP n="15"> --><dp n="d15"/>decrow(i，j)int i，j；{int k； for(k＝0；k＜M2；k++) [Mx[i][k]-＝Mx[j][k]；Mi[i][k]-＝Mi[j][k]；}}submultdrow(i，factor，j)int i，j；float factor；{int k； for(k＝0；k＜M2；k++) {Mx[i][k]-＝factor*Mx[j][k]；Mi[i][k]-＝factor*Mi[j][k]；}}/* returns the difference in microseconds between two timeval structures. */diff(tp1，tp0)struct timeval tp0，tp1；{return(tp1.tv_sec-tp0，tv_sec)*1000000+tp1.tv_usec-tp0，tv_usec；}/* this makes an ascii file which can be cut and pasted into spreadsheet */write out data_for_msexcel(filename)char *filename；{int i； FILE *fp； fp＝fopen(filename，″w″)； if(fp) {for(i＝0；i＜M；i++) fprintf(fp，″％4.3f\t％4.3f\t％4.3f\t％4.3f％c\n″，　　 20*log10(xmags[i])，　　 20*log10(xwmags[i])，　　 20*log10(xsmags[i])，　　 20*log10(xpmags[i])，　　 13)； fclose(fp)；}else{fprintf(stderr，″Could not open output file.\n″)； perror(filename)；}}