CN1790916B

CN1790916B - 使用滤波器乘积的线性校正器的校准***

Info

Publication number: CN1790916B
Application number: CN 200510137332
Authority: CN
Inventors: K·R·斯拉文
Original assignee: Tektronix Inc
Current assignee: Tektronix Inc
Priority date: 2004-11-04
Filing date: 2005-11-04
Publication date: 2012-05-02
Anticipated expiration: 2025-11-04
Also published as: CN1790916A; CN1815893B; CN1815893A

Abstract

公开了一种用于使用滤波器乘积的和来校准线性校正器的校准***，以及一种校准该线性校正器的方法。该校准***包括用于将测试信号引入例如ADC的信号处理***的第一和第二信号发生器。提供采集存储器和处理器用于采集和分析该信号处理***的输出，并且然后编程该滤波器系数到该线性校正器中。该校准方法分析从该信号处理***采集的互调和谐波分量，然后为该滤波器寻找幅度和相位响应。然后使用该幅度和相位响应来确定一组滤波系数。

Description

使用滤波器乘积的线性校正器的校准***

交叉引用

本申请要求2004年11月4日提交的申请号为60/625372的美国临时申请的优先权。

依照37 C.E.R.1.71(e)，申请人表示，本公开的一部分包含属于和要求版权保护的内容，例如但不限于源码列表、屏幕截图(screen shots)、用户界面、或者用户指令，或者提交的任何其他方面，对这些方面在任何权限中版权保护是或者可能是有用的。本版权所有人对于任何一个专利文件或专利公开文本的传真复制没有异议，例如出现在专利商标局专利文件或记录中。所有其他权利被保留，并且基于本申请或其任何部分的内容、公开显示和公开执行的派生作品的所有其它复制、分发和创建都被适用的版权法所禁止。

技术领域

本发明涉及线性误差校正，尤其涉及一种使用滤波器乘积(filter products)减少或消除由信号处理***例如模数转换器(ADC)产生的失真的线性校正器。

背景技术

减少ADC所产生的失真增加无寄生动态范围(spurious-free dynamic range)(SFDR)，其对于使用ADC来获得数据的***是有用的，所述***是例如谱分析器和其他电子测量仪器。

现代高速ADC设计使用深时钟流水线(deep clock pipeline)来帮助将模拟输入通过一系列改进的步骤精确地转换为采样的数字表示。ADC设计者付出了很大努力来去除模拟处理电路中的明显的非线性源。然而，一般很难去除所有的误差源。设计者试图去除该电路中的最明显的问题直到计算机化建模例如SPICE建模显示该转换器符合技术要求。可以通过使用例如减少非线性装置的动态范围的技术或者使用其周围的反馈来改善线性。然而，一些电路布局具有不能被完全去除的内在失真机制。

流水线处理还为内部数字和模拟电路活动提供对内部模拟信号处理进行调制的机会。在许多这种情况中，具有自身或其自身衍生物的线性函数的输入信号的自调整产生了残留的非线性失真。这导致了一些难以消除的低水平失真。这种调制可以通过内部电源分配发生。在这种情况下，能够在该电源轨上产生电压调制的电路通道数量可以是相当高。模拟这种效应使得装置建模复杂化并且减缓了计算机模拟。对第一阶来说，对于电源调制的这些影响将会几乎线性增加，所以它们可以作为线性有限脉冲响应(FIR)滤波器建模。

在模拟信号处理中的一个或多个点发生调制，其对应于乘法。在流水线ADC中，调制一般发生在转换级之间的高增益模拟放大器中。在这种情况中，典型地该谐波和互调失真的特征在于，存在第2阶和第3阶失真项，其中发生非常少的高阶失真。

以前提出的解决方案是基于Volterra滤波器的。ADC的脉冲响应可以是许多个时钟周期，例如可以使用64个时钟周期。使用Volterra滤波器的校正***将需要相似的响应长度。在第3阶失真Volterra***中，这导致大约（N³)/6抽头的滤波器，其对于具有64响应长度的校正***将导致大约50000抽头阶。具有这样巨大数量抽头的滤波器***太复杂和太昂贵以致于当时无法在实际***中实现。

在其他与扬声器中的校正失真相关的应用中提出了另一种解决方案，其使用了与Volterra滤波器的某些方面近似的滤波器结构。图1示出了具有第1阶校正和第3阶校正的这种方案的一种型式。第一阶补偿是由滤波器12(h₁)来提供的。第三阶补偿是这样提供的，即通过使用乘法器18将滤波器14的输出与滤波器16的输出相乘，使用滤波器20对乘法器18的输出进行过滤，使用乘法器24将滤波器20的输出与滤波器22的输出相乘，最后使用滤波器26对乘法器24的输出进行过滤。通过使用加法器28将来自第一阶补偿的输出和第三阶补偿相加，可以提供线性三次补偿(cubic compensation)。图1中所示***的第3阶补偿执行下述方程式，

其被描述为一般的第3阶非线性滤波器结构。这种实现方案在乘法器18之后使用了滤波器20，以及在乘法器24之后的滤波器26。一旦获得该线性三次补偿，就从正被补偿的该未知***的输出中减去。这需要该校正器能眵访问输入到该未知***的原始信号，而当该原始信号不是数字的时候是不可用的。虽然它可以是Volterra滤波器的有用的子方案，但是它也有不适用于具有良好线性频率响应的***的缺点。紧接在该乘法器之后的滤波器不能区分原始分量和由在先乘法的非线性效应所引起的混叠分量。虽然该乘法器之后的其他滤波可以提供对于信号通道中的频率相关的幅度和相位响应的一些校正，但是当用于使用大部分Nyquist频带的应用中时，混叠不允许该滤波器对该原始分量和混叠分量之间的相位和幅度响应差进行校正。

线性补偿***中的遗留问题涉及校准。这些***会需要求解相对于该输出成非线性的滤波器系数的***。对于任何可以应用到该***的校准方案，求解更多的系数将需要更多的计算量。

在现有解决方案上的细节和改进将在以下更加详细地讨论。

发明内容

如果可以通过恢复ADC中的等效失真滤波器的系数来对该失真机制建模，该ADC输出就可以通过以基本上相同的方式使该信号失真的数字处理网络，然后，减去该失真以减少或消除该ADC失真。虽然将所有ADC失真完全消除是不可能的，但是本方法改善了ADC的无寄生动态范围(SFDR)。例如，依赖于该ADC的特性，具有SFDR为80dB的ADC可以改进15dB的因子。本改进还从先前提出的布局中去掉了一些滤波器，从而简化了具有相对平滑线性频率响应的***使用的设计。这种简化可以对校正***中的更长滤波器交替使用，从而获得具有相同处理量的改进的性能。该改进对于精确测量应用是有意义的，例如那些与频谱分析器、示波器或者其他使用ADC的测量仪器有关的精确测量应用。

为了实现这些优点，提供一种合适的校准***和方法。该校准***包括：产生第一模拟信号的第一模拟信号产生器，产生第二模拟信号的第二模拟信号产生器，以及连接到该第一模拟信号产生器和第二模拟信号产生器的、用于将该第一模拟信号和第二模拟信号相加并提供模拟输出的模拟加法器。该模拟输出被输入到ADC以获得数字输出。将执行对滤波器乘积求和的线性校正器连接到该ADC的数字输出。在正常操作期间，该线性校正器将提供与该ADC的输出相比具有减小的失真的校正输出。将采集存储器连接到该数字输出；并且将用于控制该第一信号产生器和第二信号产生器的处理器连接以读取该采集存储器数据。还将该处理器连接到该线性校正器以便为该滤波器乘积编程滤波器系数。

还提供了一种校准线性校正器的方法。该校准方法通过将第一波形和第二波形输入到具有输出的信号处理***来完成。当其中每一个从该信号处理***中输出时，该采集存储器和处理器采集和分析该第一和第二波形的互调和谐波分量。通过基于失真模型以及该互调和谐波分量来确定目标函数的最小值，寻找连接到该装置输出的线性校正器中使用的滤波器的幅度和相位响应。然后，使用该幅度和相位响应来计算一组滤波器系数。

附图说明

图1(现有技术)是用于补偿扬声器的现有技术的线性校正器的布置的框图。

图2是包括第1阶、第2阶和第3阶失真的补偿的、基于滤波器乘积的线性校正器的框图。

图3是包括第1阶和第3阶失真的补偿的、基于滤波器乘积的线性校正器的框图。

图4是包括第1阶、第3阶和第4阶失真的补偿的、基于滤波器乘积的线性校正器的框图。

图5是补偿ADC的框图。

图6是包括补偿和校准***的ADC***的框图。

图7示出了基本校准过程。

图8示出了该校准过程的附加子步骤。

图9示出了该校准过程的附加子步骤。

图10示出了该校准过程的附加子步骤。

具体实施方式

本说明书包括第21-51页的用于说明本发明的具体实施例的具体实施方式的示例源代码及其注释文件，文件名称是Slavin.txt。

如上所述，在前提出的解决方案基于Volterra滤波器。然而，由于Volterra滤波器非常大并且很难与ADC一起实现，因此，需要一种解决方案，其利用更易管理的滤波器设计，同时还能够减少一些遗留的主要失真。以Volterra滤波器作为起始点，广义的非线性滤波器***能够在数学上定义为：

(方程1)

其中，N是该滤波器的脉冲响应长度，k是该滤波器阶指数(order index)。

例如，如果n＝3，则我们有DC值(h₀)、k＝1时的线性FIR滤波器项、k＝2时的2阶失真滤波器以及k＝3时的3阶滤波器的和。因而，对于n＝3，该Volterra滤波器可以表示为：

(方程2)

Volterra滤波器系数与输出y是成线性的，所以理论上，可以利用训练数据来找到一组h。一些乘积项刚好是相同组输入采样的置换(permutations)，所以对于每个阶指数k在该组h中的区别值的数量由下式给出，该数量对应于滤波器抽头的数量：

(方程3)

遗憾的是，流水线ADC***的该脉冲响应会是相当大的，所以N可以很大，对于k＝3产生很大数量的抽头。例如，如果该流水线ADC***的脉冲响应是64个时钟周期，从而N＝64，那么所需要的抽头数将是大约44000。附加的抽头被需要用于其他阶滤波器，如果存在的话。

现有的ADC线性校正器的实施例依赖于Volterra滤波器***的子集。该Volterra滤波器***的子集的特征在于：

(方程4)

对于1≤k≤n，***阶n定义了一组乘积阶。

尽管抽头的数量和值是未知的，但假定失真模型是这种形式。校正模型具有相同的形式，除非基于具有特定ADC构造的实验，事先选择阶和滤波器的长度。校准则包含找出滤波器抽头。注意到，通常滤波器抽头对于每个滤波器是不同的。对于***阶n＝3，忽略h₀(DC)项，我们得到：

(方程5)

这种结构的实施例其特征可以在于，具有使用滤波器实现的每个线性卷积的N抽头滤波器的乘积。

图2示出了用于实现方程式5的线性校正器100的实施例。信号处理***例如ADC的输出作为输入被提供给该线性校正器100。该线性卷积中的每一个，如方程5中所提供的，使用滤波器102到112来实现。该滤波器可以作为FIR滤波器来实现。第一阶项对应于滤波器102。在可选实施例中，通过用近似等于其他滤波器长度的一半的固定延迟代替滤波器102来获得该第一阶项。在另一实施例中，通过用固定延迟和滤波器的组合代替滤波器102来获得该第一阶项，从而使该总延迟近似为其他阶滤波器长度的一半。通过使用乘法器120将滤波器104和滤波器106的输出相乘以产生二阶滤波器乘积来实现第二阶项。通过使用乘法器122将滤波器108、滤波器110和滤波器112的输出相乘以产生三阶滤波器乘积来实现第三阶。然后，使用加法器124将来自滤波器102的输出与乘法器120和乘法器122的输出相加以提供该滤波器乘积的简单和(simple sum)作为输出。这里使用的术语“简单和”表示将该乘法器的值相加而不需在该乘法器和加法器124之间进行额外滤波的操作。这种简单和通过将该乘法器直接连接到加法器来获得。这里使用的术语直接连接(或者直接被连接)的含义是，在该通道上没有滤波器或其它处理元件，在该通道中可以有不改变通道上信号数据的本质特性的寄存器或其他元件。该输出现在是具有由该信号处理***例如ADC产生的减少了的非线性的补偿信号。应当注意的是，本发明的实施例消除了现有技术的解决方案中提供的、乘法器之后的滤波器。虽然与图1所示现有技术相比这将需要使用具有额外抽头的滤波器，但是它使得能够使用该滤波器补偿某些失真，而不需要来自该乘法器之后的滤波器的额外补偿。例如，如果图1的现有技术使用了一半时钟周期延迟的全通输出滤波器(所谓sm(x)/x或sinc(x)滤波器)，则该滤波器可以被结合到该乘法之前的滤波器中。通过仅在乘法器之前使用滤波器，可能的是该滤波器乘积***就能够更好地区分原始分量和混叠分量。

虽然在一些实施例中可以通过使用不同长度的滤波器来减少计算量，但是更长长度滤波器的使用会增加在校准期间需要求解的变量的数量，这将会使该校准算法变慢。对于硬件实现，更长长度的滤波器也会需要额外的延迟以匹配该滤波信号延迟。因此，在该校正器100的实施例中，所有滤波器长度都相等。

图3示出了设计成用于补偿第1阶和第3阶失真而不补偿第2阶失真的线性校正器100的实施例。在一些应用中，第二阶失真显著不到足以证明包括第2阶补偿是有道理的。如图3所示，通过使用乘法器122将滤波器108、滤波器110和滤波器112的输出相乘以产生第三阶滤波器乘积来提供第三阶补偿。然后该第3阶滤波器乘积和第1阶滤波器乘积的简单和可以提供具有减少或者消除了的第1阶和第3阶失真的补偿信号。

如图4所示，能够提供包括第4阶补偿的校正器100的实施例。利用方程4的一般形式，对于***阶n＝4，并且忽略该h₀(DC)项，我们就有：

(方程6)

