非均匀采样范文

非均匀采样范文（精选7篇）

非均匀采样 第1篇

传统的信号处理方法都是建立在数字信号均匀采样理论基础上的, 然而为提高采样速率进行多通道数据采样[1,2]时, 由于采样时间存在偏差或者各通道之间的参数不一致, 所产生的周期性非均匀采样问题将会影响后期的信号处理性能;在合成孔径雷达[3]以及合成孔径声纳[4]中, 经常采用多接收阵元技术来解决测绘速率和测绘距离之间的矛盾, 其方位向空间采样也可能会面临非均匀采样, 从而影响目标方位向聚焦性能, 使得旁瓣增大甚至出现虚假目标。因而对非均匀采样信号进行重构处理具有非常重要的意义。

本文首先介绍周期性非均匀采样信号模型, 然后推导了非均匀信号频谱和均匀采样信号频谱之间的关系, 得出非均匀采样序列频谱是均匀采样序列频谱加权和的结论, 再利用加权系数和频率无关的原理[5], 从而将非均匀信号重构问题转化为加权系数的求解问题。最后, 采用非均匀采样的正弦信号和线性调频信号对文中所研究的方法进行了验证。

1 信号模型

分析周期性非均匀采样信号最基本的方法[6]就是将其看成一个多通道数据采样问题, 如图1所示。

对于如图1所示的周期性非均匀采样序列, 其采样时刻为:

其中T为采样周期, 为一个周期为M的序列, 如果令n=k M+m, m的取值范围为[0, M-1], 那么非均匀序列的采样时刻可以变换为:

其中表示以采样周期T为参考的非均匀采样时间偏移量。

显然, 第m个通道的采样序列sm的采样率已经变成了原始采样率的1/M倍。为了恢复原始采样数据s, 将各通道采样点之间分别插入M-1个零点, 即:

将第m个通道补零后的数据延迟m个采样时间单位, 即:

其中z-1表示单位时间延迟操作。

于是, 将各通道延迟后的数据叠加后, 可得到原始采样数据, 即

2 非均匀信号频谱重构

2.1 傅氏频谱

假定被采样的原始模拟信号是一个带限信号, 即其频带范围为, 对非均匀序列分解后的各均匀信号进行傅里叶变换, 即:

原始模拟信号的谱以及泊松公式, 将上式变换为:

将各通道信号的谱叠加后即可得到视为均匀采样信号的非均匀信号的傅氏谱, 即

其中。由此可见, 非均匀信号的傅氏谱为原模拟信号的频谱加权和, 但是由于这个权值与频率有关, 所以较难重构成均匀采样信号的谱。

2.2 数字频谱

根据傅里叶变换的原理, 直接对非均匀采样后的信号进行傅氏变换, 可得

与推导傅氏谱方法类似, 经过简单代数变换之后, 可得

观察式 (14) , 不难发现, 非均匀信号的谱仍然是原始模拟信号的谱的加权和, 并且此时的权值已经与频率无关, 从而提供了一种信号重构方法。

2.3 非均匀信号重构

于是, 可以得到频带范围内频点为处的谱值。

根据式 (15) 和式 (16) , 可到频带范围内频点为处的谱值。

将式 (17) 写成矩阵形式, 可到

其中为一个M1的矩阵

其中上标T表示转置。

A表示一个MM的矩阵, 其表达式为:

3 仿真实验

本节将采用仿真实验来验证本文方法的有效性, 首先采用一个频率为1Hz的实正弦信号, 采样周期T=0.1 s, 采样通道数M=4, 4个通道的采样时间偏差分别为0 s, 0.0033 s, 0.0066s, 0.0099 s。为了对比分析, 图2 (a) 先给出了理想情况下均匀采样序列的频谱。图2 (b) 为将非均匀信号视为均匀信号时所获得的频谱, 图2 (c) 为按照傅里叶变换的方式, 直接对非均匀信号进行频谱变换的数字谱, 图2 (d) 为按照2.3节方法重构后的谱。

当采用调频率为0.3418的复基带线性调频信号时, 采样周期为T=0.05 s, 采样通道数仍为4, 4个通道的采样时间偏差分别为0 s, 0.1429 s, 0.2857 s, 0.4286 s。为了对比分析, 图3 (a) 先给出了理想情况下均匀采样序列的频谱。图3 (b) 为将非均匀信号视为均匀信号时所获得的频谱, 图3 (c) 为按照傅里叶变换的方式, 直接对非均匀信号进行频谱变换的数字谱, 图3 (d) 为按照2.3节方法重构后的谱。

观察图2和图3, 可以发现如果直接利用非均匀信号的频谱进行处理, 将会严重影响后续的信号处理性能, 而采用本文方法可以精确重构原始信号的频谱, 从而大大降低实际中多通道采样的要求。

4 结语

非均匀信号将会严重降低信号处理性能, 本文基于此推导了非均匀信号和其等价的多通道均匀采样频谱之间的关系式, 从而将信号重构转化为求解与频率无关的权值问题。并采用正弦信号和线性调频信号验证了本文方法的有效性。

参考文献

[1]Papoulis A.Generalized sampling expansion.IEEE Transactions on Circuits and system.1977, CAS-24 (11) :652-654

[2]Brown J.Multi-channel sampling of low-pass signals.IEEE Transactions on Circuits and system.1981, CAS-28 (2) :101-106

[3]Gebert N, Krieger G, Moreira A.SAR signal reconstruction from non-uniform displaced phase center sampling in the presence of perturbation[J].IEEE International Geoscience and Remote Sensing Symposium, 2005 (IGARSS’05) , Seoul, Korea, 2005, 2:1034-1037.

[4]杨海亮, 唐劲松, 张森.基于非停走停及方位非均匀采样的多接收阵合成孔径声纳CS成像算法[J].高技术通讯, 2009, 19 (9) :939-945

[5]Jenq Y C.perfect construction of digital spectrum from nonuniformly sampled signals.IEEE Transactions on Instrumentation and measurement, 1997, 46 (3) :649-652

非均匀采样 第2篇

人眼启发的图像非均匀映射变换模型构造算法

人眼的视觉信息采集是非均匀的.这种特征和注意力机制确保了人类视觉系统在信息处理方面的优先性.首先介绍了对数极坐标映射,然后基于人类视觉系统的`构造,提出了非均匀映射模型的构造规则.基于此规则,以线性模型为例进行了应用性研究.扩展了非均匀映射模型的范围,促进了空间变分辨率视觉理论的进一步发展.

作 者：王琪 李言俊 张科 WANG Qi LI Yan-jun ZHANG Ke 作者单位：西北工业大学航天学院,西安,710072刊 名：宇航学报 ISTIC PKU英文刊名：JOURNAL OF ASTRONAUTICS年，卷(期)：200627(6)分类号：V557关键词：非均匀映射 对数极坐标变换 模型构造 Non-uniformed mapping Log polar transform (LPT) Model construction

非均匀采样 第3篇

关键词： 光纤光学； 光纤陀螺； 法拉第效应； 球面磁场； 磁敏感性

中图分类号： TN 253文献标志码： Adoi： 10.3969/j.issn.1005-5630.2014.05.012

引言

光纤陀螺是一种新型的惯性测量器件，有着传统陀螺难以比拟的优点。作为一种全固态陀螺，光纤陀螺具有动态范围大、响应速度快、体积小、抗冲击振动、启动时间短、工艺简单以及易于大批量生产等优点[1-2]。

光纤陀螺自1976年问世以来，经过30多年的发展，其在军用和民用领域取得了巨大的成就，被公认为取代机械陀螺的下一代陀螺。光纤陀螺磁敏感性误差作为光纤陀螺的主要非互易误差源之一[3]，是评价光纤陀螺性能的重要参数。本文通过分析光纤陀螺内部电路辐射磁场的分布特征，建立了球面非均匀磁场中的光纤陀螺磁敏感误差模型。探讨了球面非均匀磁场对光纤陀螺磁敏感性误差的影响。

1光纤陀螺磁敏感性机理

磁光法拉第效应是当线偏振光通过处于磁场作用下的透明介质时，其线偏振光的偏振角会发生旋转，产生磁场作用下的一种旋光现象[4-5]。由于磁光法拉第效应，在单模光纤中磁场改变了构成入射线偏振光的左、右圆偏振光的相位，导致两束反向传播的线偏振光的偏振面产生一个夹角，使光在光纤环中传输时产生一个非互易相位差[6]。由于这一误差无法与光纤陀螺的Sagnac效应区分，因此产生法拉第效应误差，导致光纤陀螺具有磁敏感性。

2光纤陀螺内部电路辐射磁场

光纤陀螺的内部磁场主要由陀螺内部电路板产生。利用电磁场分布扫描系统，对光纤陀螺的电路板进行了电磁场辐射特性扫描。

2.1测试设备

利用瑞典Detectus公司生产的RX644EH型电磁场分布扫描仪，Agilent公司的E4440A频谱分析仪，以及计算机控制与显示软件构成的电磁场分布扫描系统，对光纤陀螺的电路板进行电磁场辐射特性扫描。图1为电磁场分布扫描系统框图。光学仪器第36卷

针对测试对象，XYZ轴的移动距离（被测设备最大尺寸）为600×400×200 mm，最小移动步径为1 mm，定位精度±0.3 mm。扫描测试时，选择低频磁场探头LFB-3进行测试，扫描频带为0～50 MHz。

2.2扫描结果与分析

陀螺主控电路正面辐射强度最大的8.761 MHz的辐射特性如图2所示。从图2可以看出，随着探头与主控电路板上表面之间距离的增加，该电路板辐射强度逐渐下降。这说明电路板辐射磁场分布更接近球面分布，而非强度与探头高度无关的柱面分布的匀强磁场。

3球面非均匀磁场中的光纤陀螺磁敏感误差模型

若假定光纤环置于平行于光纤环平面的磁感应强度为B0的均匀磁场中，θ为磁场方向与光纤环所成的角度，θ0为磁场方向相对基准轴的角度，τ（θ）为光纤扭转率，V为费尔德常数，Δβ为光纤双折射率，r为光纤环半径，m为光纤环在横向上缠绕的层数，则该磁场所产生的径向法拉第相位误差可以表示为[7]：光纤陀螺磁敏感误差与磁场源的位置、大小，以及光纤环的高度、半径有关。磁场源距离光纤环越近，光纤陀螺的磁敏感相位误差越大。2） 在光纤环外部，光纤陀螺磁敏感误差随着球面磁场距离光纤环上表面高度H0的增大而减小，随着球面磁场距离光纤环侧边缘距离R的增大而减小。3） 当磁场源位于光纤环中心轴附近时，光纤陀螺磁敏感误差较小。实际上，由于扭转量不是均匀分布的，各点产生的法拉第相位效应也不会完全抵消。这个结论说明，将光纤陀螺内部辐射较大的元件置于光纤环中心轴附近，可减小其对陀螺输出的影响。4） 球面磁场源在光纤环内部中心轴向移动时，光纤陀螺磁敏感误差基本不变。

5实验结果

通过设置光纤陀螺内部电路距离光纤环上表面的高度，从而改变球面磁场源距离光纤环上表面的高度，采集陀螺的输出数据，如图7所示。改变球面磁场源距离光纤环侧边缘的水平距离，采集陀螺的输出数据，如图8所示。改变球面磁场源处于光纤环内部的高度，采集陀螺的输出数据，如图9所示。图7表明，在光纤环外部，光纤陀螺磁敏感误差与球面磁场距离光纤环上表面高度H0成反比。图8表明，光纤陀螺磁敏感误差与球面磁场距离光纤环侧边缘距离R也成反比。图9表明，光纤陀螺的光纤环对于在其内部中心轴向移动的球面磁场源而言，磁敏感性强度基本不变。分析3组实验数据，发现图7、图9中的陀螺输出数据较小，说明当磁场源位于光纤环中心轴附近时，光纤陀螺磁敏感误差较小。实验结果验证了数值模拟分析的正确性，同时也表明了所建立的球面非均匀磁场中的光纤陀螺磁敏感误差模型的合理性。

6结论

基于光纤陀螺内部电路辐射磁场的扫描结果及其磁敏感性机理，着重研究了光纤陀螺在球面非均匀磁场中的磁敏感性特征。1）利用电磁场分布扫描系统，得到光纤陀螺内部电路辐射磁场为呈球面分布的非均匀磁场；2）在建立了光纤陀螺在球面非均匀磁场中的磁敏感误差模型的基础上，对处于不同球面非均匀磁场下的光纤环法拉第效应进行了数值模拟分析；3）实验验证了数值模拟分析的正确性。得到球面磁场源距离光纤环越近，光纤陀螺磁敏感误差越大；球面磁场源位于光纤环中心轴附近时，光纤陀螺磁敏感误差较小；球面磁场源在光纤环内部中心轴向移动时，光纤陀螺磁敏感误差基本不变。基于以上结论，为了减小光纤陀螺自身内部电路辐射磁场对光纤环的影响，在设计陀螺内部结构时，应考虑将内部电路的主要辐射磁场源集中于光纤环中心轴附近。全面分析了光纤陀螺在球面非均匀磁场中的磁敏感性特征，为以后光纤陀螺内部结构的优化设计和光纤陀螺磁敏感性抑制方法的研究提供了一定的参考。

参考文献：

[1]UDD E，LEFEVRE H C，HOTATE K.Fiber optic gyroscope：20th anniversary conference[J].SPIE，1996，2837：1.

[2]谭曦，刘军，殷建玲.正方形亥姆霍兹线圈的磁场均匀性[J].光学仪器，2012，34（1）：39-44.

[3]HOTATE K，TABE K.Drift of an optical fiber gyroscope caused by the Faraday effect：e xperiment[J].Lightwave Technology，1987，5（7）：997-100.

[4]王夏霄，宋凝芳，张春熹，等.光纤陀螺磁敏感性的实验研究[J].北京航空航天大学学报，2005，31（10）：1116-1120.

[5]董宇，张悦，华文深，等.基于磁致旋光效应的光电装备隐身技术[J].光学仪器，2012，34（6）：80-85.

[6]李坚，宁提纲.干涉型光纤陀螺中磁光Faraday效应的研究[J].光电子·激光，2007，18（4）：400-403.

基于非均匀采样重构的图像插值算法 第4篇

图像插值是一种常用的图像缩放方法，尤其对于运算复杂度有限制的消费类电子，插值是最通常的方法。邻近插值算法所消耗的运算和内存需求是所有插值算法中最小的，但是其视觉效果很差，存在马赛克效应。不同的插值算法区分的方法是其所使用的插值基函数。传统的插值方法所使用的基函数是基于多项式。在数学上，通常使用多项式对连续曲线进行逼近。另一类插值方法是基于加窗的sinc函数，现有文献已经讨论过不同的窗函数。在托马斯的论文中[1]，介绍了加Blackman-Harris窗的辛克函数(sinc)以及二阶连续的立方卷积函数是两个在医疗图像插值中最好的两个插值基函数。其他较为优秀的基函数有Lanczos[2]，它被MatlabTM中的图像处理工具包所使用。

本文所讨论的插值算法是基于Kim和Bose的非均匀采样重构算法[3]。他们的算法的提出解决了从信号的非均匀采样点重构原始信号的问题。他们的算法对于原始信号有条件限制，信号必须为带限和周期的。尤其当周期性条件不符合时，重构出的信号是病态的。然而，当选择偏移栅格点一段距离的均匀采样点时，此算法可以以较高的精度重构出位于栅格点的信号。

因此，从此算法在均匀采样点上的应用可以推导得到一个新的插值基函数。同其他的基函数比较，本文的基函数在视觉和峰值信噪比指标上均表现出最好的效果。

1 重构算法

1.1 原始重构算法

在Kim和Bose的算法中，信号能够从其非均匀采样中重构。重构算法可用图1表述。图1中的曲线是原始信号y=32sin(t)+16cos(0.5t-0.11)+22sin(0.25t)。在重构中的参数周期T为8π，采样点为32个，这两个参数满足了周期和带限的前提条件。在图1中的三角标注了原始信号的非均匀采样点，而十字则标注了重构点。图1说明原始的重构算法能够完全重构出原始信号，只要原始信号满足周期和带限条件。在图1中发现，只有一半的原始信号被采样了，但是算法仍能够无错地重构出整个原始信号。

原始重构算法可以通过如下公式描述，

公式(1)可以简化为

其中，(i,j)=exp(jtj(-N/2-1+i)wx)，i,j=1,2，，N,f向量为非均匀采样得到的信号，F向量为重构位置信号的离散傅里叶系数，F还不是标准离散傅里叶系数(F(0)，F(1)，F(N-1))，傅里叶系数存在这样一个关系(F(a)=F(a+kN)，k可以为任何整数)，通过调整F向量元素的顺序可以得到标准化的离散傅里叶系数向量F'。通过对F'做傅里叶反变换运算就能够得到重构栅格位置的原始信号值fr，这个过程可用如下公式表述

图1已经展示了Kim的算法可以用以从非均匀采样点重构出周期和带限的原始信号。然而原始算法的这两个条件中的周期性条件大大限制了它在实际中的应用，因为实际中的图像绝大多数不满足周期性条件。图2可以说明这个问题，将算法中的周期参数T设为8π+1，从而同原始信号的周期8π不匹配。在这种情况下，原始信号就不能通过非均匀信号完美地重构。重构得到的信号同原始信号误差非常大，因此无法将Kim的算法应用到实际图像中。

在研究中发现，如果采样点位置均匀地分布在整个信号区间内，那么此算法重构出的信号和原始信号间的误差非常小，这个结果使我们考虑将此算法应用于图像插值中。在之后的章节中本文将说明改进之后的重构算法将比现有的图像插值算法的性能更好。

1.2 改进的重构算法

首先建立原始采样点和重构点之间的数学关系。这种数学表述能够帮助我们了解并克服原算法对于信号的周期性要求，同时能够简化运算。公式(1)可以写成

其中，f为原始采样，F为N点的非正常序DFT变换系数。为了得到重构位置的值，必须先将F重构成常序DFT系数列F'。在式(2)中ti可以通过下式表述

其中，ti代表了第i个采样点。Δt则代表了采样点和对应的栅格点之间的偏移，这种位置关系可以通过图3表述。

上述的整个过程可以通过式(6)描述

其中，

式(6)中的矩阵A,B和C分别为

从式(6)可以看出向量f'是Δt的函数。而Δt则描述了重构点和采样点之间的偏移。式(6)可以写成

矩阵G中的元素gnk为

通过式(7)可以建立起采样点和重构点间的数学关系。

向量f'中的每个元素都可以写成如下表达式

从式(9)中可以看出每一个重构信号均受到每个原始采样点的影响，因此每个采样点在重构过程中存在一个响应函数，这个函数可以用式(8)表述，每个采样点的响应是式(8)偏移各自响应的位移。而gnk(Δt)是一个周期函数，因此原算法限制了信号必须为周期信号。图4给出了当采样4个点，即N=4的情况下的各个采样点的单位响应曲线。图中的响应曲线的阴影部分表明了每个采样点的周期性影响，即原算法不能应用于实际图像的区域。在图4中最后一个响应合成图中可以看到1.0～2.0的区间内没有阴影部分，这说明此处重构得到的信号不受周期条件约束，能够重构带限信号。当带限条件不满足时，重构得到的信号也接近原始信号。因此，此区间内的重构信号可以应用于图像插值中。此区间对应向量f'中的元素为f'(N/2)。f'(N/2)可以通过式(10)描述:

2 插值基函数

传统的图像插值算法可以写成以下的形式

在上式中，u(x)代表了插值基函数，而ck则是原始图像值。本文提出的算法中，重构点的计算方法如下:

通过同式(11)比较可发现式(12)中的aN/2,k(x)和(11)式中的u(x-xk)作用相同，而aN/2,k(x)=aN/2,0(x-k)。因此aN/2,k(x)是一个新的插值基函数，它可以用下式描述

其中，N代表了插值的阶。

同其他插值基函数，如双线性、双立方、以及Lanczos等比较，本文提出的基函数(13)在频率响应、峰值信噪比和视觉感观上更优越。图像插值是一种从离散的采样点得到连续信号的方法，因此对于能够恢复的连续信号有一定的要求。奈奎斯特采样定理限定了所能够从离散信号恢复的连续信号的带宽约束。插值对于离散信号而言是一个低通滤波器。如果一个信号是带限的，那么它能够由一个理想低通滤波器完全由离散采样恢复出来。然而理想低通滤波器由于其时域上的无限脉冲响应而无法实现。所有的插值及函数都是对于理想低通滤波器的逼近，因此模糊很混淆将无可避免地引入到插值过程中。

图5给出了本文提出的基函数和Lanczos基函数的频率响应曲线，从图中可以看到在相同的基函数长度下，本文提出的基函数的截止带上的频率响应斜率更高，更接近理想低通滤波器。

3 实验结果

实验中，使用峰值信噪比指标来分析不同长度不同类型的基函数的性能。选择5张标准图片作为实验样本。实验过程是先将每幅图片用各个基函数缩小2倍，然后再将缩小后的图像用相同的基函数放大2倍。经过处理的图像同原始图像一起计算出每个基函数对应的峰值信噪比。实验结果列在表1中。从表1可以看出，基函数的长度越大则峰值信噪比性能更好。原始图像变化越圆滑，峰值信噪比越高，其包含的高频信息越少，在插值过程中的损失越少。

4 结束语

本文提出的插值基函数是从非均匀采样重构算法推导得到的。这个基函数基于余弦函数和，同其他基函数不同。给出了插值基函数的频率响应曲线，从曲线可以看出本文的基函数更接近理想低通滤波器，因而其插值性能更好。而在标准图像的实验所得到的峰值信噪比表1更加肯定了这个结论，在相同长度的基函数中，本文的基函数在每幅图像下均得到了最高的峰值信噪比。

摘要：在大规模集成电路实现图像缩放时,通常使用的算法是插值。基于信号的非均匀采样重构算法,提出了一种新的图像插值算法[1]。该插值算法所使用的基函数是基于余弦和函数。本算法同其他算法相比,在插值基函数长度相同的情况下,本算法的性能最好。相较其他算法,本算法保持了图像中更多的高频分量,同时减少了插值过程中的图像混淆效应。同时并没有增加插值所需的运算资源。

关键词：图像,插值,非均匀采样重构

参考文献

[1]Th_evenaz P,Blu T,Unser M.Interpolation revisited[J].IEEETrans.Med.Imag.,2000,19:739-758.

[2]Claude E Duchon.Lanczos Filtering in One and Two Dimensions[J].Applied Meteorology,1979,18:1016-1022.

[3]Kim S P,Bose N K.Reconstruction of2-D bandlimited discrete sig-nals from nonuniform samples[J].Proc.Inst.Elec.Eng.,June1990,137:197-204.

[4]Th_evenaz P,Blu T,Unser M.Interpolation revisited[J].IEEETrans.Med.Imag.,2000,19:739-758.

[5]Hou H S,H C Andrews.Cubic splines for image interpolation anddigital filtering[J].IEEE Trans.Acoust.,Speech,Signal Processing,1978,ASSP-26(6):508-517.

[6]HarrisFJ.Ontheuseof windowsfor harmonicanalysiswiththediscreteFourier-transform[J].Proc.IEEE,1978,66:51-83.

[7]Lehmann T M,Gonner C,K Spitzer Survey:Interpolation methodsin medical image processing[J].IEEE Trans.Med.Imag.,1999,18:1049-1075.

[8]Keys R G.Cubic convolution interpolation for digital image proce-ssing[J].IEEE Trans.Acoust.,Speech,Signal Pro-cessing,1981,ASSP-29(6):1153-1160.

非均匀采样 第5篇

SD-OCT系统中, 在波数域测量值通常经过非均匀采样得到。这主要受测量装置受到波长的影响。另一方面, 在图像重建过程中, 采用相对于波数k, k = 2π/λ 的傅里叶变换实现。因此测量数据实际上是波数域内的非均匀采样值。通常在对SDOCT数据进行成像之前需要采用插值的方法进行重采样以将其变为均匀采样的数据[1，2]。然而, 这种方法面临极大的挑战, 比如增加了系统对深度方向数据变化的敏感性, 增加了系统的局部SNR, 抑制虚假目标的旁瓣等[3，4]。

重采样处理的一种替代方法就是采用近来提出的用于SDOCT信号重建的非均匀离散傅里叶 ( NUDFT) 变换。NUDFT是DFT的一种改进, 主要用于处理非均匀采用的数据[5，6]。文献[4]表明采用NUDFT进行SD-OCT数据处理能够降低其对深度方向数据变化的敏感性。然而这是以牺牲计算效率为代价的。

本文提出了一种逆成像方法[7], 将SD-OCT系统建模为一系列线性方程组, 然后通过求解一个逆问题直接求解样本的互相关而不是单独的对每一个样本求解。由于逆问题的病态特性, 需要对其附加限制条件以得到唯一解。本文采用全变差 ( TV) 作为限制条件, 从而能够保留图像的边缘信息和降低背景噪声的影响[8 - 10]。

1 SD-OCT系统建模

图1 所示为SD-OCT的系统结构示意图。一个低相关的SLD被用作光源。大多数光源通过couper’均匀地分裂到参考支路和采样支路。然后由参考支路和采样支路的反射镜反射回的光线再次通过couper’汇聚到检测支路以实现信号的检测。在检测支路, 线性扫描的CCD器件被用来收集样本数据。在检测器曝光器件, CCD的每一像素采集对应于一个信号源波长的样本数据。由CCD的所有像素可以得到沿某一特定方向的数据。为了得到所有数据, 需要采用标号为GM1 和GM2 的透镜扫描样本的过度平面。

假设为系统测量值, G ( k) 为光源能量谱, 则一个A射线的SD-OCT系统测量值可以表示为:

其中, pr为参考点反射镜的反射率, ps ( z) 为采样支路随深度而变化的反射率。式 ( 1) 包含了三个部分。第一部分与参考支路的反射光相对应; 第二部分对应于参考支路和采样支路的互扰; 第三部分对应于不同反射层的反射光线的互扰。在所有三部分中, 第二部分最重要, 由第二部分可以重建采样信号ps ( z) 。

通常假设ps ( z) 关于z = 0 对称, 因此ps ( z) 可以由测量值通过FFT估计得到。 另一方面, 注意到, 如果信号源谱G ( k) 能够独立估计得到, 就能从式 ( 1) 中去除G ( k) pr2, 因为pr= 1, 并且将最后一部分表示为误差e ( k) 。这样, SD-OCT系统的测量值简化为:

在离散域, 式 ( 2) 可以表示为:

其中, 将轴向数据分割为N段, m ∈ { 0, 1, , M - 1} 表示第m个像素。N的值不需要和M相同, 但是N的最大值需要受到最大可检测深度的限制。将式 ( 3) 表示为矩阵形式, 可得到:

其中Hm，n= 2G ( km) cos ( 2kmzn) 。

为了更简洁地表示式 ( 4) , 使用y = [I ( k0) , I ( k1) , , I ( kM -1) ]T, x = [ps ( z0) , ps ( z1) , ps ( zN -1) ]T和e = [e ( k0) , e ( k1) , , e ( kM -1) ]T分别表示SD-OCT测量向量, A-line采样信号和误差向量。矩阵H定义为SD-OCT系统的冲击响应矩阵。则式 ( 4) 可以表示为:

则L条A-line的SD-OCT测量向量可以表示为:

2 基于TV正则化的信号重建

从式 ( 6) 可以看出, 本文任务是由Y和H重建信号X 。这是一个逆问题, 其求解可以通过式 ( 7) 实现:

其中, α > 0 为正则化参数, ‖X‖TV为X的全变差, 定义为:

为了求解式 ( 7) , 本文采用扩展的拉格朗日方法, 在式 ( 7) 中加入中间矩阵U , 则变为:

其中, α1> 0, α2> 0 为预先确定的正惩罚参数。式 ( 9) 求解可用交替最小化方法实现, 每一次迭代中包含如下两个步骤:

共轭梯度方法和Chambolle投影方法可以分别用来求解X ( i) 和U ( ) i [11，12]。矩阵U的引入并不会改变式 ( 7) 所示的最小化问题, 但是会加快收敛速度。除此之外, 采用X ( i) 和U ( i) 之间的不同作为停止迭代的条件。参数 α1> 0, α2> 0 用来平衡测量误差和TV正则化分量。当测量误差较大时, α1, α2的值较大。然而, 如果 α1, α2的值设置得过大, 会使得重建的图像严重平滑化。

3 仿真验证

SD-OCT实验系统装备有SLD源, 中心波长850 nm, 带宽50nm, 具有Gaussian功率谱。系统的分辨率为6. 38 μm, 最大观测深度为2. 42 mm。参考臂的反射镜直径为2 英寸。系统测量数据由通常使用的分光计测量得到。如图1 所示, 照相机和反射镜由计算机控制。

在第一个实验中, 测量对象为内部具有气孔的的金属盘, 测量数据包含800 条A-line, 且每一条A-line包含2 048 个样本。图2 ( a) 所示为原始测量数据。为了重建样本信号, 首先需要从原始数据中估计出源功率谱G ( km) 。注意到, 原始测量数据与高频分量调制的具有Gaussian形状的低频信号相似。低频信号来源于SLD源功率谱G ( km) , 而高频分量为式 ( 1) 中第二和第三部分表示的干涉信号以及观测器噪声。

为了减轻测量噪声, 首先将800 条A-line的数据平均。由于高斯函数的Fourier变换仍然是高斯的, 因此将低通滤波器的截止频率选为平均Fourier变换的第一队零点。然后由滤波后数据的逆傅里叶变换可以估计出源功率谱Gest ( km) 。图2 ( b) 和 ( c) 分别表示平均傅里叶变换和估计的源功率谱。将估计得到的源功率谱从原始测量数据中去除, 得到调整后的A-line测量数据I ( km) , 如图2 ( d) 所示。对每一条A-line, I ( km) 可以认为是式 ( 3) 所示的ps ( zn) 的线性组合。

进一步注意到由于G ( km) 是高斯分布的, I ( km) 的包络具有高斯分布的特性。然而其包络并不是关于零点对称的, 这似乎与H中的余弦函数相悖。事实上这是由几个原因引起的。主要原因是不同采样层ps ( zn) 之间的干涉造成的, 如式 ( 1) 中的第三部分所示。可以将其重写为:

其中, ps ( z) 为实值。这个非负分量将测量信号的包络向正极移动。其他造成不对称的原因包括系统噪声、量化噪声以及系统的非线性响应等。所有这些因素都会对重建性能造成影响。

在这点上, 比较三种不同的重建采样信号的方法: FFT、NUDET以及基于TV正则化的方法。设置参数 α1= 20, α2= 100 , 总的迭代步数为80。图3 所示为采用不同方法的重建结果, 水平方向和垂直方向分别表示A-line和深度。

图 3 不同方法的重建结果 ( a) FFT; ( b) NUDET; ( c) TV 正则化

可以看出, 与FFT方法相比, NUDET方法的重建质量更高, 尤其是在图中虚线所示具有较高观测深度的区域。然而, 采用TV正则化方法所得到的重建图像很好地保留了采样信号, 所包含的残留噪声更少, 比其他两种方法所得图像具有更好的SNR。这里的SNR表示期望信号与其他所有重建误差的比值。为了能够直观地观测各种方法的重建结果, 本文将第150 条A-line的信号用图像显示出来, 如图4 ( a) 所示; 图4 ( b) 所示的深度为n = 425 。TV正则化方法所得结果具有较好的平滑性, 这是因为该方法更好地抑制了信号的波动。

图 4 采用不同方法的重建信号 ( A-line) ( a) 轴向 ( b) 横向

在第二个实验中验证了TV正则化方法对SD-OCT系统的另一个特性深度灵敏性的作用。所使用的数据为一个银盘的横断面数据。设置参数 α1= 20, α2= 100 , 总的迭代步数为100。图5 所示为采用不同方法的重建结果以及对应区域的放大结果。

图5不同方法的重建结果 (a) FFT (c) NUDET (e) TV, (b) (d) (f) 分别为对应区域的放大观测结果

可以看出, 随着深度zn的增加, 斜线向左的趋势就更加明显, 在重建图像中就会出现图中箭头标注出的误差。这是因为波数k与波长 λ 成反比, 因此在k较大的区域, 均匀分布的 λ 就会对应较少的测量值, 在较大深度处就会造成较大的误差。在NUDET和TV方法中, 这个问题不是很明显。比较图5 ( c) 和 ( e) , 可以看出采用TV正则化方法所得图像的残留噪声较少, 与第一个实验所得结果类似。

为了研究系统的深度敏感性, 采用图5 ( b) 、 ( d) 、 ( f) , 且在斜线处于不同深度的位置选择3 条A-line。选择的A-line深度分别为l = 9, l = 601, l = 930。在这些图像中所选择的区域都远离斜线所在的位置, 比如对l = 9, 选择的区域为660 < n < 760;对l = 601, 选择的区域为460 < n < 560; 对l = 930, 选择的区域为350 < n < 450 。图6 所示为不同方法的3 条A-line。

图 6 不同深度处的重建信号 ( a) FFT ( b) NUDET ( c) TV

在图6 ( c) 中, 在3 条A-line所示的不同深度处均能观测到两个峰值, 表明具有两层。而在FFT和NUDET方法中, 在深度较深时就很难观测到两个峰值。本文对图6 ( c) 中的峰值进行定位。表1 所示为采用三种不同方法得到的峰值处对应的ps ( zn) 值。其中, n = 384, 497 和704 对应第一层, 而其他的对应第二层。正如期望得到的结果一样, 当n增加或者深度扩大时, ps ( zn) 减小。如在第一层, n从384 增加到704, 在FFT方法中归一化ps ( zn) 值从- 5. 60 d B下降到- 21. 68 d B; 在NUDET方法中, 从0 d B下降到- 16. 23 d B; 在TV方法中, 从- 1. 68 d B下降到- 12. 91 d B。然而, 在三种不同方法中, ps ( zn) 的最大值与最小值之间的差异不同。

4 结语

文中采用一种线性系统来刻画SD-OCT测量样本, 通过求解一个逆问题来实现从非均匀采样数据重建图像。采用TV正则化方法实现问题的求解。通过仿真验证可看出采用该方法能够有效地抑制噪声, 且能够得到比传统的基于FFT的方法和NUDET方法更好的重建质量。此外, 采用TV正则化方法能够降低系统的深度敏感性, 调整模型参数, 这种逆成像方法能够用于其他的OCT系统。

摘要：在频域光学相干层析成像 (SD-OCT) 中, 数据通常通过波数域的非均匀采样得到。在进行快速傅里叶变换得到图像之前需要进行一个数据重采样过程。数据重采样会造成额外的计算负担, 而且会带来一定的数据误差。采用一种逆成像算法得到SDOCT图像, 将SD-OCT系统建模为一系列线性方程组, 通过求解一个逆问题以实现SD-OCT系统的成像。利用全变差 (TV) 作为限制条件以保留图像的边缘信息, 然后由SD-OCT测量数据直接估计样本的二维互相关。仿真结果表明, 和传统方法相比, 该算法所得到的噪声残余量更低。同时还验证了利用TV限制条件以抑制SD-OCT中对衰落的敏感性的可能性。

关键词：SD-OCT,非均匀采样,全变差,重采样

参考文献

[1]Brezinski M.Optical Coherence Tomography:Principles and Applications[M].[S.l.]:Elsevier, 2006.

[2]Watanabe Y, Maeno S, Aoshima K, et al.Real-time processing for fullrange Fourierdomain optical-coherence tomography with zero-filling interpolation using multiple graphic processing units[J].Appl.Opt, 2010, 49:4756-4762.

[3]Zhang K, Kang J U.Graphic processing unit accelerated non-uniform fast Fourier transform for ultrahigh speed, real-time Fourier-domain OCT[J].Opt.Express, 2010, 18:23472-23487.

[4]Wang K, Ding Z, Wu T, et al.Development of a non-uniform discrete Fourier transform based high speed spectral domain optical coherence tomography system[J].Opt.Express, 2009, 17:12121-12131.

[5]Chan H K, Tang S.High-speed spectral domain optical coherence tomography using non-uniform fast Fourier transform[J].Biomed.Opt.Express, 2010, 1:1309-1319.

[6]Jeon M, Kim J, Jung U, et al.Boppart.Full-range k-domain linearization in spectral domain optical coherence tomography[J].Appl.Opt, 2011, 50:1158-1162.

[7]Lam E Y, Zhang X, Vo H, et al.Three-dimensional microscopy and sectional image reconstruction using optical scanning holography[J].Appl.Opt, 2009, 48, H113-H119.

[8]Xu Z, Lam E Y.Image reconstruction using spectroscopic and hyperspectral information for compressive terahertz imaging[J].J.Opt.Soc.Am.A, 2010, 27:1638-1646.

[9]Zhang X, Lam E Y.Edge-preserving sectional image reconstruction in optical scanning holography[J].J.Opt.Soc.Am.A, 2010, 27:1630-1637.

[10]Ke J, Zhu R, Lam E Y.Image reconstruction from nonuniformly-spaced samples in Fourier domain optical coherence tomography[C]//Proc.SPIE, 2012.

[11]Huang Y, Ng M K, Wen Y.A fast total variation minimization method for image restoration[J].Multiscale Modeling&Simulat, 2008, 7:774-795.

非均匀采样 第6篇

电网谐波污染对电能质量以及电网的运行安全造成了严重影响。准确地对谐波进行分析和测量是治理谐波污染的基础和前提。从设计集成电路对算法进行硬件实现的角度来看,快速傅里叶变换(FFT)由于算法的运算量较小、运算方式较为规则、可以进行同址运算等特点,具有较小的硬件实现代价,是理想的谐波分析算法,得到了广泛应用。

传统的FFT谐波分析方法建立在均匀采样的基础上,由于采样电路的限制不能实现理想的同步采样,并导致频谱泄露等问题。其解决方案一是对变换过程进行改进或者校正,使频域输出结果接近真实结果[1];二是对采样得到的数据进行同步化处理,使处理后的采样数据近似同步采样得到的数据[2,3]。这些方法中,有些运算量较大,有些电路结构过于复杂。

目前,对非均匀采样数据的处理可以分为2种类型。第1种类型是利用非均匀采样的采样时间信息,采用数值计算离散时间傅里叶变换或者改进的离散傅里叶变换获得频域信息[3,4,5];这种类型的方法具有高混叠频率特性,且由于是对连续频域波形的近似,还具有较高的频域分辨率,但是这些方法不能使用FFT进行计算,运算量很大。第2种类型是把非均匀采样的数据视为均匀采样的数据,通过去除频谱中由于非均匀采样带来的噪声重建频谱或者时域波形[6,7];这种类型的方法主要利用了非均匀采样的高混叠频率特性,由于需要额外的运算来消除采样噪声,运算量也较大。

本文提出了一种基于非均匀同步过采样技术的谐波测量算法。根据频率测量单元输出的结果产生非均匀的同步过采样时钟,通过合理设计过采样和降采样2个阶段非均匀采样时钟频率的概率分布函数,使过采样阶段的采样噪声处于降采样结果的频带之外,并使降采样阶段的采样噪声幅度小于模数转换(ADC)的量化噪声幅度。降采样得到的数据可以视为均匀采样的结果直接采用FFT进行分析,不用进行额外的消除非均匀采样噪声或者频谱泄露的运算。

采样时钟的非均匀特性大幅简化了延时锁定环路(DLL)时钟产生电路的设计。由于同步采样,FFT的结果所表示的谐波在频谱上的位置不再随基波频率变化而波动,可以根据算法带来的增益衰减进行固定补偿。

1 算法结构

算法包括过采样ADC、降采样低通滤波、DLL非均匀时钟产生、基波频率测量4个部分,如图1所示。过采样ADC装置、降采样低通滤波器、基波频率测量电路都工作在DLL生成的频率变化的时钟下。过采样ADC装置对模拟电压信号进行采样保持,并量化成1位的码流。降采样低通滤波器对接收到的1位码流进行低通滤波和降采样,同时完成ADC装置输出从1位码流到二进制补码并行数据的转换。FFT处理单元接收降采样之后的并行数据,进行谐波分析与测量。基波频率测量电路直接接收ADC装置输出的码流数据,使用基于正交信号恢复环路的算法,准确测量出当前的基波信号频率,作为非均匀同步采样时钟分频后的频率期望值控制DLL生成所需的非均匀时钟。

2 算法的实现和误差分析

2.1 过采样ADC

过采样ADC采用二阶Sigma-Delta调制结构,并采用1位比较器,输出1位码流。

其z域中的噪声传递函数HN(z)为:

${\rm H}{}_{\text{Ν}}(z)=\frac{(z-1){}^2}{z{}^2-1.4z+0.6}(1)$

设f0是输入信号的基波频率,ADC装置的平均采样频率fs=32 768f0。在0～3.2 kHz带宽内的量化噪声幅度可以达到-120 dB。

2.2 降采样低通滤波

由于64次谐波的带宽约为3.2 kHz,所以降采样低通滤波器中总的降采样率R=256可以使降采样后信号带宽为64f0,数据采样率为128f0,且带宽内的量化噪声幅度满足2.1节所述。

本文采用降采样率为64的3阶级联积分梳状(cascaded integrator comb,CIC)滤波器,级联2级半带滤波器(half band filter,HBF)的结构实现降采样低通滤波。在小于48f0的带宽范围内,滤波器的幅频增益最大衰减小于-0.45 dB,在对增益准确度要求不高的情况下,可以不进行补偿。

滤波器的传递函数如式(2)式(4)所示,其中,↓表示降采样操作。CIC滤波器的实现只包含加法和减法运算,HBF的系数经过分解后都可以用移位和加法的组合实现。

HDecim(z)=HCIC(z)(↓64)HHBF(z)(↓2)HHBF(z)(↓2) (2)

式中:HCIC(z)和HHBF(z)分别为CIC滤波器和HBF的降采样传递函数。

2.3 DLL非均匀时钟产生

设参考时钟频率fref=1.638 4 MHz,NDLL为延时单元个数,则DLL输出时钟的频率为:

$f{}_{\text{D}\text{L}\text{L}\text{o}\text{u}\text{t}}(l)=\frac{{\rm N}{}_{\text{D}\text{L}\text{L}}}{{\rm N}{}_{\text{D}\text{L}\text{L}}\pm l}f{}_{\text{r}\text{e}\text{f}}(5)$

式中:l为本次输出信号的延时单元位置与上次输出信号的延时单元位置之差,l∈{0,1,,NDLL-1},根据l取值方式的不同,DLL可以输出均匀时钟或者非均匀时钟。

fDLLout/32 768对应基波频率f0,根据式(5),其频率分辨率fres=62 Hz/NDLL。如果DLL输出的是均匀同步时钟,为了减少采样的频谱泄露,在要求fres<0.002 Hz的情况下,NDLL>31 000。如此大规模的延时单元阵列在电路中很难实现。

在本文中,DLL输出非均匀时钟,并送给ADC装置作为其采样时钟。此时fDLLout(l)为离散值,对位于fDLLout(l)和fDLLout(l+1)之间的平均采样频率值fs,DLL无法直接输出,只能得到在M个连续周期内统计意义上平均采样频率为fs的非均匀时钟。即

$f{}_{\text{s}}=\frac{{\rm M}}{{\displaystyle \sum _{i=0}^{{\rm M}-1}\frac{1}{f{}_i(d)}}}(6)$

式中:d∈{0,±1,,±(NDLL-1)}。

在式(6)条件下,DLL输出时钟的瞬时频率分辨率下降。为了降低非均匀采样噪声,根据3.1节中的要求,M<256。取M=250,可以将对fres的要求下降为fres<0.5 Hz,且保持M个周期内平均频率分辨率仍为0.002 Hz。此时,延时单元的个数NDLL>124,取NDLL=125。与输出均匀时钟相比,所需延时单元的个数大幅减少,易于电路实现。

在式(6)的条件下,根据式(5),d的期望E(d)必须满足式(7)所示的条件以达到同步采样的目的。

$E(d)=\left(\frac{f{}_{\text{r}\text{e}\text{f}}}{f{}_{\text{s}}}-1\right){\rm N}{}_{\text{D}\text{L}\text{L}}(7)$

d的不同取值方式有不同的分布类型,进而带来不同的采样噪声。从第3和第4节可知,对于每一个E(d)值,可以通过控制符合相邻二值分布的2个d出现的概率得到,在电路实现中,采用符合累加溢出原理的Fraction-N算法产生。

2.4 基波频率测量

针对采样时钟的非均匀特性,采用正交信号恢复环路的方法测量基波频率。环路由降采样滤波器、正交基波信号恢复、频率计算、积分型误差控制器构成,并使用2.1节和2.3节中的过采样ADC以及DLL,如图2所示。

降采样滤波器的降采样率为4 096,降采样后环路中平均采样率为fs/4 096=8f0′,其中,f0′为基波频率的测量值。

LUTS和LUTC分别是深度为8的sin和cos查找表,生成与fDLLout/4 096同频的sin和cos本地振荡。本地振荡与输入信号正交相乘,使用零极点相消的梳状滤波器(即图2中的Comb)滤除与直流和谐波相乘产生的和频与差频分量,仅保留与基波相乘得到的2个差频分量。再与本地振荡相乘,利用三角函数和差分解的逆运算[8],可以恢复出正交的基波信号s0和s*0。设k为降采样后的采样序列号,则

$\{\begin{array}{l}s{}_0=\frac{A{}_1}{2}\sin \left(2\text{π}\frac{f{}_0}{8f{}_{0^{\prime }}}k+\phi {}_0\right)\\ s{}_0^{*}=\frac{A{}_1}{2}\cos \left(2\text{π}\frac{f{}_0}{8f{}_{0^{\prime }}}k+\phi {}_0\right)\end{array}(8)$

式中:A1为基波的幅值;φ0为基波的初相位。

频率计算电路采用Cordic迭代算法对式(8)计算arctan(s0/s*0);对结果再进行差分运算,得到:

$c{}_{\text{ω}}=2\text{π}\frac{f{}_0}{8f{}_{0^{\prime }}}(9)$

积分型误差控制器的原型为仅保留积分环节的比例积分微分(PID)控制器;为了降低电路实现代价,使用累加器代替积分器,并引入反馈系数4/π以便从cω中得到f0′。根据PID参数的Z-N(Ziegler-Nichols)试验确定方法,最终确定积分参数ci1=0.125。积分型误差控制器的输出就是频率测量结果f0′,作为频率期望值输出给DLL。

在直流到64f0带宽范围内,以0.25f0为间距存在谐波的条件下,基波频率稳态测量误差小于0.000 001 Hz;在信噪比(SNR)为40 dB的噪声环境下,频率测量误差为0.000 8 Hz;被测量基波频率的动态范围为(50±4)Hz。

3 非均匀采样噪声分布

3.1 过采样ADC阶段

在非均匀采样序列为M个周期的情况下,文献[9]把非均匀采样序列分解为M个独立的均匀采样序列,插值叠加后得出了非均匀采样序列的频谱与原模拟信号频谱的关系。

如2.3节所述,fDLLout所确定的过采样间隔以M个采样点为周期进行变化,周期为MT,T=1/fs,即T为M个采样点内的平均采样时间间隔,fs为M个采样点内的平均采样频率。按照文献[9],非均匀过采样序列x(tn)的频谱Xd与原模拟信号频谱Xc的关系为:

$X{}_{\text{d}}(w)=\frac{1}{{\rm T}}{\displaystyle \sum _{k=-\infty }^{+\infty }A}(k,w)X{}_{\text{c}}\left(w-\frac{2\text{π}}{{\rm M}{\rm T}}k\right)(10)$

当把x(tn)作为均匀采样序列x(kT)进行频域变换$X{}_{\text{d}1}={\displaystyle \sum _{n=-\infty }^{+\infty }x}(t{}_n)\text{e}{}^{-\text{j}wn{\rm T}}$的时候,有

式中:tm为每一次的采样时刻。

单独考虑模拟输入信号中的每一次谐波分量。以基波ejw0t为例,非均匀采样序列x(tn)的频谱以fs为周期进行重复,在0 ～fs之间会以w0/(2π)为基准,以fs/M为间隔分布M个由于非均匀采样而产生的额外频率分量(即非均匀采样噪声),这些噪声的增益被A(k,w)调制。

此时,在满足M<R(R=256,为降采样低通滤波器总的降采样率)的条件下,可以保证由于使用DLL产生的非均匀时钟进行采样而产生的采样噪声位于降采样后的信号带宽64f0之外。

式(6)中d的均方差σ(d)越小,|rm|越小。考虑到实现代价,令d服从相邻二值分布,此时rm<0.011,fs>1.4 MHz((1±10%)fref),可得w0rm/fs<2.310-6。基波信号幅值被|A(0)|调制,有

$|A(0)|=\left|\frac{1}{{\rm M}}{\displaystyle \sum _{m=0}^{{\rm M}-1}\text{e}}{}^{\frac{-\text{j}w{}_{0}r{}_m}{f{}_{\text{s}}}}\right|\approx 1(13)$

根据Parseval定理,此时式(11)满足${\displaystyle \sum _{k=0}^{{\rm M}-1}|}A(k)|{}^2\approx 1$,所以可以近似认为:

$|A(k)|\approx |A(k+1)|\approx \sqrt{\frac{1-|A(0)|{}^2}{{\rm M}-1}}k\ne 0(14)$

可见,其都位于降采样后的信号带宽(64f0)之外。

模拟输入信号中其他谐波分量${\displaystyle \sum _{h=0}^{64}\text{e}}{}^{\text{j}w{}_{h}t}$(h为谐波次数)产生非均匀采样噪声的情况与基波的情况类似,总的频谱为各次谐波产生的采样噪声的累加。对于小于基波频率的间谐波分量,即wh<w0时,其对应的非均匀采样噪声落在63f0之外;对于大于基波频率的谐波分量,即wh>w0时,其对应的非均匀采样噪声落在64f0之外。所有的谐波分量的调制信号幅值$|A{}_h(0)|=\left|\frac{1}{{\rm M}}{\displaystyle \sum _{m=0}^{{\rm M}-1}\text{e}}{}^{\frac{-\text{j}hw{}_{0}r{}_m}{f{}_{\text{s}}}}\right|\approx 1$,h64。

3.2 降采样阶段

降采样低通滤波器对ADC装置输出的过采样数据进行降采样,这个过程依然是非均匀采样过程。式(6)中d满足独立同分布,降采样后信号的频谱为:

式中:Mds为降采样后采样间隔的变化周期,即降采样后Mds个采样点的平均采样频率为fs/R。

以基波ejw0t为例,非均匀采样序列降采样后的频谱以fs/R为周期进行重复。

在R mod M=0的情况下,Mds=0,即降采样后的采样序列为均匀采样序列,在0 ～fs/R之间不会产生非均匀采样噪声。

在R mod M≠0的情况下,在0～fs/R间会以w0/(2π)为基准,以fs/(RMds)为间隔分布Mds个由于非均匀降采样过程而产生的额外频率分量(即非均匀采样噪声),这些噪声的增益被Ads(k,w)调制。

由于rm,ds=rm/R<rm,所以将式(13)和式(14)中的A(k)替换为Ads(k),M替换为Mds,式(13)和式(14)在此阶段依然成立。

在式(6)中的d服从相邻二值分布的情况下,有

Mds=M (17)

总的频谱为各次谐波产生的额外频率分量的累加。

3.3 对滤波器的影响

从归一化频率和同步采样的角度来看,滤波器的特性与输入基波频率的相对关系不发生变化,即2.2节所述的增益衰减只需要进行固定补偿。

对滤波器特性的影响来自采样时钟的非均匀特性。非均匀特性对通带内衰减的影响等效于改变了滤波器的阶数,导致滤波器的增益变化。在式(6)中的d服从相邻二值分布的情况下,低通滤波器的截止频率fc变为$f{}_{\text{c}}\left(1\pm 1/{\rm N}{}_{\text{D}\text{L}\text{L}}\right)$,当NDLL=125的时候,第64次谐波会受到非均匀采样的影响,而第63次谐波不受影响。

4 仿真分析

为了保证由于使用fDLLout进行采样而产生的非均匀采样噪声位于降采样后的信号带宽64f0之外,需要满足M<R,取M=250。若要求fs/32 768的频率分辨率小于0.002 Hz,需要fDLLout/32 768的瞬时频率分辨率小于0.5 Hz,即需要NDLL>124,取NDLL=125。DLL的参考时钟频率fref=1.638 4 MHz。

式(15)要求式(6)中的d为独立同分布,式(13)和式(14)成立要求d的均方差σ(d)尽可能小,考虑到硬件实现代价,d满足d和d+1的二值分布可以满足上述要求,此时Mds=M。过采样ADC装置以及降采样低通滤波器2级非均匀同步采样对各次谐波的增益为|A(0)Ads(0)|=|A(0)|2≈1。

在上述条件下,使用DLL产生的非均匀时钟进行采样,可以达到近似理想的同步采样效果。2个非均匀采样阶段带来的非均匀采样噪声远远小于ADC的量化噪声,从而可以忽略,如图3所示。

图3中,基波信号的幅度为(20lg 0.5)dB,频率为51.112 Hz。获得图3的方法如下:①ADC装置使用过采样模型,并使用均匀采样时钟,获得理想量化噪声曲线;②ADC装置使用理想的采样保持模型,在非均匀采样时钟的条件下,得到非均匀采样噪声曲线;③ADC装置使用过采样模型,且使用非均匀采样时钟,得到量化噪声叠加采样噪声后的曲线。由于非均匀采样噪声远远小于理想量化噪声,③获得的曲线与①获得的曲线基本重合。

基波频率测量的范围决定了整个谐波测量算法中所允许的基波波动范围。在基波频率波动±8%(即46～54 Hz)的范围内,若ADC装置的测量满量程定为1,则得到的幅度测量误差如表1所示。不失一般性,在奇数次谐波误差之外,表中列出了直流、分数次间谐波以及偶数次谐波的误差。

表1中的数据根据降采样低通滤波器的增益衰减曲线对每一次谐波进行了固定的补偿。设基波幅度为220 V时取值为0.8(ADC装置的测量满量程定为1,下同),63次谐波幅度为0.22 V时取值为0.000 8,则基波的最大相对误差为3.7510-9,63次谐波的最大相对误差为0.157%。文献[10]给出了另外一种同步采样方法,其基波频率波动范围为±5%,最高可测量31次谐波。按照它的输入信号条件,本文算法对基波和31次谐波的测量误差与文献[10]相比分别小了4个数量级和2个数量级。

与文献[2,3,4,5,6,7]相比,本文提出的算法可以对采样数据直接进行FFT运算,无需消除非均匀采样噪声以及对数据进行重新处理或者消除频谱泄露等额外操作。在满足频率分辨率为f0/4且谐波相对误差小于0.1%的情况下,本文算法所需要的乘法运算量小于文献[1]中最简单的加窗插值算法所需乘法运算量的8%。与传统的基于锁相环(PLL)[8]的模拟同步采样方法相比,本文算法具有更稳定可靠、成本更低的优势。

5 结语

本文从减少运算量,易于超大规模集成电路硬件实现的角度,提出了一种基于非均匀同步过采样技术的谐波测量算法,简化了时钟产生电路的设计,减小了非均匀采样噪声,降低了对运算能力的要求,提高了测量精度。

参考文献

[1]邱海锋,周浩.非同步采样下电网谐波分析方法的探讨[J].继电器,2008,36(1):57-62.QIU Haifeng,ZHOU Hao.Study on the approaches ofelectrical harmonic analysis for asynchronous sampling[J].Relay,2008,36(1):57-62.

[2]李根,庞浩,徐建飞,等.一种基于改进锁相环的谐波分析逻辑电路[J].电力系统自动化,2007,31(21):82-85.LI Gen,PANG Hao,XU Jianfei,et al.Logic circuit forharmonic analysis based on enhanced PLL[J].Automation ofElectric Power Systems,2007,31(21):82-85.

[3]蔡菲娜,左伍衡.改进的双速率同步采样法及其傅里叶变换[J].浙江大学学报:工学版,2005,39(3):414-417.CAI Feina,ZUO Wuheng.Improved double speed synchronoussampling method and Fourier transform[J].Journal of ZhejiangUniversity:Engineering Science,2005,39(3):414-417.

[4]BLAND D M,LAAHO T I,TARCZYNSKI A.Analysis ofalgorithms for nonuniform time discreet Fourier transform[C]//Proceedings of IEEE International Symposium on Circuits andSystems,May 12-15,1996,Atlanta,GA,USA:453-456.

[5]JENQ Y C.Perfect reconstruction of digital spectrum fromnonuniformly sampled signals[J].IEEE Trans onInstrumentation and Measurement,1997,46(3):649-652.

[6]STOICA P,LI J,HE H.Spectral analysis of nonuniformlysampled data:a new approach versus the periodogram[J].IEEETrans on Signal Processing,2009,57(3):843-858.

[7]TAO R,LI B Z,WANG Y.Spectral analysis andreconstruction for periodic nonuniformly sampled signals infractional Fourier domain[J].IEEE Trans on Signal Processing,2007,55(7):3541-3547.

[8]KARIMI H,KARIMI G M,IRAVANI M R.Estimation offrequency and its rate of change for applications in powersystems[J].IEEE Trans on Power Delivery,2004,19(2):472-480.

[9]JENQ Y C.Digital spectra of non-uniformly sampled signals:fundamentals and high speed waveform digitizers[J].IEEETrans on Instrumentation and Measurement,1988,37(2):245-251.

非均匀采样 第7篇

在工业生产中存在着这样一类系统,它们的输入采样和(或)输出刷新呈现出不等时间间隔的非均匀的特点,这类系统被称为非均匀采样系统[1,2]。这类系统有的是由于硬件设备的限制、经济条件的制约或环境因素的影响等原因造成的被动非均匀采样系统,如产品组分和浓度的检测[3,4]、熔融指数的人工采样[5,6]和Kappa值的测量[7,8]等。还有的是为了在同样采样点数的情况下获得更多的有用信号而根据信号幅值变化情况主动调整采样间隔的主动非均匀采样系统[9]。近年来,非均匀采样系统的研究受到了国内外专家和学者的普遍关注,特别是作为状态估计和控制理论基础的系统辨识领域[10,11,12,13,14,15,16,17,18]。

在白噪声的前提下,传统算法,如最小二乘算法等,都能给出非均匀采样系统的参数无偏估计。而在工业实践中的扰动都是有色噪声,传统的最小二乘等算法不能满足无偏估计的要求[19]。为了获得参数的无偏估计,研究人员提出了一系列的算法。例如,基于滤波的算法[20]、最大似然算法[21,22]、基于辅助模型的算法[23]等。基于滤波的算法首先采样一个动态的滤波器来白化有色噪声,然后使用传统算法,如最小二乘算法等,获得无偏估计。但是当工业过程输出信号噪信比较大或者参数较多时这类算法会出现估计偏差[24]。并且最大似然类算法给出的估计值都是渐近无偏的,也就是说最大似然算法并不能保证给出的所有估计都是无偏的。此外,最大似然算法的计算量比较大,也在一定程度上限制了其应用[24]。基于辅助模型的算法,如基于辅助模型的递推最小二乘算法(AM-RLS)和基于辅助模型的随机梯度(AM-SG)等,可以同时对系统模型和噪声模型的无偏估计,是一类比较有应用前景的算法。但在辨识过程中,AM-RLS和AM-SG都没有把系统噪声的方差考虑进辨识算法以提高算法的辨识精度,所以本文的基本思路是考虑结合噪声的方差,基于辅助模型的思想和递推最小二乘的形式提出一种递推贝叶斯RB(recursive Bayesian)辨识算法。并用输入非均匀采样广义输出误差模型NUSID-GOE(non-uniformly sampled input data for generalized output error model)验证了算法的有效性。

1 问题描述

考虑如图1所示的输入非均匀更新和输出定期采样的NUSIDGOE模型。图中,包括零阶保持器Hτ、线性时不变过程Pc和采样器ST三个部分。其中,u(t)和x(t)分别是输入信号和无噪输出信号;y(t)为有噪输出;w(t)是有色噪声;y(k T+T)是采样时刻t=k T+T时刻的输出采样值,其中T为输出采样周期,也称框架周期。

假定输入非均匀采样时间间隔为{τ1,τ2,…,τr},它们满足如下的条件:T=τ1+τ2+…+τr。定义t0=0,ti=ti-1+τi,其中i=1,2,…,r,那么输入u(t)可以写为:

其中k=0,1,2,…。

如图1所示,若Pc和w(t)采用文献[20]的形式,则NUSIDGOE模型的连续输出y(t)可以表示为:

其中:

定义参数向量如下:

其中:

ns=na+rnb+1,nn=nc+nd,n=ns+nn;下标s和n分别表示过程模型和噪声模型。

类似地,分别定义信息向量φ如下:

对于t=k T时刻,等式可以改写成:

其中:

假定系统的模型结构是已知的,即(na,nb,nc,nd)和输入采样次数r均是已知的,那么算法的目的就是如何利用观测数据D(k)估计参数向量θ。

2 递推贝叶斯参数辨识算法

贝叶斯方法的主要观点就是把待估计参数作为一个随机变量,通过观测与该变量相关的其他随机变量来推测该变量的值[25]。在提出的算法中,以最大化后验概率密度作为估计准则,即:

根据贝叶斯公式,参数的后验概率密度函数表示为:

其中p(y(k)|θ,D(k-1))为在给定参数θ和数据条件下输出y(k)的先验概率密度。

但是式并不能直接用于计算p(θ|Dk),因为数据的先验概率密度和p(y(k)|θ,D(k-1))和参数的概率密度p(θ|D(k-1))均是未知的。对于后者,一个可行的假设是p(θ|D(k-1))为均值为方差为P(k-1)的正态分布,其概率密度函数为:

对于p(y(k)|θ,D(k-1))未知的情况,假定噪声模型输入v(t)为均值为零方差为σv2的白噪声,即v(k)~N(0,σv2),那么根据式(13),可以得到y(k)~N(φT(k)θ,σv2),因此数据的先验概率密度函数为:

于是参数的后验密度函数p(θ|Dk)的表达式为:

又因为式等价于,求导后整理得到:

其中:

这样就得到了递推贝叶斯(RB)算法的简洁形式:式(22)和式(23)。为了避免计算P(k)的逆矩阵以减小算法的计算量,对式(23)运用矩阵求逆反演公式,得到:

为方便,定义增益向量,RB算法的常用形式如下所示:

标准的RB算法并不能用于参数的估计,原因是:信息向量φ(k)中的x(k T-i T)项和噪声项w(k T-i T)及v(k T-i T)、噪声的方差项σv2均是未知的。可以采用辅助模型的思想[16]:用辅助模型的输出代替x(k T-i T),分别用w(k T-i T)和v(k T-i T)的估计值代替信息向量中的相应项,定义:

其中:

考虑到图1的模型,根据式(28)-式(30),可以表示为:

噪声方差可以用下式进行估计:

于是,RB算法可以整理如下:

RB算法的计算步骤如下:

Step1

收集输入输出数据u(k T+ti)和y(k),i=0,1,2,…,r,k=1,2,…。

Step2

令u(k T+ti)=0,i=0,1,2,…,r-1,y(k)=0,

Step3

为参数和矩阵P设置初值如下:

其中1n为元素全为1的n维列向量。

Step4

分别根据式(29)、式(30)和式(28)构造。

Step5

利用式(34)估计,利用式(36)计算增益矩阵L(k)。

Step6

分别使用式(35)和式(37)计算参数向量和协方差矩阵P(k)。

Step7

从参数向量中分离出过程参数,并分别利用式(31)-式(33)估计。

Step8

令k=k+1,转step4。

3 收敛性分析

定理1对于系统式(13)和RB算法式(22)、式(23),假设{v(t),Ft}是定义在概率空间{Ω,F,P}上的鞅差向量序列,其中{Ft}是截至时刻t之前所有观测数据生成的σ代数序列[26];如果噪声序列{v(t)}满足如下假设:

并且存在正常数c1,c2和t0使得下面的持续激励条件成立[27]:

那么算法给出的参数估计值收敛于真值,即:

其中θ0为参数真值。

证明在式两边同时减去θ0并用φT(k)θ0+v(k)代替y(k),得到:

其中:

众所周知,,w.p.1,并考虑到式(38),可得为式(41)的不变集[19]。那么问题式(40)转化为差分方差的稳定性问题。

在式(23)两边同乘以P(k),则Q(k)可以写成:

令λ为矩阵Q(k)的一个特征值,p为其对应的特征向量,那么下式成立:

考虑到式(42)、式(43)可以写成:

在式(44)的两边同乘以P-1(k),考虑到式(23),并整理得:

在式(45)两边同乘以pT并整理得:

考虑到P(k)的对称正定性和式(39),P-1(k-1)和φ(k)φT(k)均是正定的,所以均大于0。从而式(46)中的项(1-λ)和λ是正负同号的,也即,这也就是说Q(k)的特征值满足0<λ<1。根据差分方差稳定性理论,差分方差式(43)是稳定的,即,,因此定理1成立。

4 仿真实例

考虑如下的NUSIDGOE模型:

其中:

仿真时,输入序列{u(t)}采用零均值单位方差的不相关随机信号,噪声序列{v(t)}采用零均值方差为σv2=0.82白噪声。分别使用RB算法和AM-RGELS算法[20],对系统进行了辨识,辨识结果和估计误差见表1、表2所示,其中参数估计误差定义为。图2为两种算法参数估计值和估计误差的对比情况。

根据表1、表2的数据和图2的仿真结果可以看出:

(1)数据长度增加时,两种算法的估计误差都在变小;

(2)对同样的数据长度,RB算法比AM-RGELS算法的估计精度高;

(3)在相同噪声水平下,RB算法比AM-RGELS算法的收敛速度要快。

5 结语

本文针对输入非均匀采样广义输出误差模型提出了一种递推贝叶斯参数估计算法。它以通过最大化参数的条件后验概率密度函数作为估计准则,由于不仅考虑了参数的先验概率密度而且还考虑了数据的先验概率密度。因此该算法与单纯考虑数据的概率密度的极大似然类参数估计方法相比,精度更高。同时,该算法的辨识精度也会比传统最小二乘算法略高,因为它在辨识过程中考虑了噪声的方差。综上所述,该算法把给定数据集下最可能出现的参数值作为参数的估计值。收敛性分析表明,所提算法的参数估计值能收敛到真值,仿真结果表明了算法的有效性和较高的估计精度。

摘要：针对传统最小二乘算法在辨识过程中没有考虑噪声的协方差和参数的先验概率密度的问题,提出一种递推贝叶斯算法。该算法以最大化参数的后验概率密度函数为准则进行参数估计。实验结果证明所提算法可以获得更高精度的参数估计值。收敛性分析表明,该算法给出的参数估计值收敛于参数真值。该算法综合考虑了噪声方差、数据的先验概率分布和参数的先验概率分布,可以获得比最小二乘法更高的精度的估计值。