如图4所示，该第4阶项可以通过使用乘法器148将滤波器140、滤波器142、滤波器144和滤波器146一起相乘来执行。此外，使用加法器124将该滤波器乘积与其他阶补偿直接求和，而在该乘法器148和加法器124之间没有任何中间滤波器。如从前面例子中可以清楚的，一个本领域普通技术人员将能够求解方程4，以便针对任何预期阶而将其以与方程5和方程6相似的形式放置，其然后可以如这里一般教导地使用滤波器乘积的简单和来实现。如图4所示，第2阶补偿被去除以考虑如下所详细讨论的校准。这种配置将在第4阶具有较大重要性的那些环境中提供补偿。类似地，能够提供对于第4阶和第5阶的补偿并且消除第2阶和第3阶以考虑如下所述的校准。

对于所提出的通常的Volterra形式的分解的有效性以及相应的滤波器乘积结构的证明在于理解对于每个乘积阶k存在单一的自调制机制，而不在于任何随机Volterra滤波器***能够以这种方式被分解的可能性上。

线性校正器100的各种实施例可以利用专用的硬件执行，例如FPGA或者ASIC，或者利用运行软件的通用处理器执行。当前，尽管运行在通用处理器上的软件对于后采集校正是有用的，但是FPGA或者ASIC对于执行实时校正是有用的。在未来，也有可能为了实时校正利用在通用处理器上的软件。

虽然在一些实施例中，该线性校正器100用于补偿从ADC输出的信号，而在其他实施例中，线性校正器100的结构可以集成在与该ADC相同的封装(packing)中，或者可能集成在与该ADC相同的芯片上，以便形成补偿ADC。该补偿ADC190在图5中示出。它包括ADC模块192，该ADC模块192包含将模拟信号转换成数字信号的各种电路。该ADC模块192的数字输出被输入到可以上述教导实现的线性校正器100。该线性校正器的输出是具有减少了失真的补偿输出。这种组合结构提供了补偿ADC。

为了适当地优化上述线性校正器，有必要对该线性校正器进行校准以对每个滤波器确定合适的滤波器系数。与一般Volterra滤波器不同，图2-4所示的校正器的滤波器乘积输出不是关于其系数成线性的，所以在一般情况下，该滤波器系数的计算是一个非线性优化问题。

因此，提出了一种具有校准的线性校正器***和一种校准方法。如图6所示，配备有线性校正器100和校准电路的ADC 200。提供模拟选择开关202以控制进入该ADC的输入。该模拟选择开关202由处理器204控制。在校准期间，处理器204控制该模拟选择开关202以允许校准频率信号进入该ADC 200的输入。处理器204还控制数字选择开关206，以控制在校准期间直接来自该ADC 200的数字输出是否进入采集存储器208。处理器204还控制第一校准频率发生器210和第二校准频率发生器212。在校准期间，该第一校准频率发生器210和第二校准频率发生器212的输出在在输入到ADC 200之前，在求和电路214中组合。在本校准方法的实施例中，进入ADC的该模拟校准信号和模拟求和与该被校准的ADC相比具有显著少的失真。

在校准期间，该采集存储器208采集未校正的ADC输出用于校准。处理在校准频率对的范围上多次采集的数据。该采集的数据由处理器204进行处理。

在正常操作期间，通过模拟选择开关202将该模拟输入输送到ADC 200，并使用该线性校正器100对该ADC的输出进行校正。该校正器的基本操作已经结合图2-4在上面进行了描述。具有该线性校正器的滤波器可以通过处理器204加载。然后可以使用第二处理器216提供其他处理，例如可以就像对该ADC的直接输出那样处理该校正的输出。图6所示的实施例示出了用于校准和其他处理的采集存储器208。在可选实施例中，可以使用单独的存储器用于其他处理。类似地，图6所示的实施例示出了第二处理器216，在其他实施例中，处理器204能够提供这种额外的处理和控制该校准过程。

如图7所示，校准线性校正器中的滤波器乘积内的一组滤波器，有四个基本步骤。对于第1阶以上的每个乘积阶，使用该乘积阶作为参数执行该步骤，例如第2阶或第3阶。步骤310采集和分析输入频率对的互调（IM)和谐波分量。步骤340使用步骤310的结果寻找该校正网络中的滤波器幅度和相位响应。步骤360使用在步骤340中找到的幅度和相位响应，根据所需的幅度和相位响应使用非线性相位滤波器设计算法来设计每个滤波器。步骤380通过设定该滤波器抽头来编程该线性校正器100内的滤波器以对ADC失真进行校正。

步骤310的基本信号采集和分析可以通过选择如图8中的步骤410所提供的基数据记录长度来完成。在步骤420中选择一组校准频率。步骤430生成两个基于额定校准频率的正弦输入，以及在步骤420中选择的校准指数。步骤440为每个校准中心频率指数计算一组谐波和互调频率。步骤450测试可能引起频率分量重叠的混叠。步骤460或者i)恰在所计算的谐波和互调频率执行一组离散傅里叶变换(DFT)，或者ii)执行为完整的一组可能频率生成结果的快速傅里叶变换(FFT)，然后，选择在该计算的谐波和互调频率上的输出结果。

在步骤410，选择用于DFT分析的基数据记录长度L。为了提高效率，在该校准方法的某些实施例中，将L选择为2的幂。对L的选择对应于基本分析频率F＝f_s/L，其中f_s是采样率。然后将所有频率f_m表示为基本分析频率F的谐波(整数倍)m，即f_m＝mF。角频率被表示为：

ω_m＝2πf_m＝2πmF＝2πmf_s/L (方程7)

在步骤420，选择将用于随后步骤以设计滤波器响应的一组校准频率。这些频率被选择为F的谐波倍数，并且规则间隔在输入频率可以发生的整个尼奎斯特(Nyquist)频带中，但总是排除0(DC)的极限(extreme)和f_s/2。

在步骤430，对于具有额定校准频率ω_c的每个校准指数c，以角频率ω_c-1＝ω_c-ω₁和ω_c+1＝ω_c+ω₁生成两个正弦输入，其对应于：

x_c(t)＝p_c-1sin(ω_c-1(t+λ_c-1))+Q_c+1sin(ω_c+1(t+λ_c+1)) (方程8)

其中λ_c是在频率指数c的延迟。在该校准方法的实施例中，幅度｛P_c-1,Q_c+1}是相似的以改善IM校准结果。任何微小差别可以在之后的校正阶段中去除。在该校准方法的实施例中，选择该校准频率以使得它们使用该ADC输入动态范围的大部分。然后，采集整数倍M组L个相邻的A/D输出样本。每个第M个样本被平均以获得低噪声y(t)，t＝0...L-1，其中L是上面步骤410中提供的该基数据记录长度。这种类型的简单平均是非常快的，允许对ML的大量值(可以是百万样本或更多)进行快速处理。L应当被选择成使得它不会太小或者使分析频率会分隔的太远。大的L的值可以生成比所需要的更多的校准频率数据，导致后面的算法运行更慢。L＝2¹⁰＝1024的值已被成功地使用。

在步骤440，对于每个校准中心频率指数c(其中0<c<L)，产生一组谐波和互调频率。DC项被忽略，因为DC误差会由于许多其他原因而发生。例如，第二阶***在{ω_c-1-ω_c+1，2ω_c-1，ω_c-1+ω_c+1，2ω_c+1，}具有IM和谐波失真频率，对应的指数：m＝{2，2c-2，2c，2c+2}。事实上，能够只使用和频(sum-frequency)项m＝{2c-2，2c，2c+2}来校准。仅使用频率的和项而不使用差项的优点在于，前者不是由任何低阶失真滤波器乘积生成。这包括可以称为第1阶(线性)“失真”***的原始输入信号。

应当注意，在采样的***中，IM和谐波失真会包含混叠频率分量。因此，在步骤450，执行测试以处理混叠导致频率分量重叠的情况，因为该测量将与该失真模型不一致。对于给定c，生成一组可能的谐波和IM项(DC除外)。在一些实施例中，该组将包含除DC项之外的所有可能的谐波和IM项。对于来自步骤440的列表m中的指数m_i的每个频率，由于采样理论不能区分频率间隔为L的混叠分量，所以它们相互卷绕如下：

w(m_i)＝m_imodL (方程9)

而且，在该卷绕间隔L内，混叠不能在w(m_i)和L-w(m_i)之间区分。该混叠的指数重新映射(aliased index remapping)是较低范围的解释(interpretation)，由下式给出：

alias(m_i)＝min(w(m_i)，L-w(m_i)) (方程10)

其中min()是其输入值的最小值。因此0≤alias(m)<L/2。m_i到alias(m_i)的这种重新映射仅用于对于给定整数m_i和L检测任何失真频率冲突。冲突检测大多很容易通过生成所有谐波和IM频率指数来实现，然后如上面在方程9和方程 10中所述重新映射每个频率指数，并且对得到的混叠频率指数值进行数字排序。然后，会发现具有相同重新映射频率指数值的相邻频率指数。如果发现有冲突，在该c值处不进行测量。如果L是2的幂，并且c是偶数，那么{c-1，c+1}总是互素的。从而，冲突仅发生在混叠频率。即使这样，通常冲突也是很少的。例如在第3阶***中，冲突仅发生在c＝L/4时。在这种情况下，互调分量位于2(L/4-1)+(L/4+1)＝3L/4-1，所以其混叠频率是在L-(3L/4-1)＝L/4+1，该频率也是输入频率。类似地，(L/4-1)+(2L/4+1)＝3L/4+1，其混叠到上/4-1，也是输入频率。

在本校准方法的一个实施例中，在所选择的谐波和互调频率处执行一组L-样本的实数输入DFT，如步骤460所示。在每个频率m处的实数和虚数分量使用对于该平均实数样本序列的单个频率DFT来得到：

X_m＝Re(Z_m)

Y_m＝Im(Z_m) (方程11)

其中Re和Im是复输入的实部和虚部。方程11的形式是对于来自结合上述步骤410提供的方程7的正规化采样频率。在可选实施例中，可以使用快速傅里叶变换(FFT)代替多重DFT来获得在每个频率m处的实和虚分量。

在该失真***中，每个滤波器在这个阶段(stage)具有未知的幅度和相位响应。为了找到这些响应，提供用于该***的失真的数学模型。基于提供的该模型，可以执行非线性优化。用于第N阶乘积的自调制模型将n个未知滤波器响应表示为每个乘积项的修改输入：

(方程12)

这种自调制模型可以以可选的一般形式重写为：

(方程13)

其中：

(方程14)

并且

B_i，c+1＝A_i，c+1/A_i，c-1 (方程15)

假定有足够大的L，如果用于校准的两对频率C-1和C+1足够接近，则对于给定滤波器，在每个频率的延迟可以假定为近似相同，那么方程15可以简化为：

（方程16)

对于该第n阶失真***，变量数等于2n+1。对于使用的每一对校准频率，到ADC的正弦波输入幅度都严密匹配。否则，该***不能区分在两个不同频率处的它们的幅度之间的差别和每个滤波器的幅度响应中的差别。

如果在该两个不同的校准频率处通过每个滤波器的幅度响应是近似相同的，那么更简化的自调制模型可以提供如下：

(方程17)

通过选择适合这些失真模型中的一个的ADC，如方程13、方程16或方程17所表示的，对于每个滤波通道在每个频率可以找到一组相位和幅度解。表示该失真模型的任何方程都可以被扩展以提供一组互调和谐波频率分量，其对应于上面在步骤440所找到的那些频率分量。该减少可以表示为在每个频率的余弦(实部)和正弦(虚部)项。一旦获得频率项的扩展和收集，连同在步骤460确定的那组DFT一起，那么就可以以符号表示的形式(symbolically)得到在这些表达式的实部和虚部与DFT结果之间的差的平方和，其可以从如上所述的DFT或FFT得到。这导致了对于第2阶失真的大的方程，以及对于每个更高阶的甚至更大的方程。

因为相位和延迟λ通过角频率ω相关的，即通过：

(方程18)

代入来自方程7的ω_m，并且使用正规化频率(即f_s＝1)，使得模型方程式中的任何方程13、方程16或方程17以频率指数的形式重写。例如，方程16可以重写为：

(方程19)

例如以阶数n＝2分解成实部和虚部，我们扩展方程19并且收集项以获得在四个非DC互调和谐波频率指数的项：{2，2c，2c-2，2c+2}：

(方程20)

(方程21)

(方程22)

(方程23)

其中该组系数{C，S}被确定为：

对于方程20：

对于方程21：

对于方程22：

对于方程23：

现在，对这些系数和它们相对应的来自方程11的DFT或FFT分量之间的差的平方和进行估计：

(方程24)

这里方程24可以表示该目标函数或者误差。对于第2阶的情况，要寻找的未知量可以用向量表示为：

x_c＝{K_c，B_1，c+1，B_2，c+1,λ_1，c,λ_2，c} (方程25)

一旦已经提供了该目标函数，就可以确定每个滤波器的幅度和相位响应。现在参照图9，示出了本校准方法的一个实施例。如步骤540中所示，首先将该目标函数最小化。在最小化该目标函数之后，对幅度执行解后校正(post-solution correction)，如步骤560所示。该校正解现在被分解到为每一阶提供的滤波器之间，如步骤580所示。

为了最小化该目标函数，如步骤540所示，将需要求解非线性方程的***。可以使用一类基于牛顿(Newton-based)的算法来对每个预期的阶求解这些非线性方程和最小化该目标函数。可能的基于牛顿的算法包括修正牛顿(MN)、高斯一牛顿(GN)和准牛顿(Quasi-Newton)(QN)。在每个搜索迭代k，这些方法都估计该目标函数的数值梯度向量(g_k)和该目标函数的Hessian(S_k)的数值逆，或者对于GN和QN方法是该目标函数的Hessian函数的逆的近似。该目标函数的梯度向量和该目标函数的Hessian的逆可以描述如下：

在x＝x_k (方程26)

在x＝x_k (方程27)

通过估计该梯度向量和逆Hessian，能够获得可用于最小化该目标函数的搜索方向(d_k)，其中该搜索方向如下给出：

d_k＝-S_kg_k (方程28)

该目标函数

的梯度向量定义为这样一个向量，其第i个分量是该函数关于其第i个未知变量的偏导数。因此，它是这个阶段的符号结果而非数值。对于对应于方程16的具有5个未知量的第二阶情况，该目标函数由方程24提供：

Hessian矩阵的第(i，j)个元素是该目标函数对于其第i个和第j个变量的偏导数。再次，对于第2阶情况，该Hessian矩阵如下给出：

因此该矩阵元素是该未知变量x的项。该偏导运算符是可交换的，所以Hessian矩阵总是平方对称的(square-symmetric)。

对于该目标函数可以找到相应的梯度和Hessian矩阵。虽然求解这种复杂表达式的导数是非常繁复的，但是可以通过使用符号数学软件工具例如在Mathematica中的符号导数函数来执行计算机化符号计算。可选地，可以使用众所周知的近似来数值地求解任何函数f(x)的导数f′(x)：

(方程29)

对于适当小的值Δ，假定该计算精度可用。可以使用符号或者数值微分。符号微分可以通过符号数学软件来提供，通过这种软件提供的这种处理对处理器204使用通用处理器。应当注意，如果已经得到了符号微分，在本校准方法的一些实施例中，可以使用较低级语言例如C语言将这些符号结果编码为函数，以便改善在处理器204中的这种操作的性能。在可选实施例中，可以由处理器使用数值微分代替提供该符号解的编码实施。

如上所述，可以通过使用基于牛顿的算法直接或叠代求解逆Hessian。例如，在本校准方法的一个实施例中，使用准牛顿算法叠代确定该Hessian矩阵的逆的近似值。有多个QN叠代的子范畴可以使用。例如可以使用BFGS(Broyden，Fletcher，Goldfarb，Shanno)叠代公式。在本校准方法的一个实施例中，该BFGS叠代公式与弗莱彻线搜索(Fletcher’s Line Search)算法一起使用。该逆Hessian S_k+1的第(k+1)次BFGS叠代如下给出：

(方程30)

其中列向量p_k、标量q_k和方矩阵r_k如下给出：

p_k＝S_kγ_k

(方程31)

S_k(方矩阵)是该逆Hessian的第k次近似值，而S₀＝I为单位矩阵。应当注意，上标T是转置，例如：是δ_k的转置。并且列向量：

γ_k=g_k+1-g_k (方程32)

其中g_k是较早从(方程26)中获得的梯度向量。对于下一个叠代的该新的解向量估计x_k+1由应用到该前一个x_k的校正列向量δ_k给出：

x_k+1＝x_k+δ_k (方程33)

为了寻找每个叠代的未知的δ_k，使用沿着搜索方向d_k的约束搜索。一旦使用方程33找到x_k并且使用方程28和方程30找到d_k，然后，就找到了非负标量α_k，其最小化：

e₂(x_k+1)＝e₂(x_k+δ_k)＝e₂(x_k+α_kd_k) (方程34)

这是α_k中的一维优化问题，而不论该原始问题的维数如何。因而，对这种特定类型方程的求解已知为线搜索。

当前的解可估计向量x_k，且其相应的搜索方向d_k被给定作为弗莱彻线搜索算法的输入，以给出标量值α_k。于是，对于方程33中的对x_k的列向量校正如下给出：

δ_k＝x_k+1-x_k＝α_kd_k (方程35)

众所周知，在一般情况下，寻找目标函数的全局最小值是一个困难问题。所找到的最小值总是依赖于该目标函数，但是也可以依赖于该初始条件x₀°如果该失真模型能够相当好地符合，那么通常收敛会迅速和可靠。然而，如果该真实的失真机制包括多个失真乘积或者其他一些不遵循该失真模型的机制，则该解会发散。在这种情况下，向一个解前进，并且该误差函数方程24在每次叠代中减少，但是向量x中的未知量之一按指数规律向无穷发散。对BFGS叠代数量的限制可以可靠地检测这种问题。然而，到达这种限制的过程也减慢了该校准算法。

如果在噪声减少平均后，看到发散，那么可以设置发散失败标记以指示对于进一步诊断的需要，或者该ADC的失败。

然而，一旦逼进最小解，则算法收敛是二次的，所以能迅速收敛到可重复的结果。这种特性可以用于探索方程16的模型的解空间。如果一个解收敛到在先找到的解，则可以作为重复解而拒绝。该过程会继续，直到已经获得不能再找到最小值的统计上符合要求的判断。

在本校准方法的一些实施例中，在具有π的邻近整数倍的角度处可能有其他邻近解。这些重复解由于微不足道地相关或“无兴趣”而被辨别并被排除。

找到解的一个方法是在方程19中以B的随机值开始x，使用：

B＝e^{u(Random-1/2)} (方程36)

Random是从0到1的随机实数值。u＝0.1的值对于实数ADC数据可以有很好的效果。每个角

可以在以下范围中随机产生：

(方程37)

然后可以使用方程18在已知频率将角转变为方程19中的延迟。

在另一实施例中，通过将用作校准指数c的函数的随机角的初始范围限制到在先前校准中的起作用的范围来寻找已知ADC设计的解。对初始范围的这种限制会加速该解的搜索。

即使该开始值至少被方程37约束，该最终值也可以存在于这些范围之外而进入最接近的相邻重复解中。这是因为该开始点可以存在于在该重复最小值之间的重复的“峰”的另一侧。然而，由于在该间隔π的重复的函数“形状(landscapes)”的相似性，即使更远也不会找到重复解。在本校准方法的一个实施例中，找到了最近的相邻解并且排除了重复解。

对于每个校准频率，一旦找到第一解angle，就可以很容易地找到中心解(即角最接近0)。在本校准方法的一个实施例中，所找到的角直接被验证为中心解(在±π/2内)。在可选实施例中，使用multiple＝floor(angle/π+0.5)来寻找距该中心解区域最接近的π的整数倍，然后将该π的倍数从该第一解角(solutionangle)中减去以使该角进入该中心区域中。在一些应用中，接近±π/2的解可以随着在另外具有校准频率的角度平滑变化中的在+π/2或-π/2之间的噪声波动。这些波动可以通过在之后的处理中查找角的连续性而得到校正。然后，该BFGS算法的第二次运行可以通过使用被替换的第一解作为好的开始点而迅速找到更准确的中心解。

从该中心解，可以***地找到距离该中心解最近的相邻解。如上面寻找该中心解的过程，可以通过将{-π，0，π}的所有组合加到每个角并且使用作为结果的向量作为BFGS解搜索的开始点来找到该重复解。再次，因为每个开始点非常接近该正确值，BFGS算法非常迅速地收敛。对于第2阶搜索有3²-1＝8个相邻解，对于第3阶搜索有3³-1＝26个相邻解。

找到这些解之后，可以使用在其他随机开始位置的搜索来寻找任何更多的解，如果它们存在的话。如果一个解与先前找到的解中的一个匹配，则拒绝它。如前面所述，通过由具有符号梯度的BFGS算法找到的解的高重复性而使得验证重复(duplicate)的过程变得更简单。

如果一个新解与迄今所找到的任何一个都不匹配，则可以找到其中心解，并且使用与前述相同的方法找到其所有相邻解。在本校准方法的一个实施例中，富有成效的搜索结果重设超时计数器(tmeout counter)，而寻找新解的失败则减少该超时计数器。如果该超时(timeout)到达零，则停止对其他解的搜索。选择根据方程24给出最好的误差测量的中心解作为对于该校准频率指数的最终解。

在本校准方法的一个实施例中，如果已经在先前的校准频率指数c_prev找到一个解，则可以使用它作为在c的下一次搜索的开始点。如果它不能收敛，则使用前述的慢的随机搜索方法。

如果对于特定的c没有找到解，可以使用在c_next的解作为开始点进行第二次扫描(pass)。如果仍然没有找到解，那么可以使用相邻值之间的修补(patching)或内插。如果对太多的校准值指数都没有能找到解，则该ADC是有故障的，可以设定一个故障标记。

如果对于所有校准频率已经找到可以接受的一组懈，则进入该自调制***的乘积的每个通道的幅度和相位响应的特征在于该组选定的离散的校准频率。利用这些频率的均匀分布，对于每个通道的响应可以由如图2-4中的FIR滤波器所表示。在本校准方法的一个实施例中，现在对每个滤波器的系数进行叠代寻找。

在本校准方法的另一实施例中，对来自方程25的该组向量x_c在将它们传送到该滤波器设计之前进行校正，如步骤560所示。在本校准方法的实施例中，该校准输入的幅度在合理的校准范围上将不会显著影响该滤波器响应。从方程14可以看出，这种依赖性通过用除K_c而可以被消除，其中n是该乘积阶指数。而且，在一些实施例中，在方程8的两个校准正弦输入的幅度中的任何差别还应当对该滤波器响应具有可以忽略的影响。在这种情况下，不仅该失真滤波器的相对响应，而且该校准正弦输入的相对幅度都会影响该测量结果。因此，在这些实施例中，使用相对幅度P_c-1和Q_c+1来进行校正。该比率：

R_c＝Q_c+1/P_c-1 (方程38)

给出了在x中B参数被高估的量。因此对方程25中的x的该参数针对第n阶校准进行校正以给出新向量X_c如：

(方程39)

在任何校准频率，滤波器乘积***的组合幅度响应是每个滤波器的幅度响应的函数。因此，我们可以使得一个滤波器的系数变大两倍，而另一个缩小两倍，而具有相同的结果。***幅度响应的变化在校准频率之间可以是重要的。我们可以通过使得乘积中的每个滤波器对于第n阶失真校正***具有相同的来自方程39中第一项的第1/n次方根响应，从而在滤波器乘积***中的所有滤波器中平均分散跟踪幅度响应变化的负荷，以使：

(方程40)

对于第n阶失真***，在频率指数c的每个向量X_c中有n对幅度和频率响应。对n个滤波器中的一个任意指定一对幅度和相位响应，对作为频率的函数的每个滤波器可能会产生非平滑幅度响应、相位响应，或者两者。非平滑响应可以产生不能对选择的滤波器长度很好地跟踪预期响应的滤波器设计。因而，该幅度响应和相位响应被分配到每一阶可用的滤波器之间，如步骤580所示。因为一组足够密集的校准频率通常会导致对于相邻校准指数的延迟响应之间的高的相关度，所以在本校准方法的一个实施例中，通过以该校准循环周期(1/f_c-1)的一半为模，测试相对于在频率指数c_prev处的前一组的新组可能的置换，来获得在该新组延迟和先前一组延迟之间的更好的配合。在奇数个滑移(slips)的情况下，将对应于半个循环的延迟应用到该滤波器之一。在本校准方法的实施例中，例如采用方程16所述的模型的实施例，每个相位值具有相应的幅度，所以重新排列相位也重新排列相应的幅度。如果两个相位交叉，它们可以被分组为对一个滤波器的最小相位值和对另一滤波器的最大相位值。在每个相位图上导致的绞缠(kink)不会像不连续性一样对所得到的滤波器质量不利，而是产生可行的滤波器设计，如果提供足够的滤波器抽头的话。

通过前面的校准步骤，对于每个滤波器路径都获得一组适当分组的幅度和相位响应。图6中的滤波器设计步骤360使用在步骤340计算的该幅度和相位响应设计每个滤波器，以给出对该计算的响应的近似。在射频(RF)应用的情况，使用该ADC的上Nyquist频带的中间部分，而不对该频带中的外部频率区域进行校准。这意味着该外部频率是不受约束的。如图10所示，在本校准方法的实施例中，使用叠代过程来完成该滤波器设计和寻找该滤波器抽头。在可选实施例中，可以使用简单的逆FFT。

如步骤610所示，对每个通道选择一个滤波器长度。在本校准方法的一个实施例中，对于所有通道使用相同的滤波器长度。在可选实施例中，对于部分或所有通道使用不同的滤波器长度。在本可选实施例中，可以需要额外的延迟以补偿该不同的滤波器长度。该滤波器长度可以至少部分依赖于该预期的性能和被校准的ADC。以示例的方式，使用了全部63抽头的滤波器长度或者全部127抽头的滤波器长度。

在步骤620，选择FFT长度并且将该幅度和相位响应转换为复直角坐标的滤波器频率响应向量。该FFT长度应当比在步骤610选择的滤波器长度要长。例如，在本校准方法的实施例中，选择512个样本的实数FFT长度。因为该FFT比滤波器长度长，所以这里将对应于滤波器抽头的频率之外的值初始化为零。

在步骤630，执行该滤波器频率响应向量的逆FFT以获得一组实数滤波器抽头，但是还有扩展到该FFT长度例如512的抽头。在时域中，超出该选择的滤波器长度的抽头被强制为零。

在步骤640，执行前向FFT以回到频域。该频率响应向量的一些区域应当是不受约束的。在所关心的频率的FFT槽(bins)被重写回对于该预期幅度和频率响应的原始值，使得该不受约束的频率槽(frequency bins)不变。

重复步骤630和640以叠代收敛到预期解。一旦在步骤630提供的该滤波器抽头已经从先前的叠代中变得稳定，就找到了该滤波器抽头。

一旦从步骤630获得该滤波器抽头，就可以如图3中的步骤380所提供的那样对硬件编程。该编程滤波器乘积***现在能够对ADC谐波和互调失真进行校正。

本领域普通技术人员将会清楚，可以对本发明的上述实施例细节进行许多改变，而不脱离其基本原则。因此，本发明的范围应当由所附的权利要求来确定。

Slavin.txt

Mathematica code for ADC Correction Algorithm.

(＊by Keith Slavin＊)

(＊COPYRIGHT NOTIFICATION

＊Pursuant to 37 C.F.R.1.71(e)，applicant notes that a portion of this disclosure

contains material that

＊is subject to and for which is claimed copyright protection，such as，but not

limited to，source code listings，

＊screen shots，user interfaces，or user instructions，or any other aspects of this

submission for which copyright

＊protection is or may be available in any jurisdiction.The copyright owner has no

objection to the facsimile

＊reproduction by anyone of the patent document or patent disclosure，as it appears

in the Patent and Trademark office

＊patent file or records.All other rights are reserved，and all other

reproduction，distribution，

＊creation of derivative works based on the contents，public display，and public

performance of the application

＊or any part thereof are prohibited by applicable copyright law.

＊)

BeginPackage[″correction`″，″Global`″]

(＊frequency Indexing:i1＝i-1，i2＝i+1＊)

(＊

＊design filters to correct ADC output.lowF，hiF are normalized

＊frequencies from 0 to 0.5(of fs).′data′is a 2-Darray of acquired

＊ADC data in/login/keith/adc_orrection/xcal_data6/summary.To load，

＊type:<<Z:adc_correction\xcal_data6\summary；The result is in″data″.

＊Type:ProcessData[data，0.1，0.4，127，10^-12]to get the filter

＊coefficients in the form:

＊{{2nd_orde_f1_list，2nd_order_f2_list}，

＊{3rd_order_f1_list，3rd_orde_f2_list，3rd_order_f3_list}}

＊)

ProcessData[caldata_，calFreqLo_，calFreqHi_，correctFilterLen_，error_]:＝

Module[

{order，responses2，correctionFilters2，responses3，correctionFilters3，

correctionFilters}，

order＝2；

Print[″Generate′rsponses2′:″]；

responses2＝ProcessCalData[ca]data，order，calFreqLo，calFreqHi，error]；

(＊design the correction filters from the cal responses＊)

print[″Generate′correctionFilters2＇:″]；

correctionFilters2＝DesignFilters[responses2，correctFilterLen]；

order＝3；

Print[″Generate′responses3′:″]；

responses3＝ProcessCalData[caldata，order，calFreqLo，calFreqHi，error]；

Print[″Generate′correctionFilters3′:″]；

correctionFilters3＝DesignFilters[responses3，correctFilterLen]；

Return[{correctionFilters2，correctionFilters3}]；

]

Page1

Slavin.txt

(＊

＊ Run a calibration on randomly generated distortion filters，and then correct

＊ the same distortion on a multi-sine input signal.

＊ (＊Note：0.2...0.3 range puts 3rd order IM in 0.1...0.4fs

＊ and 2nd order IM from 0...0.1 and from 0.4...0.5fs.＊)

＊ testFreqLo＝0.2；

＊ testFreqHi＝0.3；

＊ testFreqStep＝0.004； (＊step in fs units from testFreqLo to testFreqHi＊)

＊ calFreqLo＝0.1；(＊low calibration range in fs units＊)

＊ calFreqHi＝0.4； (＊high calibration rage in fs units＊)

＊ distortFilterLen＝63；(＊odd length of all distorting filters＊)

＊ spread＝4.0；(＊spread(samples)of random distorting filters＊)

＊ analysisLen＝1024；(＊samples in DFT analysis of cal frequencies＊)

＊ freqPairs＝128；(＊frequency pairs in cal(must divide analysisLen)＊)

＊ correctFilterLen＝63；(＊correction filter lengths(must be odd)＊)

＊ testLen＝32768；（＊samples used in test signal analysis＊)

＊

＊ these three parameters are in LSB units.They make

＊ noise+distortion+original tones with similar amplitudes to the

＊ AD6645 ADC:

＊

＊ distortionAmpl＝0.04；

＊ tonePairAmpl＝15000.0；

＊ noiseAmpl＝2.0；

＊)

Demo[]:＝

Module[

{}

MainTest[0.2，0.3，0.004，0.1，0.4，63，4.0，1024，128，63，32768，

0.04，15000.0，0.01]；

]

MainTest[testFreqLo_，testFreqHi_，testFreqStep_，calFreqLo_，calFreqHi_，

distortFilterLen_，spread_，analysisLen_，freqPairs_，

correctFilterLen_，testLen_，

distortionAmpl_，tonePairAmpl_,noiseAmpl_]:＝

Module[

{margin，testdata，order，

filter2，filter3,filters，

caldata，error，

responses2，responses3，

correctionFilters2，correctionFilters3，correctionFilters，

distorted，corrected}，

Print[″Version0.3″]；

(＊generate second-order random distortion filters＊)

print[″Generate＇filter2＇:″]；

order＝2:

filter2＝GenFilters[distortFilterLen,order，spread]；

(＊generate third-order random distortion filters＊)

Print[″Generate＇filter3＇:″］；

order＝3；

filter3＝GenFilters[distortFilterLen，order，spread]；

(＊combine filter data＊)

filters＝{filter2，filter3}；

(＊generate multi-tone test data record＊)

Page2

slavin.txt

Print[″Generate′testdata′:″]；

margin＝(distortFilterLen-1)+(correctFilterLen-1)；

testdata＝GenTest[testLen，mrgin，

testFreqLo，testFreqHi，testFreqStep]；

(＊generate a set of cal data tone pairs＊)

print[″Generate＇caldata＇:″]；

cal data＝DistortCal[filters，analysisLen，freqPairs，

distortionAmpl，tonePairAmpl，noiseAmpl]；

(＊process the cal data to get　filter responses＊)

error＝10^-12；

correctionFilters＝

ProcessData[caldata，calFreqLo，calFreqHi，correctFilterLen，error]；

(＊distort the test data using the original filters＊)

print[″Generate′distorted＇:″］；

distorted＝DistortData[testdata，filters，

distortionAmpl，tonePairAmpl，noiseAmpl]；

(＊plot the distorted signal＊)

PlotunwindowedFFT[Take[distorted，testLen]]；

Print[″Generate′corrected＇:″]；

corrected＝CorrectData[distorted，correctionFilters]；

(＊plot the corrected signal＊)

PlotUnwindowedFFT[Take[corrected，testLen]]；

]

(＊aliases data set　by first finding data mod n，and then

wrapping data>n/2 to n-data＊)

AliasF[data_，n_]:＝

Module[

{len，modData，i}，

len=Length[data]；

modData=Mod[data，n]；

Return[Table[Min[modData[[i]]，n-modData[[i]]，{i，1，len}]]；

]

(＊Expand distortion product＊)

AllTerms[i_，order_]:＝

Module[

{t，tt，j，q0，q1，q2，q3，q3a，q4，q5，q6，q7，q8，terms，len}，

If[order＝0，

Return[{0}]；

]；

q0＝TrigReduce[(Cos[(i-1)tt］+Cos[(i+1)tt])^order]；

(＊extract tt variable outside as a common factor for each cos term＊)

q1＝q0/.Cos[a_]:>Cos[Apart[a]]；

(＊expand into seperate terms＊)

q2＝Expand[q1]；

q3＝Table[q2[[j]]，{j，1，Length[q2]}]；

(＊find the cosine terms only，discard their coefficients，and convert

any DC term to 0＊)

q4＝Table[If[NumericQ[q3[[j]]]，0，q3[[j]]/(q3[[j]]/.tt->0)]，

{j，1，Length[q3]}]；

(＊eliminate the Cos[]arount each term＊)

q6＝q4/.Cos[a_]:>Apart[a]；

(＊form(k＊i+m)＊tt->k＊i+m for integer k，m＊)

q7＝q6/.a_＊tt:>a；

Page3

Slavin.txt

q8＝q7；

(＊make i term positive:if k is -ve，then its derivative is -ve＊)

terms=Table[If[D[q8[[j]]，i]<0，-q8[[j]]，q8[[j]]]，{j，1，Length[q8]}]；

Return[Expand[terms]]；

]

(＊

＊＊

＊S is an approximation to the inverse of the Hessian，G.

＊

＊BFGS iteration is:

＊ T T T T

＊ ga S ga dx dx dx ga S+S ga dx

＊ S[k+1]＝S+(1+--------)------ -(------------------)

＊ T T T

＊ dx ga dx ga dx ga

＊ T T T T T

＊ (dx ga+ga S ga)dx dx dx ga S+S ga dx

＊＝S+----------------------- -(------------------)

＊ T 2 T

＊

(dx ga) dx ga

＊

＊ T T T

＊ If:v1＝s ga，v2＝dx ga，v3＝v1 dx ＝S ga dx

＊ then:

＊ T T T

＊ (v2+ga v1)dx dx (v3+v3)

＊ S[k+1]＝S+----------------- - --------

＊ 2

＊ (v2) v2

＊

＊Note:SetPrecision[]is used to compensate for the gradual loss in

＊accuracy that Mathematica imposes on iterative calculations.

＊)

BFGS［x0_,vars_，gradients_，objectiveFunc_，epsi_，maxIterations_]:＝

Module[

{len，k，x，S，g，err，ds，a，ax，dx，xnew，g1，ga，v1，v2，v3，part0，part1，part2，

eTable，rTable，dxdxT，max，subst，residual}，

len＝Length[x0]；

k＝0；

(＊convert 1-D list to a 2-D column vector＊)

x＝Transpose[SetPrecision[{x0}，39]]；

(＊S[k]is an approximation to the inverse of the modified Hessian of

e2[x]at x＝x[k]＊)

S＝IdentityMatrix[len]；(＊S is a square matrix and always symmetric＊)

g＝GradFunc[x，vars，gradients]；(＊returns a gradient column vector＊)

err＝Norm[g]；(＊returns sqrt of sum of squares of elements of g＊)

eTable＝Table[0，{maxIterations}]；

rTable＝Table[-40.0，{maxIterations}]；

While[(err>epsi)&&(k<maxIterations)，

ds＝-S.g；(＊ds is a direction column vector＊)

(＊scalar′a′line search result＊)

ax＝FletcherLineSearch[x，ds，vars，gradients，objectiveFunc]；

a＝SetPrecision[ax，39]；

dx＝a＊ds；(＊dx is a column vector＊)

xnew＝x+dx；(＊new column vector result＊)

g1＝GradFunc[xnew，vars，gradients]；(＊column vector result＊)

ga＝g1-g；(＊also a column vector＊)

(＊dx.Transpose[dx] is a square symmetric matrix＊)

dxdxT＝dx.Transpose[dx]；

v1=s.ga；(＊also a column vector＊)

Page 4

Slavin.txt

v2[＝(Ttanspose[ga].dx)[[1]][[1]]；(＊v2 is scalar＊)

If[v2＝＝0.0，

Print[″BFGS error:v2＝″，v2]；

Print[″Iteration:″，k]；

Print[″v2＝″，v2]；

Print[″g＝″，g]；

Print[″S＝″，S]；

Print[″Determinant[S]＝″，Det[S]]；

Print[″Inverse[S]＝″，Inverse[S]]；

Print[″ds＝″，ds]；

Print[″a＝″，a]；

Print[″dx＝″，dx]；

Print[″xnew＝″，xnew]；

Print[″g1＝″，g1]；

Print[″ga＝″，ga]；

Abort[]；

]；

v3＝v1.Transpose[dx]；(＊v3 is a square matrix＊)

(＊Transpose[ga].v1 is a scalar，so part0 is scalar＊)

part0＝(v2+(Transpose[ga].v1)[[1]][[1]])/(v2^2)；

part1＝par0＊dxdxT；(＊part1 is square symmetric＊)

(＊any matrix added to its transpose is always square symmetric＊)

part2＝(v3+Transpose[v3])/v2；(＊part2 is square symmetric＊)

(＊S(old and new)is square symmetric＊)

S＝SetPrecision[S+part1-part2，39]；

x＝SetPrecision[xnew，39]；

g＝SetPrecision[g1，39]；

subst＝-Table[vars[[j]]->x[[j]][[1]]，{j，1，len}]；

residual＝Log[objectiveFunc/.subst]；

err＝Norm[dx]；

(＊Print[″iteration＝″，k，″，dx＝″，err]；＊)

(＊Print[″x＝″，Flatten[x]]；＊)

k+＝1；

eTable[[k]]＝err；

rTable[[k]]＝residual；

If[residual<-40.0，

Break[]；

]；

If［k＝＝maxIterations，

(＊Print[″BFGS Failed to converge.Final vector:″，x]；＊)

Return[]；

]；

(＊ Print[″Total iterations＝″，k，″，Line error＝″，err]；＊)

(＊ Print[″Finished:x＝″，Flatten[x]]；＊)

Return[Flatten[x]]；

]

(＊＊)

CorrectData[data_，filters_]:＝

Module[

{maxOrder，len，c，sum，j，nrSamples，samples，t}，

maxOrder＝Length[filters]；

len＝Length[filters[[1，1]]]；

c＝(len-1)/2；

sum＝0；

For[j=1，j<＝maxorder，++j，

sum+＝GenDistortion[data，filters[[j]]]；

]；

Page5

Slavin.txt

nrSamples＝Length[sum]；

Print[″nrSamples＝″，nrSamples]；

samples＝Take[data，{c+1，nrSamples+c｝]-sum；

Return[samples]；

]

(＊make real tap filter have the required DC＇delay′(in sample units，

where 0 is the first tap or sample.Usually DCdelay＝Length[tans]/2

is used＊)

DCadjust[realtaps_，DCdelay_]:＝

Module[

{newtaps，taps，i，sum1，sum2，offset，div}，

taps＝Length[realtaps]；

sum1＝Sum[irealtaps[[i+1]]，{i，1，taps-1}]；

sum2＝Sum[realtaps[[i+1]]，{i，0，taps-1}]；

div＝taps＊(DCdelay-(taps-1)/2)；

If[div＝＝0.0，

Print[″DCadjust divide by zero.″]；

Abort[]；

]；

offset＝(sum1-sum2＊DCdelay)/div；

Print[″offset＝″，offset]；

newtap＝realtaps+offset；

Return[newtaps]；

]

(＊find the′delay＇of a filter at DC＊)

DCdelay[realtaps_]:＝

Module[

{sumc，sumci，i，taps}，

(＊calculate the dc delay without Arg[]aliasing.＊)

taps＝Length[realtaps]；

sumc＝Sum[realtaps[[i+1]]，{i，0，taps-1}]；

If[sumc＝＝0.0，

Print[″DCdelay divide by zero.″]；

Abort[]；

]；

sumci＝Sum[realtaps[[i+1]]i，{i，1，taps-1}]；

Return[N[sumci/sumc]]；

]

Delays[list_，n_]:＝

Module[

{i，L，len，reals，angle，diff，offset，items，complex，

delay，radiansPerSample，prevAngle，newDelay}，

L＝Length[list]；

len＝n/2；

items＝len+1；

reals＝Join[list，Table[0，{i，1，n-L}]]；

complex＝Take[Fourier[reals，{FourierParameters->{-1，1}}]，items]；

delays＝Table[0，{items}]；

delays[[1]]＝DCdelay[list]；

For[i＝1，i<items，++i，

If[Abs[complex[[i+1]]]<10^-12，

delays[[i+1]]＝delays[[i]]，

(＊else＊)

radiansPerSample＝N[Pi＊i/len]；

delay＝delays[[i]]；

prevAngle＝delay＊radiansPerSample；

angle＝Mod[Arg[complex[[i+1]]]，pi，prevAngle-Pi/2]；

Page6

Slavin.txt

newDelay＝angle/radiansPerSample；

delays[[i+1]]＝newDelay；

]；

(＊ possible oversampling to interpolate＊)

For[i＝1，i<len，++i,

If[Abs[complex[[i+1]]]<10^-12,

delays[[i+1]]＝(delays[[i]]′+delays[[i+2]])/2；

]；

Return[delays-(L-1)/2]；

]

(＊

＊ design filters to correct ADC output.

＊ DesignFilters[] works on the output of ProcessCalData

＊)

DesignFilters[data_，correctionFilterLen_]:＝

Module[

{fi，samples，freqs，start，end，skip，i，j，freqIndex，ampl，temp，cycle，

refDlyA，refDlyB，dly,qq,K，logsize,logFreqs，mags，parts，filters,

nrVars，gainA，gainB,delays，totalGains，xx，sortList，sorted

prev，next，permute，undef，modified，odd}，

xx＝data；

freqs＝Length[xx]-1；(＊input is oflength 2＾m+1＊)

logFreqs＝Log[2，freqs]；

If[！IntegerQ[logfreqs]，

Print[Length[xx]-1 is not a power of 2.″]；

Abort[]；

]；

logSize＝logFreqs+3；

samples＝2^logSize；

skip＝samples/(2＊freqs)；

Print[″skip＝″，skip]；

Print[″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊ correction phase ＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊″]；

undef＝-2；

start＝end＝undef；

For[fi=0，fi<＝freqs，++fi，

i＝fi+1；

If[Length[xx[[i]]]＝＝0，

If[(start！＝undef)&&(end＝＝undef)，

end＝fi-1；

]；

，(＊else＊)

If[start＝＝undef，

start＝fi；

]；

order＝0；(＊default is set later＊)

Print[″start＝″，start，″，end＝″，end]；

For[fi＝start，fi<＝end，++fi，

i＝fi+1；

freqIndex＝fi＊skip；

cycle＝N[samples/freqIndex]；

If[fi＝＝start，

nrVars＝Length[xx[[i]]]；

order＝(nrVars-1)/2；

page7

Slavin.txt

odd＝False；

，(＊else＊)

(＊a previous result exists＊)

prev＝Table[xx[[i-1，order+1+j]]，{j，1，order}]；

next＝Tabl[xx[[i，order+1+j]]，{j，1，order}；

{permute，modified，odd}＝MatchMod[prev，next，cycle]；

xx[[i]]＝Flatten[{xx[[i，1]]，

Table[xx[[i，1+permute[[j]]]]，{j，1，order}]，

modified}]；

]；

If[xx[[i，1]]<0.0

xx[[i，1]]＝-xx[[i，1]]；

odd＝！odd；

]；

If[odd，

If[xx[[i，order+2]]>0.0，

xx[[i，order+2]]-＝cycle＊0.5，

xx[[i，order+2]]+＝cycle＊0.5

]；

(＊ Print[″end:xx[[″，i，″]]＝″，xx[[i]]]；＊)

]；

(＊initialize uncalibrated regions to allow later matrix transpose＊)

For[fi＝0，fi<start，++fi，

i＝fi+1；

xx[[i]]＝Table[0，{2＊order+1}]；

]；

For[fi＝end+1，fi<＝freqs，++fi，

i＝fi+1；

xx[[i]]＝Table[0，{2＊order+1}]；

]；

qq＝Transpose[xx]；

K＝qq[[1]]；

Print[″K:″]；

PlotFreq[K，start，end]；

For[j＝1，j<＝order.++j，

Print[″Amplitude″，j]；

PlotFreq[qq[[1+j]]，start，end]；

]；

For[j＝1，j<＝order，++j，

Print[″Delay″，j]；

PlotFreq[qq[[order+1+j]]，start，end]；

]；

(＊gainA^(order-1)＊gainB＝K，gainA applies to allfilters except one

which uses gainB＊)

totalGains＝Take[K，{start+1，end+1}]；

gainA＝Abs[totalGains]＾(1/order)；

gainB＝totalGains/(gainA^(order-1))；

delays＝Table[Take[qq[[j+order+1]]，{start+1，end+1}]，{j，1，order}]；

filters＝Table[Fj]terDesign[If[j＝＝order，gainB，gainA]，delays[[j]

freqs，start，correctionFilterLen]，{j，1，order}]；

For[j＝1，j<＝order，++j，

ListPlot[filters[[j]]，{Plotjoined->True，PlotRange->All}]；

]；

mags＝Table[Magnitude[filters[[i]]，samples]，{i，1，order}]；

parts＝Table[Take[mags[[j]]，{1，1+freqs＊skip，skip}],{j，1，order}]；

Print[″Original Combined amplitude response:″]；

PlotFreq[K，start，end]；

Print[econstructed combned response:″]；

Page8

Slavin.txt

ampl＝Product[parts[[i]]，{j，1，order}]；

PlotFreq[amp]，start，end]；

Print[″Amplitude response difference:″]；

PlotFre[K-amp]，start，end]；

refDlys=Table[GenDelay[fi]ters[[j]]，2＊freqs]，{j，1，order}]；

For[j＝1，j<＝order，++j，

Print[″Reconstructed delay response″，j，″:″]；

PlotRawDelay[refDlys[[j]]]；

]；

Return[filters]；

]

(＊single frequency DFT with gain normalization，returns a complex value＊)

DFT[data_，m_]:＝

Module[

{len，sum，t}，

len＝Length[data]；

sum＝0；

For[t＝0，t<len，++t，

sum+＝(Exp[I＊2＊pi＊m＊t/len] data[[t+1]])；(＊1-based array indexing＊)

]；

Return[2sum/len]；

]

(＊

＊Generates a broad-spectrum set of distorted calibration data.

＊′filters′＝{filter2，filter3，...}

＊where filter2＝{{<taps>}，{<taps>}

＊wherefilter3＝{{<taps>},{<taps>}，{<taps>}}

＊and each list<taps>has the same length′len′(must be odd).

＊For example，before runnning Distortcal[]，with

＊′distortFilterLen＝63；′and′spread＝4.0；′

＊filter2＝GenFilters[distortFilterLen，2，spread]；

＊filter3＝GenFilters[distortFilterLen，3，spread]；

＊filters＝{filter2，filter3}；

＊′sx′is the output data record length for each test frequency pair.

＊′freqs′is one more than the number of output data frequency

＊sets，e.g.for 127 sets，freqs＝128.

＊Note sx/freqs must be integer.

＊′distortion is the distortion level in output lsb units

＊′signal′is the two-tone amplitude in output lsb units.

＊′noise′is the added noise standard deviation in output lsb units.

＊Values that match the records acquired for the AD6645 are；

＊DataGen[filters，1024，128，0.04，15000.0，2.0]；

＊)

DistortCal[fiters_，sx_，freqs_，distortion_，signal_，nois_]:＝

Module[

{error，order，q，span，data，i，j，len，qfreq，f1，f2，phase1，phase2，

sines，samples，t，c，sum，maxOrder}，

maxOrder＝Lenqth[filters]；

len＝Length[filters[[1，1]]]；

error＝False；

If[！IntegerQ[sx]，

Print[″Analysis length is not integer.″]；

error＝True；

]；

If[！IntegerQ[freqs]，

Print[″Frequency range is not integer.″]；

error＝True；

Page9

Slavin.txt

]；

If[！IntegerQ[len]，

Print[″Filter lengths are not integer."]；

error＝True；

]；

If[Mod[sx，4]！＝0，

print[″Analysis length is not a multiple of 4.″]；

error＝True；

]；

If[Mod[sx，freqs]！＝0，

Print[″Analysis length is not a multiple of the number of″

″frequencies.n″]；

error＝True；

]；

If[！OddQ[len]，

Print[″Filter lengths are not odd.″]；

error＝True；

]；

If[error，

Abort[]；

]；

c＝(len+1)/2；

q=2.0＊Pi/sx；

span＝sx/(2＊freqs)；

data＝Table[0，{freqs-1}]；

For[i＝1，i<freqs，++i，

qfreq＝span＊i；

f1＝qfreq-1；

f2=qfreq+1；

phase1＝Random[Rea]，{0，2＊pi}，10]；

phase2＝Random[Rea]，{0，2＊pi}，10]；

(＊Length[sines]＝sx+len-1＊)

sines＝Table[(Cos[q＊f1＊t-phase1]+Cos[q＊f2＊t-phase2])，

{t，0,sx+len-2}]；

(＊Length[samples]＝Length[sines]+1-len＝sx＊)

sum＝0；

For[j＝1，j<=maxOrder，++j，

sum+＝GenDistortion[sines，filters[[j]]]；

]；

samples＝distortion＊sum；

samples+＝Table[Gaussian[0，noise]，{t，1，sx}]；

samples+＝signal＊Take[sines，{c，sx+c-1}]

(＊output highest frequencies first＊)

data[[Freqs-i]]＝samples；

]；

Return[data]；

]

(＊＊)

DistortData[data_，filters_，distortion_，signal_，noise_]:＝

Module[

{maxOrder，len，c，sum，j，nrSamples，samples，t}，

maxOrder＝Length[filters]；

len＝Length[filters[[1，1]]]；

c＝(len-1)/2；

sum＝0；

For[j＝1，j<＝maxOrder，++j，

sum+＝GenDistortion[data，filters[[j]]]；

]；

nrSamples＝Length[sum]；

Page10

Slavin.txt

Print[″nrSamples＝″，nrSamples]；

samples＝distortion＊sum；

If[noise！＝0，

samples+＝Table[Gaussian[0，noise]，{t，1，nrSamples}]；

]；

samples+＝signal＊ Take[data，{c+1，nrSamples+c}]；

Return[samples]；

]

(＊

＊′symbVars′is the list of symbolic variables；

＊ Flatten[{k，Table[a[j]，{j，1，order}]，Table[d[j]，{j，1，order}}]；

＊′ix′is the index symbol to use(passed in at a higher level).

＊′sx′is the symbolic time scale factor.

＊′dft′is a symbolic variable that is expanded to a list，used internally

＊ and also returned as the list′dfts′.

＊Returns；

＊′symbolicBins′:a list of IM and harmonic frequences as a functi on of the

＊ frequency index′ix′.

＊ ′e2′: a symbolic error equation as a function of ix，sx，symbVars

＊ and′dfts′.

＊′g＇:a list of symboli gradient equations that are derivatives of

＊ ′e2′with respect to each of the′symbVars＇variable list.

＊′dfts′:a list of symbolic variables:{dft[1]，dft[2]，...}

＊)

DoInit[symbVars_，ix_，sx_，dft_]:＝

Module[

{e2，g，j，nrAnalysisBins，symbolicBins，nrVars，dfts，dft，p}，

(＊

＊Create the symbolic expressions for the objective function e2

＊and all its partial derivatives with respect to ＇vars′.

＊Note:only really needs to be called once for multiple analysis

＊frequencies，as the results are purely symbolic.

＊)

{symbolicBins，p}＝Terms[ix，sx，symbVars]；

nrAnalysisBins＝Length[symbolicBins]；

dfts＝Table[dft[j]，{j，1，nrAnalysisBins}]；

(＊e2 is the symbolic sum of the squares of the differences between the

predicted vlues′p′and′dfts＇at each analysis bin＊)

e2＝Sum[(p[[j]][[1]]-Re[dfts[[j]]])^2+

(p[[j]][[2]]-Im[dfts[[j]]])^2，{j，1，nrAnalysisBins}]；

(＊calculate symbolic partial derivatives＊)

nrVars＝Length[symbVars]；

g＝Table[D[e2，symbVars[[j]]]，{j，1，nrvars}]；

Return[{symbolicBins，e2，g，dfts}]；

]

(＊

＊′e2′is symbolic error equation.

＊′g′is symbolic gradiente equation.

＊′dfts′is array of symbolic dft results.

＊′symbVars′is list of veriables that BGGS is attempting to solve.

＊′y′is list of dft results.

＊′x′is set of start values for BFGS.

＊′sx′is samples symbol ′s′is actual sample value.

＊′ix′is frequency index symbol，′i′is actual index value.

＊′acc′ is BFGS accuracy goal.

＊′maxIterations′limits BFGS search time.

Page11

Slavin.txt

＊)

Dosolve[e2_，g_，dfts_，symbVars_，y0_，x_，sx_，s_，ix_，i_，acc_，

maxIterations_]:＝

Module[

{e2n，gn，j，nrAnalysisBins，subst，substx，xnew，subst2，residual，len，y}，

y＝SetPrecision[y0，39]；

nrAnalysisBins＝Length[y]；

(＊evaluate the symbolic expressions with the known constants＊)

substx＝Table[dfts[[j]]->y[[j]]，{j，1，nrAnalysisBins}]；

subst＝Flatten[{sx->s，jx->i，substx}]；

e2n＝e2/.subst；

gn＝g/.subst；

xnew=BFGS[x,symbVars，gn，e2n，acc，maxIterations]；

len＝Length[xnew]；

subst2=Table[symbVars[[j]]->xnew[[j]]，{j，1，len}]；

residul＝e2n/.subst2；

Return[{xnew，residual}]；

]

ExpFunc[q_]:＝If[q<-100，0，Exp[q]]

(＊takes sparce gain and delay lists，places them in a′size＇length list

at″start′(zero based)，for iterative　filter calculations＊)

FilterDesign[gain_，dlyin_，size_，start_，requiredtaps_]:＝

Module[

{l1，l2，ampl，delay，dly，n，taps，nrtaps，offset，i，j，

fftsize，deviation，avgdev，min，bestTaps}，

l1＝Length[gain]；

l2＝Length[dlyin]；

If[l1！＝l2,

Print[″Missmatched gain and delay list lengths.″]；

Abort[]；

]；

If[EvenQ[requiredtaps]，

Print[″Number of taps must be odd.″]；

Abort[]；

]；

If[(start<＝0)||((start+l1)>＝size)

Print[″start and size values illegal for″，l1，″gains and delays.″]；

Abort[]；

]；

dly＝dlyin+size；(＊center the delays＊)

n＝size+1；

ampl＝Flatten[{Table[0，{start}]，gain，Table[0，{n-start-l1}]}]；

delay＝Flatten[{Table[size，{start}]，dly，Table[size，{n-start-l2}]}]；

If[Length[ampl]！＝n，

Print[″Internalerror:Length[ampl]＝″，Length[ampl]，″，n＝″，n]；

Abort[]；

]；

If[Length[delay]！＝n，

Print[″Internal error:Length[delay]＝″，Length[delay]，″，n＝″，n]；

Abort[]；

]；

fftsize＝2＊size；

offset=Floor[reqiredtaps/2]；

min＝Infinity；

For[j＝0，j<100，++j，(＊note:Break[] for small j in loop＊)

taps＝Taps[ampl，delay]；

nrtaps＝Length[taps]；

(＊set coeffcients outside desired filter size to zero

Page12

Slavin.txt

For[i＝0，i<fftsize，++i，

If[(i<(size-offset))||(i>＝(size+offset))，

taps[[i+1]]＝0；

]；

{ampl，delay}＝Response[taps]；

(＊ overwrite gain and delay with original values in critical

region＊)

deviation＝0.0；

For[i＝start，i<start+l1，++i，

deviation+＝Abs[delay[[i+1]]-dly[[i-start+1]]]；

]；

avgdev＝deviation/(l1+1)；

If[avgdev<min，

min＝avgdev；

bestTaps＝taps；

]；

Print[″Loop″，j，″:average absolute delay deviation＝″，avgdev]；

If[j＝＝25,

Break[]；

]；

For[i＝start，i<start+l1，++i，

ampl[[i+1]]＝gain[[i-start+1]]；

delay[[i+1]]＝dly[[i-start+1]]；

]；

For[i＝0，i<start，++i，

delay[[i+1]]＝delay[[start+1]]；

]；

taps＝Table[bestTaps[[i+1]]，{i，size-offset，size+offset}]；

dcdelay＝DCdelay[taps]；

margin＝Floor[0.45 requiredtaps]；

If[dcdelay<margin，

Print[″limiting lower，″]；

taps＝DCadjust[taps，margin]；

{ampl，delay}＝Response[taps]，

(＊else＊)

If[dcdelay>＝requiredtaps-margin，

Print[″limiting upper，″]；

taps＝DCadjust[taps，requiredtaps-margin-1]；

{ampl，delay}＝Respons[taps]；

]；

Return[taps]；

]

(＊x and ds are 2-D column vectors，′vars′is a list of variables in the

gradients(derivatives of the objective function).Return the objective

function value at vars＝x＊)

FletcherLineSearch[x_，ds_，vars_，gradients_，objectiveFunc_]:＝

Module[

{m，tau，rho，sigma，gamma，mhat，epsilon，xk，s，f0，gk，deltaf0，dk，aL，aU，fL，dfL，

a0，a0hat，a0Lhat，a0Uhat，gtemp，df0，deltaa0，alpha，div}，

m＝0；

sigma＝0.1；

tau＝Min[sigma，0.1]：(*tau<＝sigma and tau<＝0.1＊)

rho＝sigma；(＊rho<＝slgma＊)

gamma＝0.75；

mhat＝400；

epsilon＝1.0＊10^-10；

xk＝x；(＊xk is a column vector＊)

s＝ds；(＊′s′a column vector＊)

Page13

Slavin.txt

(＊step 1: evaluate fk and gk at k＝0.f0 is a scalar，g0 is a column

vector＊)

f0＝ObjectiveFunc[xk，vars，objectiveFunc]；

gk＝GradFunc[xk，vars，gradients]；(＊gk is acolumn vector＊)

m+＝2；

deltaf0=f0；

(＊step2:initialize the search＊)

dk＝s；(＊dk is a column vector＊)

aL＝0.0；

aU＝1.0＊＊10^99；

fL＝f0；

dfL＝(Trranspose[gk].dk)[[1]][[1]]；

If[Abs[dfL]>epsilon，

If[div＝＝0.0，

Print[″FletcherLineSearch divide by dfL＝0.0″]；

Abort[]；

]；

a0＝-2.0＊deltaf0/dfL；

，

a0＝1.0；

]；

If[(a0×10^-9)||(a0>1.0)，

a0＝1.0；

]；

(＊step3: ＊)

While[True，

deltak＝a0＊dk；(＊deltak is a column vector＊)

f0＝ObjectiveFunc[xk+deltak，vars，objectiveFunc]；

m+＝1；

If[(f0>(fL+rho＊(a0-aL)＊dfL))&&

(Abs[fL-f0]>epsilon)&&(m<mhat)，

If[a0<aU，

aU＝a0；

]；

div＝2.0＊(fL-f0+(a0-aL)＊dfL)；

If[div＝＝0.0，

Print[″FletcherLineSearch divide by div＝0.0″]；

Abort[]；

]；

a0hat＝aL+(dfL＊(a0-aL)^2)/div；

a0Lhat＝aL+tau＊(aU-aL)；

If[a0hat<a0Lhat，

a0hat＝a0Lhat；

]；

a0Uhat＝aU-tau＊(aU-aL)；

If[a0hat>a0Uhat，

a0hat=a0Uhat；

]；

a0＝a0hat；

(＊else＊)

(＊gtemp is a vector＊)

gtemp＝GradFunc[xk+a0＊dk，vars，gradents]；

df0＝(Transpose[gtemp].dk)[[1]][[1]]；

m+＝1；

If[(df0<sigma＊dfL)&&(Abs[fL-f0]>epsilon)&&

(m<mhat)&&(dfL！＝df0)，

div＝dfL-df0；

deltaa0＝(a0-aL)＊df0/div；

If[deltaa0<＝0，

a0hat＝2＊a0；

，(＊else＊)

a0hat＝a0+deltaa0；

]；

Page14

Slavin.txt

a0Uhat＝a0+gamma＊(aU-a0)；

If[a0hat>a0Uhat，

a0hat＝a0Uhat；

]；

aL＝a0；

a0＝a0hat；

fL＝f0；

dfL＝df0；

，(＊else＊)

Break[]；

]；

If[a0<1.0＊10＾-5，

alpha=1.0＊10^-5；

，(＊else＊)

alpha＝a0；

]；

Return[alpha]；

]

FTSize[len_，min_]:＝

Module[

{}，

Return[Max[min，2^Ceiling[Log[2，len]]]]；

]

Gaussian[mean_，sd_]:＝

Module[

{result}，

Return[sd＊Sqrt[2.0]*InverseErf[Random[]]＊Sign[Random[]-0.5]+mean]；

]

(＊finds sample delay as a function of frequency.

uses ALL′list′taps-no symmetry is assumed＊)

GenDelay[list_，min_:16]:＝

Module[

{len，delays，endindex，nx}，

len＝Length[list]；

nx＝FTSize[len，min]；

endindex＝nx/2；

Return[Take[Delays[list，nx]，{1，endindex+1}]]；

]

(＊

＊distort the data for a particular filter order Length[taps].

＊)

GenDistortion[data_，taps_]:＝

Module[

{len，filters，flen，delay，output，i，file}，

len＝Length[data]；

filters＝Length[taps]；

flen＝Length[taps[[1]]]；

If[EvenQ[flen]，

Print[″Filter lengths are not all odd.″]；

Abort[]；

]；

For[i＝1，i<＝taps1，++i，

len1=Length[taps[[i]]]；

Page15

Slavin.txt

If[len1！＝flen，

Print[″Filter lengths are not all equal.″]；

Abort[]；

]；

delay＝(1+flen)/2；

output＝Table[1.0，{len-flen+1}]；

For[i＝1，i<＝filters，++i，

filt＝taps[[i]]；

If[Length[filt]！＝flen，

Print[″Filter lengths are not all the same.″]；

Abort[]；

]；

output＊＝ListConvolve[filt，data]；

]；

Return[output]；

]

(＊generate a set of filters for a particular order of distortion＊)

GenFilters[len_，order_，spread_]:＝

Module[

{error，taps，i，j，c}，

error＝False；

If[！IntegerQ[len]，

Print[″Filter lengths are not integer.″]；

error＝True；

]；

If[！OddQ[len]，

Print[″Filter lengths are not odd.″]；

error=True；

]；

If[error，

Abort[]；

]；

c＝(len+1)/2；

taps＝Table[Gaussian[0，1]＊ExpFunc[-Abs[(j-c)＊2/spread]]，

{i，1，order}，{j，1，len}]；

Return[taps]；

]

(＊

＊′samples′gives nr of output sampels(plus extra)，and also fundamental

＊length，so all generated sines are harmonics of the fundamental.

＊′extra′gives the number of extra samples in the output record that may be

＊needed to compensate for convolution data loss due to filters.

＊′freqLo′is the low normalized frequency for sine output:

＊0<freqLow<0.5.The nearest higher harmonic of the′samples′frequency

＊is the lowest output frequency.

＊′freqHi′is the high normalized frequency for sine output:

＊0<freqLow<0.5.The nearest lower harmonic of the′samples′frequency

＊is the higheest output frequency.

＊′freqInterval′is the normalized frequency step，also a harmonic of the

＊′samples′frequency.

＊)

GenTest[samples_，extra_，freqLo_，freqHi_，freqInterval_]:＝

Module[

{kInterval，kStart，kEnd，sum，x，j，t}，

kInterval＝Floor[freqInterval＊samples+0.5]；

Page16

Slavin.txt

kStart＝Ceiling[freqLo＊samples]；

kEnd＝Floor[freqHi＊samples]；

sum＝0；

x=2.0＊Pi/samples；

For[j＝kStart，j<＝kEnd，j+＝kInterval，

sum+＝Table[sin[x＊j＊t]，{t，0，samples+extra-1}]；

]；

Return[sum]；

]

(＊′x′is a 2-D column vector，′vars′is a list of variables in the list of

′gradients′functions which are the partial derivatives of the objective

function with respect to each variable.＊)

GradFunc[x_，vars_，gradients_]:＝

Module[

{g，subst，j，len}，

len＝Length[x]；

subst＝Table[vars[[j]]->x[[j]][[1]]，{j，1，len}]；

g＝gradients/.subst；

Return[Transpose[{g}]]；(＊return a column vector as a 2-D matrix＊)

]

IsSame[xx_]:＝

Module[

{nrVars，order，indices，permArray，perms，same，i，j，perm，v0，v1}

nrVars＝Length[xx]；

order＝(nrVars-1)/2；

indices＝Table[i，{i，1，order}]；

permArray＝Permutations[indices]；

perms＝Length[permArray]；

same＝False；

For[j＝1，j<＝perms，++j，

perm=permArray[[j]]；

v0＝perm[[1]]；

v1＝pem[[2]]；

If[v1×v0，

same＝(Chop[xx[[v0+order+1]]-xx[[v1+order+1]]]＝＝0)；

If[same，

Break[]；

]；

Return[same]；

]

Log10List[list_]:＝

Module[

{minv，lminv，log，i，item，L，result}，

L＝Length[list]；

minv＝Min[list]；

If[minv<＝10^-10，

minv＝10^-10；

]；

result＝TaBle[0，{L}]；

For[i＝0，i<L，++i，

item＝list[[i+1]]；

If[item<minv，

item＝minv；

]；

result[[i+1]]＝20.0 Log[10，item]

Page17

Slavin.txt

]；

Return[result]；

]

LogFFT[list_]:＝

[Mdule[

{i，L，magfreg，logmagf，data，bins}，

L＝Length[list]；

bins＝Floor[(L+1)/2]；

magfreq＝Abs [Take[Fourier[list，{FourierParameters->{-1，1}}]，bins]]；

freq=Table[N[i/L]，{i，0，bins-1}]；

logmagf＝Log10List[magfreq]；

data=Transpose[{freq，logmagf}]；

Return[data]；

]

Magnitude[list_，n_]:＝

Module[

{i，L，reals，magfreq}，

L＝Length[list]；

reals＝Join[list，Table[0，{i，1，n-L}]]；

magfreq＝Abs[Take[Fourier[reals，{FourierParameters->{1，1}}]，n/2+1]]；

Return[magfreq]；

]

(＊

＊patt1 and patt2 numeric list lengths are equal.Find the permutation index

＊of patt2 that minimizes the sum of the squares of the differences of patt1

＊and patt2，modulo：'mod’.-Used for reorganizing phases.Also returns a

＊modified version of patt2 which is close to patt1，modulo′mod′。

＊For example:

＊MatchMod[{0.2，0.5，0.8}，{1.22，0.88，0.47}，1]

＊yields the indicesmodified patt2，and odd status:

＊{{1，3，2}，{0.22，0.47，0.88}，False}

＊)

MatchMod[patt1_，patt2_，mod_]:＝

Module[

{len1，len2，indices，permArra，perms，min，sum，perm，i，j，diff，index，x，new，

step，steps，total，odd}，

len1=Length[patt1]；

len2＝Length[patt2]；

If[len1！＝len2，

Print[″Length mismatch.″]；

Abort[]；

]；

indices＝Table[i，{i，1，len1}]

permArray＝Permutations[indices]；

perms＝Length[permArray]；

min＝Infinity；

step＝mod＊0.5；

For[j＝1，j<＝perms，++j，

sum＝0；

perm＝permArray[[j]]；

total=0；

For[i＝1，i<＝len1，++i，

diff＝patt1[[i]]-patt2[[perm[[i]]]]；

steps＝Floor[diff/step+0.5]；

total+＝steps；

sum+＝(diff-step＊steps)^2；

Page18

Slavin.txt

]；

If[sum<min，

(＊found best solution so far＊)

odd＝OddQ[total；(＊record the odd-ness of the solution＊)

index＝j；

min＝sum；

]；

perm＝permArray[[index]]；

new＝Table[x＝patt2[[perm[[i]]]]；

x+step＊Floor[(patt1[[i]]-x)/step+0.5]，

{i，1，len1}]；

Return[{perm，new,odd}]；

]

(＊′x′is a 2-D column vector，′vars′is a list of all variables in the

objective function.Return the objective function value at vars＝x＊)

ObjectiveFunc[x_，vars_，objectiveFunc_]:＝

Module[

{subst，j，len}，

len＝Length[x]；

subst＝Table[vars[[j]]->x[[j]][[1]]，{j,1，len}]；

Return[objectiveFunc/:subst]；(＊return a scalar＊)

]

PlotColor[list_，color_]:＝

Module[

{min，max，r,g，b，data}，

r＝color[[1]]；

g＝color[[2]]；

b＝color[[3]]；

data＝{0，Transpose[list][[2]]};＊always include zero＊)

max＝Max[data]；

min＝Min[data]；

Return[ListPlot[list，{PlotJoined->True，

GridLines->Automatic，

PlotRange->{min，max}，

ImageSize->600，

PlotStyle->RGBColor[r，g，b]}]]；

]

(＊startIndex and endIndex indices(0-based)apply to the list＊)

plotFreq[list_，startIndex_:0，endIndex_:-1]:＝

Module[

{i，L，delays，freq，data，idata，endi，max，min，selected}，

L＝Length[list]；

If[endIndex＝＝-1，

endi＝L-1，

endi＝endIndex

]；

freq＝Table[N[i/(2＊(L-1))]，{i，startIndex，endi}]；

selected＝Take[list，{startIndex+1，endi+1}]；

data＝Transpose[{freq，selected}]；

idata＝{0，selected}；(＊always include zero＊)

max＝Max[idata]；

min＝Min[idata]；

Return[ListPlot[data，{PlotJoined->True，

GridLines->Automatic，

PlotRange->{{0，0.5}，{min，max}}，

ImageSize->600，

Page19

Slavin.txt

PlotStyle->RGBColor[1，0，0]}]]；

]

(＊finds sample delay as a function of frequency.

uses ALL ′list′taps -no symmetry is assumed，plot in red＊)

PlotRawDelay[delays_]:＝

Module[

{i，len，freq，data}，

len＝Length[delays]-1；

freq＝Table[N[i/(2＊len)]，{i，0，len}]；

data＝Transpose[{freq，delays}]；

Return[PlotColor[data，{1，0，0}]]；

]

PlotUnwindowedFFT[list_]:＝

Module[

{}，

Return[WfmColor[LogFFT[list]，{1，0，0}]]；

]

PolarToComplex[ampls_，angles_]:＝

Module[

{l1，l2，fl，i，len，reals，angle，ampl，real，imag，out}，

l1＝Length[ampls]；

l2＝Length[angles]；

If[l1！＝l2，

Print[″Missmatched amplitude and angle list lengths.″]；

Abort[]；

]；

len＝l1-1；

fl＝2＊len；

out＝Table[0，{i，1，fl}]；

For[i＝0，i×len，++i，

angle＝angles[[i+1]]；

ampl＝ampls[[i+1]]；

real＝ampl Cos[angle]；

imag＝ampl Sin[angle]；

out[[i+1]]＝real+I imag；

If[(i ！＝0)，

out[[fl-i+1]]＝real-I imag；

]；

out[[len+1]]＝ampls[[len+1]]；

Return[out]；

]

PrAmplPhase[data_，bin_]:＝

Module[

{out，ampl，ph}，

out＝DFT[data，bin]；

ampl＝Abs[out]；

ph＝Arg[out]；

Return[{ampl，ph}]；

]

ProcessCalData[data_，order_，lowf_，hiF_，acc_]:＝

Module[

{fi，

Page20

Slavin.txt

delta，freqs，start，end，skip，i，j，t，freqIndex，overhang，

i1，i2，ampl，ph，len，bins，sBins，w1，w2，scaling，temp，cycle，

amplRet，phRef，amplRatio，q，scale，

amplitude2，relPhase2，relDelay2，amplitude3，relPhase3，relDelay3，

corrPh，turns，rph，xbin，analyze1，reference，redo，d1，d2，x，y，

dly，dlyA，dlyB,amplA，amplB,donesomething，

descr，expected，K，logSize，logFregs，amplx，

magA，magB，partAm1，partAp1，partBm1，partBp1，transpose，

delays，base，partial，de，amplI，ii，nrBins，

filter1，filter2，frBins，index，f1，f2，fs，out1，out2，out3，

e2，g，s，nrVars，symbolicBins，symbVars，

k，a，d，ix，sx，xnew，yl，allSymBins，blen，slen，numericBins，

dfts，dft，results，maxorde，gainA，gainB，delavA，delayB，

maxIterations，substx，jj，kk，random，

residual，lowestResiduals，bestResults，startPoints，xStart，gains，

solutions，found，xx，qresults，solution,smallesttResidueIndex，

sortList，sorted，state，const，maxTries，extraTries，

symbvarsW，symbolicBinsW，e2W，gW，dftsW}，

Print[″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊″]；

Print[″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊ initialization ＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊″]；

Print［″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊″]；

(＊Various″magic″parameters＊)

(＊maximum tries to find a first BFGS solution.Too large means time

may be wasted finding a solution when none can be found.Too small

and a solution may not be found.＊)

maxTries＝100；

(＊extra attempts to find more BFGS solutions after the first is found.

Too large and time may be wasted finding additional solutions when

none exist. Too small，and a better solution can be missed.The

probability of missing a second solution is about 1/2 per try，

so the probability of missing it is(1/2)^extraTries＊)

extraTries＝25；

(＊maximum BFGS iterations allowed in which to converge to a solution.

Too large and time may be wasted trying to find a BFGS solution，

too small and a solution may not be found.Typically，2nd order

needs<30 and 3rd order needs<80.＊)

maxIterations＝100；

(＊maximum number of total unique solutions at a cal frequency.

2nd order:3^2＝9 solutions per group，3rd order:3^3＝27 per group.

Too larg a value and more memory is used.Too small and there is

insufficient storage for multiple groups.I′ve never seen more than 3

groups，or 3＊27＝81 solutions for any frequency.＊)

solutions＝200；

(＊sampling frequency(MHz)＊)

fs＝100.0；

overhang＝63；

logFreqs＝7；

symbVars＝

Flatten[{k，Table[a[j]，{j，1，order}]，Table[d[j]，{j，1，order}]}]；

{symbolicBins，e2，g，dfts}＝DoInit[stmbvars，ix，sx，dft]；

If[(lowF<0)||(lowF>0.5)||(hiF<0)||(hiF>0.5)||

(lowF>＝hiF)，

Print[″low or high frequency bounds are<0 or×0.5，or″]；

Print[″low frequency bound″，lowF，″is>＝hidh bound″，hiF]；

Abort[]；

]；

Print[″sampling frequency:″，fs，″MHz″]；

logSize＝logFreqs+3；(＊added value must be>＝3＊)

Page21

Slavin.txt

delta＝1；

samples＝2^logSize；(＊power of 2＊)

n＝(2＊overhang+1)；(＊n is odd filter length＊)

len＝samples+2＊(n-1)；(＊overhang required at each end＊)

If[logFreqs>logSize-3，

Print[″logFreqs must be<＝logSize-3″]；

Abort[]；

]；

freqs＝2^logFregs；(＊frequencies to Nyguist. A power of 2＊)

start＝Max[Floor[lowF＊freqs＊2]，1]；

end=Min[Ceiling[hiF＊freqs＊2]，freqs-1]；

If[start＝＝freqs/2，

Print[″Error:start frequency evaluates to fs/2.This cannot be″

″used because it starts on an alias notch.″]；

Abort[]；

]；

If[end＝＝freqs/2，

Print[″Error:end frequency evaluates to fs/2.This cannot be″

″used because it starts on an alias notch.″]；

Abort[]；

]；

skip＝samples/(2＊freqs)；

Print[″Frequency indices:start＝″，start，″，end＝″，end]；

Print[″Range:0:″，freqs-1，″，index skip＝″，skip]；

(＊nr of different IM and harmonic frequency bins to consider＊)

nrBins＝Length[symbolicBins]；

amplitude＝Table[0，{i，1，freqs+1}，{j，1，nrBins}]；

yl＝Table[0,{i，1，freqs+1}，{j，1，nrBins}]；

bins＝Table[0，{i,1，freqs+qs+1}，{j，1，nrBins}]；

amplRatio＝Table[0.0，{i，1，freqs+1}]；

Scaling＝Table[0.0，{i，1，freqs+1}]；

descr＝Table[0.0,{i,1,freqs+1}]；

allSymBins=AllTerms[ii，orde]；

blen＝Length[allSymBins]；

dly＝0.0；

redo＝-1；

const＝2＊Pi/samples；

Print[″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊ DFT analysis ＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊″]；

Print[″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊″]；

(＊

＊default x[[i]]＝False indicates not measurable due to aliasing，

＊True indicates measurable but failed to converge on search，

＊otherwise x[[i]]equals a list of solutions.

＊)

x＝Table[False，{i，1，freqs+1}]；

For[fi＝start，fi<＝end，++fi，

i＝fi+1；(＊offset for non-zero based array indexing＊)

index＝freqs-fi；

analyze1＝data[[index]]；

freqIndex＝fi skip；

i1＝freqIndex-delta；

i2＝freqIndex+delta；

w1＝2＊Pi＊i1/samples；

w2＝2＊Pi＊i2/samples；

Page 22

Slavin.txt

Print[″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊″]；

Print[″fi＝″，fi，″，index＝″，index，″，freqIndex＝″，freqIndex]；

(＊look for alias collisions using all frequencies＊)

numericBins＝allSymBins/.ii->freqIndex；

sBins＝Sort[AliasF[numericBins，samples]]；

Print[″sBins＝″，sBins]；

For[j＝1，j<blen，++j，

If[sBins[[j]]＝＝sBins[[j+1]]，

redo＝i；

Break[]；

]；

If[redo＝＝i，

Continue[]

]；

x[[i]]＝True；(＊indicate measurable at i＊)

(＊symbolicBins represents only the frequencies used for analysis，

which excludes Dc and the ibput cal frequencies themselves＊)

bins[[i]]＝symbolicBins/.ix->freqIndex；

Print[″bins″，order，″＝″，bins[[i]]]；

(＊generate timing references from signal+small distortion＊)

reference＝analyze1^order；

temp＝Abs[DFT[analyze1，i1]]；

temp+＝Abs[DFT[analyze1，i2]]；

scaling[[i]]＝temp＊0.5；

(＊calculate amplitudes，phases and delays＊)

For[j＝1，j<＝nrBins，++j，

xbin＝bins[[i，j]]；

temp＝DFT[analyze1，xbin]；

{amplRef，phRef＝PrAmplPhase[reference，xbin]；

yl[[i，j]]＝temp＊Exp[-I＊phRef]；

]；

f1＝N[fs＊i1/samples，8]；

f2＝N[fs＊i2/samles，8]；

Print[″f1＝″，fs-f1，″(aliased＝″，f1，″)，″，

″f1=″，fs-f2，″(aliased＝″，f2，″).″]；

Print[″yl[″，i，″]＝″，yl[[i]]]；

If[redo＝＝(i-1)，

redo＝-1；

]；

Print[″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊″]；

Print[″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊ Optimization ＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊″]；

search＝1；(＊related search radius＊)

For[fi＝start，fi<＝end，++fi，

i＝fi+1；(＊offset for non-zero based array indexing＊)

If[！x[[i]]；

Continue[]；

]；

Page23

Slavin.txt

cycle＝samples/freqIndex；

(＊measurable＊)

freqIndex＝fi*skip；

Print[″fi＝″，fi，″，freqIndex＝″，freqIndex，

″，cycl delay(samples)＝″，N[cycle]]；

Print[″y[″，fi，″]＝″，yl[[i]]]；

lowestResiduals＝Table[Intinity，{solutions}]；

bestResults＝Table[{}，{solutions}]；

startPoints＝Table[{}，{solutions}]；

solution＝0；

smallestResidueIndex＝1；

kk＝maxTries；

nestStart＝{}；

random＝False；

While[--kk>＝0，

slen＝Length[x[[i-1]]]；

random＝(slen＝＝0)||(kk！＝maxTries-1)||IsSame[x[[i-1]]]；

If[random，

If[slen＝＝0，

xStart＝Flatten[

{0.1＊(Random[]-0.5)，

Table[Exp[0.1＊(Random[]-0.5)]，{order}]，

Table[cycle＊(Random[]-0.5)＊0.5，{order}]}]；

，

xStart＝x[[i-1]]

Table[Exp[0.4＊(Random[]-0.5)]，{slen}]；

]；

，

xStart＝x[[i-1]]；

]；

{results，residual}＝DoSolve[e2，g，dfts，symbVars，

yl[[i]]，xStart，sx，

samples，ix，freqIndex，acc，

maxIterations]；

If[Length[results]<＝0，

Continue[]；

]；

If[！random，

Print[″result＝″，results]；

Print[″residual＝″，residual]；

x[[i]]＝results；

Break[]；

]；

(＊center the angles as much as possible＊)

xStart＝results；

parity＝0；

For[jj＝1，jj<＝order，++]]，

xx＝Floor[[results[[order+jj+1]]＊(2.0/cycle))+1/2]；

xStart[[order+jj+1]]-＝(xx＊cycle＊0.5)；

parity+＝xx；

]；

If[OddQ[parity]，

xStart[[1]]＝-xStart[[1]]；

]；

(＊reoptimize from more centered starting point＊)

{results，residual}＝DoSolve[e2，g，dfts，symbVars，

yl[[i]]，xStart，sx，

samples，ix，freqIndex，acc，

maxIterations]；

If[Length[results]<＝0，

Continue[]；

Page24

Slavin.txt

]；

sortList＝{Table[results[[jj+order+1]]，{jj，1，order}]，

Table[results[[jj+1]]，{jj，1，order}]，

Table[xStart[[jj+order+1]]，{jj，1，order}]，

Table[xStart[[jj+1]]，{jj，1,order}]}：

(＊rank solutions in order of second sub-list (angles)＊)

sorted＝Transpose[Sort[Transpose[sortList]]]；

results＝Flatten[{results[[1]]，sorted[[2]]，sorted[[1]]}]；

xStart＝Flatten[{xStart[[1]]，sorted[[4]]，sorted[[3]]}]；

(＊see if solution found already＊)

found＝False；

For[jj＝1，jj<＝solutions，++jj，

If[lowestResiduals[[jj]]＝＝Infinity，

(＊reached end of valid solutions，nothing found＊)

Break[]；

，(＊else＊)

If[Chop[lowestResiduals[[jj]]-residual]＝＝0.0，

found＝True；

Break[]；(＊solution already found＊)

]；

If[found，

Continue[]；

]；

(＊no solution found＊)

If[jj>solutions，

(＊finding another solution when no more are allowed＊)

Print[″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊″，

″Too many solutions found.″，

″＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊″]；

Abort[]；

]；

(＊new solution，so reset loop limit＊)

kk＝extraTries；

(＊find index into data set with lowest error residual＊)

If[jj>1，

If[residual<lowestResiduals[[smallestResidueIndex]]，

smallestResidueIndex＝jj；

]；

(＊find a multiplicity of local solutions at approximately Pi

intervals on all angles＊)

lowestResiduals[[jj]]＝residual；

bestResults[[jj]]＝results；

startPoints[[jj]]＝xStart；

++solution；

Print[″Solution″，solution]；

Print[″Starting vector x＝″，xStart]；

Print[″Minimum at″，results]；

Print[″Lowest residual error＝″，residual]；

++jj；

state＝Table[-search，{order}]；

While[True，

(＊find next state＊)

xStart＝results；

(＊Print[″state＝″，state]；＊)

parity＝0；

For[xx＝1，xx<＝order，++xx，

Page25

Slavin.txt

parity+＝state[[xx]]；

xStart[[order+1+xx]]+＝(cycle＊state[[xx]]＊0.5)；

]；

If[OddQ[parity]，

xStart[[1]]＝-results[[1]]；

]；

{qresults，residual}＝DoSolve[e2，g，dfts，symbVars，

y[[i]]，xStart，sx，

samples，ix，freqIndex，acc，

maxIterations]；

If[Length[qresults]>0，

If[jj>solutions，

″Too many solutions found.″，

Abort[]；

]；

lowestResiduals[[jj]]＝residual；

bestResults[[jj]]＝qresults；

startPoints[[jj]]＝xStart；

++jj；

]；

xx＝1；

While[xx<＝order，

state[[xx]]+＝1；

If[state[[xx]]<＝search，

Break[]；

]；

state[[xx]]＝-search；

++xx；

]；

If[xx>order，

Break[]；

]；

]；(＊end of While[True，]＊)

]；(＊end of While[--k>＝0，]＊)

If[random，

If[Length[bestResults[[1]]]＝＝0，

Print[″Solutions not found.″]；

(＊leave x[[i]]as True indicates measurable but not valid＊)

x[[i]]＝bestResults[[smallestResidueIndex]]；

]；

(＊

＊for each i，x[[i]] is either a legitimate list of solutions，or is

＊boolean False indicating it is not measurable due to alias frequency

＊conflicts，or is True indicating it is measurable，but that a

＊solution hasn′t been previously tound.

＊)

Print[″Filling in holes using search:″]；

(＊fill in holes where no solution was found＊)

doneSomething＝True；

While[doneSomething，

doneSomething＝FALSE；

for[fi＝start，fi<＝end，++fi，

i＝fi+1；(＊offset for non-zero based array indexing＊)

If[(Length[x[[i]]]＝＝0)&& x[[i]]]，

(＊x[[i]]is measurable＊)

freqIndex＝fi＊skip；

Page26

Slayin.txt

If[(fi>start)&&(Length[x[[i-1]]]>0)，

xStart＝x[[i-1]]；

{results，residual}＝DoSolve[e2，g，dfts，symbVars，

yl[[i]]，xStart，sx，

samples，ix，freqIndex，acc，

maxIterations]；

If[Length[results]>0，

x[[i]]＝results；

doneSomething＝True；

Print[″x[″，fi，″]＝″，x[[i]]]；

Print[″residual＝″，residual]；

]；

If[(Length[x[[i]]]＝＝0)&&

(fi<end)&&(Length[x[[i+1]]]>0)，

xStart＝x[[i+1]]；

{results，residual}＝DoSolve[e2，g，dtts，symbVars，

yl[[i]]，xStar，sx，

samples，ix，freqIndex，

acc，maxIterations]；

If[Length[results]>0，

x[[i]]＝results；

doneSomething＝True；

Print[″x[″，Fi，″]＝″x[[i]]]；

Print[″residual＝″，residual]；

]；

Print[″Forcibly filling all renaining holes using previous valuesls:″]；

For[fi＝start+1，fi<＝end，++fi，

i＝fi+1；(＊offset for non-zero based array indexing＊)

If[Lenth[x[[i]]]＝＝0，

x[[i]]＝x[[i-1]]；

]；

Print[″Scale output values:″]；

(＊scale as separate step to avoid problems with valid handling＊)

For[fi＝start，fi<＝end，++fi，

i＝fi+1；

If[scaling[[i]]＝＝0.0，

temp＝(scaling[[i-1]]+scaling[[i+1]]＊0.5；

，

temp＝scaling[[i]]；

]；

x[[i，1]]/＝(temp^order)；

(＊Print[″output:x[″，fi，″]＝”，x[[i]]]；＊)

]；

Return[x]；

]

(＊finds amplitudes and delay (in sample units)response of even length real

tap filter＊)

Response[real]taps_]:＝

Module[

{taps，dcDelay，complex，len，l1，ampls，delays，i，ampl，

radiansPerSample，delay，prevangle，angle，diffangle，newdelay，mod，

scale，xtaps}，

taps＝Length[realtaps]；

Page27

Slavin.txt

(＊Print[″realtaps＝″，realtaps]；＊)

(＊calculate the dc delay without Arg[] aliasing.＊)

dcDelay＝DCdelay[realtaps]；

(＊extend input to scale＊taps to use a longer FFT which

should reduce phase wrapping ambiguities＊)

scale＝2；

xtaps＝Flatten[{realtaps，Table[0，{(scale-1)taps}]}]；

(＊taps input returns l1＝(taps/2)+1 items＊)

complex＝Fourier[xtaps，{FourierParameters->{1，1}}]；

(＊Print[″complex＝″，Take[complex，{1，taps＊scale，scale}]]；＊)

len＝scale＊taps/2；

l1＝len+1；

ampls＝Table[0，{l1}]；

delays＝Table[0，{l1}]；

delays[[1]]＝dcDelay；

ampls[[1]]=Abs[complex[[1]]]；

mod＝Pi；

(＊for all non-dc frequencies＊)

For[i＝1，i<＝len，++i，

ampl＝Abs[complex[[i+1]]]；

ampls[[i+1]]＝ampl；

If[ampl<10^-12，

(＊if amplitude too small for accurate angle，use previous

delay as this is needed by the delay reconstruction＊)

delays[[i+1]]＝delays[[i]]，

(＊ else＊)

(＊reconstruct delay using previous delay which should

work well as long as delay changes，when translated

to angles，are not>Pi bwtween adjacent frequencies

at the new frequency＊)

(＊conversion ratio between delays and radians＊)

radiansPerSample＝N[Pi＊i/len]；

delay＝delays[[i]]；(＊get previous delay＊)

(＊convert previous delay to angle at new frequency＊)

prevangle=delay＊radiansPerSample；

(＊get angle at new frequency＊)

anale＝Mod[Arg[complex[[i+1]]]，Pi，prevangle-Pi/2]；

(＊Print[″angle＝″，angle]：＊)

newDelay＝angle/radiansPerSample；

delays[[i+1]]＝newDelay；

]；

(＊where magnitude is too small，angles are unreliable，so

average adjacent delays to get delay instead＊)

For[i＝1，i<len，++i，

If[ampls[[i+1]]<10^-12，

delays[[i+1]]=(delays[[i]]+delays[[i+2]])/2；

]；

If[ampls[[len+1]]<10^-12，

delays[[len+1]]=2 delays[[len]]-delays[[len-1]]；

]；

Return[{Take[ampls，{1，len+1，scale}]，

Take[delays，{1，len+1，scale}]}]；

]

Page28

Slavin.txt

(＊in each list，first value is DC，last is at w＝Pi＝fs/2＊)

Taps[ampls_，delays_]:＝

Module[

{l1，l2，len，taps，angles，complex，targetDelay，realtaps}，

l1＝Length[ampls]；

l2＝Length[delays]；

If[l1！＝l2,

Print[″Missmatched amplitude and delay listl lengths.″]；

Abort[]；

]；

len＝l1-1；

taps＝2＊len；

angles＝Table[delays[[i+1]]＊i＊Pi/len，{i，0，len}]；

(＊Print[″angles＝″，angles]；＊)

complex＝PolaroComplex[ampls，angles]；

(＊Print[″complex＝″，complex]；＊)

realtaps＝InverseFourier[complex，{FourierParameters->{1，1}}]；

Return[realtaps]；

]

(＊

＊Finds symbolic IM and harmonic frequency terms using the frequency

＊index and samples symbols 'ix’and'sx'.’vars’is the list of variables

＊defining the filter product to use.It is of the form；

＊ Flatten[{k，Table[a[j]，{j，1，order}]，Table[d[j]，{j，1，order}]}]；

＊)

Terms[ix_，sx_，vars_]:＝

Module[

{r0，r1，r2，r3，r4，r5，t，tt，a，b，c，j，cost，sint，len，i1，i2，p，q，nrVars，order，

q0，q1，q2，q3，q3a，q4，q5，q6，q7，q8，terms，j1，j2，t1，t2}，

nrVars＝Length[vars]；

order＝(nrVars-1)/2；

If[！InteggerQ[order]，

Print[″Vars list must be of from:″

″{kx，ax0，ax1，...，axn，dx0，dx1，...，dxn}″]；

Abort[]；

]；

q0＝TrigReduce[(Cos[(ix-1)tt]+Cos[(ix+1)tt])^order]；

(＊extract tt variable outside as a common factor for each cos term＊)

q1＝q0/.Cos[a_]:>Cos[Apart[a]]；

(＊expand into seperate terms＊)

q2＝Expand[q1]；

(＊make into a list*)

q3＝Table[q2[[j]]，{j，1，Length[q2]}]；

(＊eliminate DC frequencies from list＊)

q4＝Table[If[NumericQ[q3[[j]]]，0，q3[[j]]/(q3[[j]]/.tt->0)]，

{j，1,Length[q3]}

(＊eliminate the Cos[] arount each term＊)

q5＝q4/.Cos[a_]:>Apart[a]；

(＊form(k＊ix+m)＊tt->k＊ix+m for integer k，m＊)

q6＝q5/.a-＊tt:>a；

(＊make ix term positive:if k is-ve，then its derivative is-ve＊)

q7＝Table[If[D[q6[[j]]，ix]<0，-q6[[j]]，q6[[j]]]，{j，1，Length[q6]}]；

(＊eliminate lower-order harmonics from term list＊)

If[order>＝2，

j1＝j2＝1；

(＊a higher order contains all the terms in the lower order＊)

t1＝Sort[Expand[q7]]；

t2＝Sort[AllTerms[ix，order-2]]；

len1＝Length[t1]；

Page29

Slavin.txt

len2＝Length[t2]；

While[j2<＝len1，

If[t1[[j1]]＝＝t2[[j2]]，

t1[[j1]]＝{}；

If[j2＝＝len2，

Break[]；

]；

j2+＝1；

]；

j1+＝1；

]；

terms＝Flatten[t1]；

，

terms＝q7；

]；

len＝Length[terms]；

(＊convert into form:Cos[(k＊ix+m)＊2＊Pi＊t/sx]＊)

cost＝terms/.a_:>Cos[Apart[a＊2＊Pi＊t/sx]]；

(＊find corresponding Sin[] terms＊)

sint＝cost/.Cos[a_]:>sin[a]；

(＊expand general result for order＊)

r0＝product[Cos[2＊Pi/sx(ix-1)(t-vars[[j+1+order]])]+

vars[[j+1]]Cos[2＊Pi/sx(ix+1)(t-vars[[j+1+order]])]，

{j，1，order}]；

r1＝TrigReduce[r0]；

r2＝r1/.Cos[a_]:>Cos[Collect[Together[a]，t]]；

r3＝r2/.Cos[b_t+c_]->Cos[b＊t]＊Cos[c]-Sin[b＊t]＊Sin[c]；

r4＝r3/.{Cos[a_]:>Cos[Apart[a]]，Sin[a_]:>Sin[Apart[a]]}；

r5＝Table[{Coefficient[r4，cost[[j]]]，Coefficient[r4，sint[[j]]]}，

{j，1，len}]；

Return[{terms，vars[[1]]＊r5}]；

]

WfmColor[list_，color_]:＝

Module[

{max，r，g，b，data}，

r＝color[[1]]；

g＝color[[2]]；

b＝color[[3]]；

data＝Transpose[list][[2]]；

max＝Max[data]；

Return[ListPlot[list，{PlotJoined->True，

GridLines->Automatic，

PlotRange->{max-150，max+5}，

ImageSize->600，

Plotstyle->RGBColor[r，g，b]}]]；

]

EndPackage[]

Page30

Claims

1.一种校准***，包括：

第一模拟信号发生器，产生第一模拟信号；

第二模拟信号发生器，产生第二模拟信号；

模拟加法器，连接到所述第一模拟信号发生器和所述第二模拟信号发生器，用于将所述第一模拟信号和所述第二模拟信号相加并提供模拟输出；

ADC，具有连接到所述模拟输出的输入和数字输出；

具有连接到所述数字输出的输入的线性校正器，实现对滤波器乘积求和；

采集存储器，连接到所述数字输出；和

处理器，用于控制所述第一模拟信号发生器和所述第二模拟信号发生器，被连接以读取所述采集存储器数据，还被连接到所述线性校正器以便为所述滤波器乘积编程滤波系数；

其中所述第一模拟信号发生器产生第一正弦波，所述第二模拟信号发生器产生第二正弦波；以及

其中所述第一正弦波和所述第二正弦波具有对应于基数据记录长度的频率的互素倍数的频率。

2.如权利要求1所述的校准***，其中提供额定频率指数，所述额定频率指数是对应于所述基数据记录长度的频率的偶倍数，所述第一正弦波对应于所述额定频率指数减一，以及所述第二正弦波对应于所述额定频率指数加一。

3.如权利要求1所述的校准***，其中所述处理器包括用于寻找目标函数最小值的专用电路。

4.如权利要求1所述的校准***，其中所述处理器是运行软件程序以寻找目标函数最小值的通用处理器。