CN118089533A - 基于频偏修正的多干涉谐波叠加信号频率及相位参数重建 - Google Patents

基于频偏修正的多干涉谐波叠加信号频率及相位参数重建 Download PDF

Info

Publication number
CN118089533A
CN118089533A CN202410133247.4A CN202410133247A CN118089533A CN 118089533 A CN118089533 A CN 118089533A CN 202410133247 A CN202410133247 A CN 202410133247A CN 118089533 A CN118089533 A CN 118089533A
Authority
CN
China
Prior art keywords
frequency
sequence
phase
spectrum
value
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202410133247.4A
Other languages
English (en)
Inventor
常林
高杰华
闵俊一
吕璟辉
崔萧管
徐子康
尹健
费云夕
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Huzhou University
Original Assignee
Huzhou University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Huzhou University filed Critical Huzhou University
Priority to CN202410133247.4A priority Critical patent/CN118089533A/zh
Publication of CN118089533A publication Critical patent/CN118089533A/zh
Pending legal-status Critical Current

Links

Landscapes

  • Spectrometry And Color Measurement (AREA)

Abstract

本发明公开了一种基于频偏修正的多干涉谐波叠加信号频率及相位参数重建方法,可以避免传统傅里叶变换频谱分析中欠采样和非同步采样等问题导致的谐波频率和幅值求解不准确,本专利主要包括三个主要的测量步骤。本发明所采样的近似熵计算法和频率幅值修正法,能够实现非同步采样问题的量化评估以及高精度地对谐波频域特征进行分析,所设计的欠采样时频谱折叠的判定方法使所设计的技术的适用范围更加广泛。本专利所设计的方法具有测量精度高、可操作性强、对测量误差鲁棒性强等特点,可用于透明多表面光学元件的波长移相干涉测量中。

Description

基于频偏修正的多干涉谐波叠加信号频率及相位参数重建
技术领域
本发明涉及一种基于频偏修正的多干涉谐波叠加信号频率及相位参数重建方法,结合近似熵量化计算了移相干涉测量中的非同步采样问题,利用序列重排技术实现了具有更低频谱泄漏的谐波频率分析,并且据此实现了多表面干涉谐波的幅值、频率修正和欠采样识别,应用在精密激光移相干涉测量领域。
背景技术
干涉测量是一种高精度的测量技术,通过将被测波形与参考波形进行干涉,进而通过一定的方法提取被测波面与参考波面的初始相位,最终实现对被测面表面的精确测量。其过程如下:首先,将被测件放置在干涉仪前方,使用干涉仪将被测波面与已知波面进行干涉,形成干涉图样,并通过CCD相机等硬件进行干涉图的存储。然后,通过测量干涉图样的变化(通常需要进行移相),可以确定被测波的波长、相位和振幅等参数。
具体而言,当被测波面与已知波面(也可称为参考波面)在干涉仪中相遇时,它们会发生干涉,形成稳定的干涉条纹图样,人为引入主动的相位变化值(也即移相)获得一系列具有相同初始相位、同时也具有不同的主动相位的干涉条纹图(也可以称为干涉强度图,通常简称为干涉图)作为数据基础。由于干涉图样的光强变化值与被测波面的波长、相位和振幅等参数密切相关,因此基于采样得到的具有不同相位值的干涉光强图,通过相应的参数解调方法,可以实现上述参数的精确求解。干涉测量术属于非接触式光学测量法,具有高精度、高灵敏度等优点,广泛应用于微纳结构识别与声、光检测等领域。
在主动相位引入方式方面,传统的硬件移相法通过改变干涉仪前方的参考镜的微观距离从而改变被测件和参考镜之间的光程差(距离值与材料折射率的乘积),可以实现主动相位的改变。但是这种移相方式存在较多问题,具体如下:
(1)相移器稳定性问题:硬件移相法需要使用压电结构推动参考镜(称为移相器),而相移器的稳定性直接影响干涉测量的精度。但是,在实际测量时,相移器容易受到温度、机械振动和老化等因素的影响,从而导致测量误差;
(2)硬件移相的实现较为复杂:硬件移相法需要进行对压电结构进行主动电压输入和输出修正等操作,复杂度较高,需要使用专门的电路和设备,增加了测量成本和测量难度;
(3)测量范围有限:硬件移相法的测量范围受到相移器移相范围的限制,对于某些需要大范围或者大口径相位调制的干涉测量,硬件移相法可能无法满足要求;
(4)难以实现高精度相位测量:受到上述不良因素的影响,硬件移相法的相位测量精度也受到限制,难以实现高精度相位测量;
(5)多表面叠加信号无法求解的问题,当被测件是透明且前后表面具有较高平行度的元件时,采集得到的干涉图是各表面干涉的叠加结果,硬件移相法无法处理该种干涉条纹图。而这种多表面元件在光学系统中应用广泛,正是本专利所针对的测量目标。
而近些年获得较快发展的波长移相法可以较好地克服上述问题。波长移相干涉测量技术是通过改变光源的波长进行主动相位的改变(波长与主动相位呈线性关系),也即移相,并在改变波长的同时进行干涉图采集。
这种移相方法的优势在于:激光器作为移相器和光源,不需要引入具有复杂结构的硬件移相装置,因此干涉仪结构较为简单,可以克服硬件移相所具有的机械迟滞效应。并且由于避免了参考镜的推动,因此该种移相方法可以实现大口径被测件的测量。最重要的是,在对具有多个表面的透明平板元件进行检测时,因为平板元件前后表面平行度较高,各被测面与参考镜中的参考件均会产生干涉,采集得到的干涉图是各个被测面表面干涉信息叠加的结果,传统的硬件移相法得到的干涉图无法进行参数解调。但通过波长移相法得到的干涉图中,各个被测面对应干涉信号(称为谐波)因为具有不同的光程差,因此可以在时域中展现出不同的谐波频率,这为多个被测表面的同时重建提供了可能性。
在干涉参数的求解方面,通过传统的离散傅里叶变换对谐波信号的频率和相位进行求解是一种常用的经典方法。但是这种方法的缺陷在于:由于波长移相过程中(也即采样过程中),波长的改变量并非完全线性的,这就会造成非同步采样问题。直观而言,这是由于谐波信号的采样长度与信号周期不成整数倍关系所导致的,完全的同步采样在实际测量中受到各种误差的影响,也是不可能存在的。而非同步采样问题的存在,会导致传统离散傅里叶变换后的频率幅值谱的频谱泄漏,致使谐波频率、相位和幅值的求解精度受到损失,非常不利于高精度检测的实现。同时,在实际测量过程中,非同步采样问题是难以量化识别的。
此外,根据多频谐波采样定理,当谐波频率大于的采样频率一半时,则在单边频率谱中,该种高频信号的频率曲线会以二分之采样频率为轴线折叠到探测到的频谱中,这种现象可以称之为欠采样频率模糊现象或频谱折叠问题,此时探测到的频率是假频,基于这种频谱得到的频率和相位等参数信息,会致使测量的失效,无法反映被测件的真实面形特征。
在目前已存在的技术和装置中,专利发明人常林等人在发明专利《一种用于光学多表面透明板夹持的可旋转夹具-申请号:202010595035.X》、《一种用于光学透镜的一体式夹持及夹紧机构-申请号:202110177033.3》以及论文《Lin Chang,Yingjie Yu,et al.Ageneral auto-shift minimal-step phase-shifting algorithm for arbitrary cavitylength[J].OPTICS AND LASERS IN ENGINEERING,2022,149:106791.(SCI二区,第一作者)》、《Lin Chang,Yingjie Yu,et al.Multi-surface phase-shifting algorithm usingthe window function fitted by the nonlinear least squares method[J]JOURNAL OFMODERN OPTICS,2022,69(3):160-171.(SCI四区,第一作者)》中所设计的新型系统装置和测量方法实现了被测件的便捷夹持、灵活位置下的各表面光学检测等目标。
但是,上述专利和论文中的创新结构和方法无法解决上面所指出的问题,具体包括:
(1)存在非同步采样问题的量化识别;(2)具有更低频谱泄漏的频率分析的实现;(3)精确的干涉光强对应频率、幅值修正;(4)精确的各测量表面的同时重建;(5)欠采样频峰模糊问题发生时的高精度识别。
发明内容
为解决上述问题,本发明设计一种基于频偏修正的多干涉谐波叠加信号频率及相位参数重建的方法,可以以更低的频谱泄漏实现高精度的谐波频率求解和各表面同时面形重建,并可以通过近似熵对非同步采样程度进行量化。
对材质透明且前后表面具有高平行度的光学元件进行检测的大致过程可以归纳为:将被测件夹持在夹具上并且放置在干涉仪前方使待测元件与干涉仪中的参考镜平行以产生强度叠加干涉图、打开干涉仪并且利用外部控制器驱动光源改变波长实现移相同时采集干涉图、对干涉图进行处理获得各表面对应的包裹相位、对包裹相位进行解包裹和消倾斜等操作最终实现各表面对应波面高度的重建。
本发明所设计的测量方法的主要流程如下:
步骤1:通过波长调谐移相干涉仪采集2N-1帧干涉图,其中N应为不小于20的整数,选取干涉图中像素点(xp、yp)的光强记为Ip进行分析,设定近似熵阈值且对Ip进行近似熵计算,步骤1主要包括(1.1)向量序列分组、(1.2)向量距离量化计算、(1.3)量化向量相似度、(1.4)近似熵量化求解、(1.5)主频率搜索5个主要过程,具体计算步骤如下:
优选地,近似熵(Approximate Entropy,ApEn)主要用于评估信号序列的复杂性和不规则性的参数。时间序列信号复杂度与近似熵的结果成正比,也即,越规则的序列则近似熵越小。这种计算方法的一个显著的优势在于它对数据量的需求较小,且计算结果稳健、可靠,这意味着大多数实际测量的时间序列都满足其要求,从而可以获得稳定可靠的结果,是对干涉测量中的光强序列规则性进行量化的一种合适的方法。
(1.1)向量序列分组:将光强Ip带入x(na)中,对于由NA个光强数据组成的时间序列x(na)=x(1),x(2),…,x(NA),近似熵的计算方法如下:
解释上述过程如下,各光强数据也即:某像素点在不同帧序数(na)下的光强变化即为一时间序列,NA即为此时对应的总帧数。
按序号组成一组维数为m的向量序列G,具体表示为其中/> 向量序列G代表从第ia点开始的m个连续的x(na)中的元素值,例如当m=2时,/>
(1.2)向量距离量化计算:优选地,定义向量与/>的距离/>为:
上式中的max为求解最大值,||为计算绝对值。
(1.3)量化向量相似度:优选地,对于一组与/>计算/>小于等于误差半径ra的数量并将结果记为ci,进而计算εm,r
其中,m和ra可以根据需求进行设定,总体误差变动范围用GA表示,数值上GA=2ra,但设定ra时应满足ra处于0.1FSD(x)和0.25FSD(x)之间的要求,其中FSD(x)为序列的x(na)的标准差;m一般设定为2。
(1.4)近似熵量化求解:优选地,增加维数到m+1,重复步骤(1.1)~(1.3),得到εm+1,r,并通过下式计算近似熵FApEn(m,r):
FApEn(m,r)=εm,rm+1,r (式-3)
(1.5)主频率搜索:对光强信号进行IP补零操作,补零长度为5N,补零后序列的总长度为7N-1且记为IP0,对IP0去均值后进行傅里叶变换并取单边幅值谱进行3个主要谐波峰值对应横坐标搜索,且将各个谐波峰值的横坐标记为主频率值。
至此,步骤1计算过程结束,如果近似熵大于预设熵值,则转至步骤2。
步骤1所设计的干涉光强序列的近似熵计算方法,其有益效果是能够有效识别非同步采样问题,可以为后续的进一步频谱分析和处理提供可靠的量化基础。
步骤2:优选地,将光强Ip带到x(n)中以后进行下一步处理,步骤2主要包括(2.1)构造重排序列、(2.2)序列重排FFT函数求解、(2.3)加窗序列频谱函数求解、(2.4)构造相位差、(2.5)频率和幅值修正五个具体步骤,流程如下:
以x(n)(此时,只截取n≥0的部分进行分析,n=0,1,…,N-1为时间序列的序号)为单频复指数信号作为示例展示序列重排傅里叶变换的计算过程:
上式中,信号的数字频率ω0表示为β倍频率间隔2π/N的形式,其中β可以是小数、θ为初始相位、N为总序列点数,j为虚数单位,则经过归一化以后,x(n)的原始FFT频谱函数(FFT即为快速傅里叶变换,是离散傅里叶变换的一种快速计算形式)可表示为:
上式中,k=0,1,…,N-1。进而可将上式展开为:
(2.1)构造重排序列:序列重排FFT是考虑包含序列中某点在所有循环移位后的分段傅里叶变换而产生的一种分析方法。
优选地,所述的序列重排是对包含着序列中的某个点的向量进行重新排序,进而进行循环抽取和排列,以x(0)为例,存在且只存在N个包含该点的N维向量:
优选地,上述公式中,x0到xN-1右下角的角标表示序列重排后的向量序号,将每个向量进行循环排序且把序列值x(0)移到序列中的第一位,则可得到另外的N个N维向量:
优选地,对准x(0)相加并取其平均值,则可得到重排后的数据向量,以xap表示如下:
(2.2)序列重排FFT函数求解:根据傅里叶变换的移位性质,以角标ka表示向量序号,则x’ka(n)的离散傅里叶变换X’ka(k)和x(n)的原始FFT频谱函数Xka(k)之间关系如下:
优选地,对上式的X′k(k)进行求和平均即为序列重排FFT函数:
(2.3)加窗序列频谱函数求解:对于长度为N的序列,若用窗函数Wf对单频复指数信号进行加窗FFT(注意,这里仍然是以n≥0为例截取长度为N的序列值进行分析和展示),以函数Fg为窗函数的频域形式,则得到的离散传统FFT的谱表达式为:
上式中,Δω=2π/N为频率分辨率单位,群延迟量为τ=(N-1)/2,根据FFT的LTI特性,对幅值为A、初相为θ、时间序列延迟为n0、角频率为ω*的复指数信号 再进行加窗FFT谱分析,则其离散谱表示为:
若对序列进行N个值和后N个值都执行加窗操作,则前、后窗函数均为Wf的双窗序列重排FFT的频谱函数表示为:
则根据序列重排FFT的线性特性,对信号进行同样的双窗序列重排FFT谱分析,其离散谱表示为:
(2.4)构造相位差:优选地,以k*代表频率幅值谱中主频率谐波峰值的主谱线序号值,进而将传统相位谱值减去序列重排后的相位谱值/>则可得到如下的相位差表达式如下:
(2.5)频率和幅值修正:优选地,测出传统FFT和序列重排的两种谱分析在主谱线上的相位大小并作差,便可计算出频谱分析得到的频率与实际频率之间的偏移值,进而通过下式对信号真实频率的估计:
此处注意,若不只取主谱线的相位,也可以绘制相位差值的曲线。
分别取(式-13)和(式-15)的主谱线上的模值,有
将(式-15)进行进一步变换,即可得到频率幅值的修正值:
由(式-17)至(式-19)就可完成光强信号IP中3个谐波的频率、幅值估计,以频率值增序频排列,则经修正后的3个谐波的频率和幅值分别记为fp1、fp2、fp3和Ap1、Ap2、Ap3
在上面的计算中,窗序列Wf可以根据用户需求进行选取,常用的窗函数有汉宁窗、海明窗、矩形窗、三角窗、布莱克曼窗等,默认选择汉宁窗。
步骤2所设计的干涉光强序列的序列重排FFT以及频率、幅值修正方法,其有益效果是能够显著降低频谱泄漏,并且可以得到精确的谐波频率和幅值信息,可以为后续的相位解调、波面重建以及欠采样的判定打下坚实基础。
至此,步骤2计算过程结束,转至步骤3。
步骤3:频谱折叠判定及相位求解:当(Ap2-Ap3)/Ap1大于折叠阈值EJS时,则判定不存在频谱折叠,若(Ap3-Ap2)/Ap1大于折叠阈值EJS则存在频谱折叠的问题,若满足|(Ap2-Ap3)/Ap1|≤EJS则无法准确判断是否发生频谱折叠,将像素点移动到(xpp、ypp)后再次执行步骤1和步骤2;
折叠阈值可根据用户需求进行设定,一般应不小于0.2且不大于0.5。
优选地,当存在频谱折叠时,FQS-1=fp1为厚度变化谐波的相位计算频率值、FQS-2=fp2为前表面谐波的相位计算频率值、FQS-3=fp3为后表面谐波的相位计算频率值;
优选地,当不存在频谱折叠时,FQS-1=fp1仍然为厚度变化谐波的相位计算频率值、FQS-2=fp3为前表面谐波的相位计算频率值、FQS-3=fp2为后表面谐波的相位计算频率值;
为各谐波的初始相位、Bp和Ap分别是各谐波对应的相位求解函数,进而干涉图中的初始相位的求解公式可以写为:
上式中U0表示干涉图中的光强,U0中包含着5N帧元素全为0的干涉图。各谐波对应的Bp和Ap分别通过下面式子确定:
上式中Re和Im分别表示求实部和求虚部,FQS-p、BP和AP中的右下角角标p=1、2、3时分别代表厚度变化信号、前表面信号和后表面信号。应注意,此时(式-20)和(式-21)中的n=0,1,2,...7N-1。
求解出各谐波的包裹相位以后,对包裹相位进行解包裹、消倾斜以后便可根据波面高度和相位之间的关系进行各谐波面形的重建。
步骤3所设计的相位解调和欠采样识别方法,其有益效果是能够高效实现干涉图中的初始相位提取,并且可以根据前面计算得到的精确的谐波频率和幅值信息,对欠采样问题进行有效识别和处理,对于精确的多表面波面重建是非常便利的。
至此,步骤3计算过程结束,测量完成。
本发明与现有技术相比较,具有如下显而易见的突出实质性特点和显著优点:
(1)结合所采用的近似熵计算法和频率幅值修正法,能够实现非同步采样问题的量化评估;
(2)所设计的欠采样时频谱折叠的判定方法使所设计的技术的适用范围更加广泛,有效克服欠采样频谱折叠问题发生时的相位解调和频率对应;
(3)所采用的频率和幅值修正方法能够显著提高谐波参数重建精度;
(4)所采用的序列重排数据组合方法能够有效降低频谱泄漏;
(5)通过本专利所采用的近似熵求解、频率幅值修正和相位解调的组合方法,可以通过一组干涉图同时对被测件的多个波面进行高精度重建。
附图说明
下面结合附图和实例过程对本发明做进一步说明。
图1为干涉仪模型示意图及样机实物图;
图2为干涉仪内部光路图;
图3为多表面干涉时的光强叠加示意图;
图4为实际采集的3帧光强叠加干涉图;
图5为近似熵计算过程示意图;
图6为传统离散傅里叶变换得到的离散频率幅值谱;
图7为补零与频谱折叠示意图;
图8为同步与非同步采样情况下近似熵判定结果实例图;
图9为频率及幅值修正示意图;
图10为频谱泄漏对比及误差求解实例;
图11为光强补零示意图;
图12为本发明所设计的方法的求解实例;
图13为本发明所设计的方法的求解误差分布;
图14为本发明实现流程图。
具体实施方式
本发明的优选实施例结合附图详述如下:
实施例一:
在本实施例中,对本发明所设计的测量方法进行详细解释和阐述。波长调谐移相干涉测量所用到的干涉仪模型示意图及样机实物图如附图中的图1所示,其内部光路图如图2所示,由光源发出的一定波长的光束在被测件和参考镜表面发生反射并产生干涉图,由干涉仪内部的相机(CCD)进行干涉图的采集,如果被测件前后表面均与参考镜的参考面平行度较高,则产生如附图中图3所示的叠加干涉图,连续3帧实采干涉图如附图中的图4所示。
从干涉图中某个像素点的角度出发,观察其在不同帧数下的光强变化可发现,这种光强变化遵循余弦信号特征;而多表面干涉光强则为3个具有不同频率的谐波的叠加,因此采用时间信号分析中的频域分析法是非常合适的。而本专利所采用的复指数信号的分析方法,能够和余弦函数分析法进行有效结合,因二者之间可通过欧拉公式进行转换,并且复指数信号和余弦信号均以分析得到信号的频率、相位和幅值为目的,频谱表现是一致的。
本发明所设计的测量方法的主要流程如下:
步骤1:通过波长调谐移相干涉仪采集2N-1帧干涉图,其中N应为不小于20的整数,选取干涉图中像素点(xp、yp)的光强记为Ip进行分析,设定近似熵阈值且对Ip进行近似熵计算,步骤1主要包括(1.1)向量序列分组、(1.2)向量距离量化计算、(1.3)量化向量相似度、(1.4)近似熵量化求解、(1.5)主频率搜索5个主要过程,具体计算步骤如下:
近似熵主要用于评估信号序列的复杂性和不规则性的参数。时间序列信号复杂度与近似熵的结果成正比,也即,越规则的序列则近似熵越小。这种计算方法的一个显著的优势在于它对数据量的需求较小,且计算结果稳健、可靠,这意味着大多数实际测量的时间序列都满足其要求,从而可以获得稳定可靠的结果,是对干涉测量中的光强序列规则性进行量化的一种合适的方法。
(1.1)向量序列分组:将光强Ip带入x(na)中,对于由NA个光强数据组成的时间序列x(na)=x(1),x(2),…,x(NA),近似熵的计算方法如下:
解释上述过程如下,各光强数据也即:某像素点在不同帧序数(na)下的光强变化即为一时间序列,NA即为此时对应的总帧数。
按序号组成一组维数为m的向量序列G,具体表示为其中/> 向量序列G代表从第ia点开始的m个连续的x(na)中的元素值,例如当m=2时,/>
(1.2)向量距离量化计算:定义向量与/>的距离/>为:
上式中的max为求解最大值,||为计算绝对值。
(1.3)量化向量相似度:对于一组与/>计算/>小于等于误差半径ra的数量并将结果记为ci,进而计算εm,r
其中,m和ra可以根据需求进行设定,总体误差变动范围用GA表示,数值上GA=2ra,但设定ra时应满足ra处于0.1FSD(x)和0.25FSD(x)之间的要求,其中FSD(x)为序列的x(na)的标准差;m一般设定为2。
解释上述过程如下,根据定义,εm,r的本质是的对数表示下的平均值,ra的设定取0.2是一种有效的量化计算方法。m的选择可以决定计算近似熵的时序列的截取长度。通常选择m=2更为合适,因为这样可以获得更多的信息。避免选择m>2,因为这需要序列长度可达到数千个点,且有可能会无法正确体现序列本身的性质。
(1.4)近似熵量化求解:增加维数到m+1,重复步骤(1.1)~(1.3),得到εm+1,r,并通过下式计算近似熵FApEn(m,r):
EApEn(m,r)=εm,rm+1,r (式-3)
(1.5)主频率搜索:对光强信号进行IP补零操作,补零长度为5N,补零后序列的总长度为7N-1且记为IP0,对IP0去均值后进行傅里叶变换并取单边幅值谱进行3个主要谐波峰值对应横坐标搜索,且将各个谐波峰值的横坐标记为主频率值。
上述近似熵计算示意图如附图中图5所示。至此,步骤1计算过程结束,如果近似熵大于预设熵值,则转至步骤2。
解释上述过程如下,可以从上述计算过程看出,近似熵的本质是序列在m维和m+1维的近似比例的对数表示下的平均值的差值,若差值越大,序列越不规则、越复杂,反之亦然。在进行干涉图采样过程中,由于严格的同步采样难以实现,因此多数情况下实采干涉图的光强序列均为非同步采样序列,而近似熵的计算可以有效地量化非同步采样的程度,对后续的参数重建是有利。近似熵判定的预设熵值可根据用户需求选择,一般而言,默认选择0.3,因同步采样情况下的近似熵值通常小于0.1。
步骤2:将光强Ip带到x(n)中以后进行下一步处理,步骤2主要包括(2.1)构造重排序列、(2.2)序列重排FFT函数求解、(2.3)加窗序列频谱函数求解、(2.4)构造相位差、(2.5)频率和幅值修正五个具体步骤,流程如下:
以x(n)为单频复指数信号作为示例展示序列重排傅里叶变换的计算过程,若x(n)为一长度为N的序列(实际上采样光强序列的长度为2N-1,但是这里截取n≥0的部分进行分析展示,这样更符合傅里叶变换的特点),且其时域形式表示为:
上式中,信号的数字频率ω0表示为β倍频率间隔2π/N的形式,其中β可以是小数、θ为初始相位、N为总序列点数,j为虚数单位,则经过归一化以后,x(n)的原始FFT频谱函数(FFT即为快速傅里叶变换,是离散傅里叶变换的一种快速计算形式)可表示为:
上式中,k=0,1,…,N-1。进而可将上式展开,可以得到传统傅里叶变换以后长度为N的序列的频谱函数:
附图中图6所示的是经典的快速傅里叶变换得到的单边频谱图,可以看到,由于数据点较少,因此可供进行频谱分析的有效频点也较少,极易受到频谱泄漏等不良因素的影响。并且也可以看到,图6所示的频谱图无法有效量化评估非同步采样问题和欠采样产生的频谱折叠问题,图7所展示的折叠后的频谱和探测到的频谱中,可以清晰地展现这种折叠问题。而使用本方法中的步骤1中的近似熵计算,可以有效地获得非同步采样量化结果,附图中图8所示。在图8中,左图做展示的光强序列中3个谐波的频率分别是f1=5.2、f2=10.2、f3=25.3,而采样频率则为50,显然,在这种情况下,各谐波频率与采样频率以及采样点数之间均不是整数倍关系,因此发生了非同步采样问题,此时近似熵结果为0.8723,显著大于右图中所展示的谐波的频率分别是f1=5.0、f2=10.0、f3=25.0所对应的近似熵值0.0691,这也充分验证了本专利所提出的近似熵方法的有效性。
解释上述过程如下,对序列长度和序列中的序号表示的不同在于分析角度不同,是为下面的频率分析的便捷性考虑的,如na在计算近似熵中表示序号、NA表示总的序列点数;而在序列重排中,以n表示序号、以N表示总的序列点数;转换到频谱函数分析时,为了加以区分,以k表示经过傅里叶变换后的频谱函数中的序号。以x(n)和x(na)作为分析对象,对于前者而言n=1时则代表序列x中的第一个元素值,对后者而言,na=0时代表序列中起始值,本质都是一个向量中的第一个元素,但是时间序列分析(尤其是频谱分析时)强调序列起始点,因此以不同的符号进行表示。
(2.1)构造重排序列:序列重排FFT是考虑包含序列中某点在所有循环移位后的分段傅里叶变换而产生的一种分析方法。
所述的序列重排是对包含着序列中的某个点的向量进行重新排序,进而进行循环抽取和排列,以x(0)为例,存在且只存在N个包含该点的N维向量:
上述式子中,向量右上角的T表示转置操作,x0到xN-1右下角的角标表示序列重排后的向量序号,将每个向量进行循环排序且把序列值x(0)移到序列中的第一位,则可得到另外的N个N维向量:
对准x(0)相加并取其平均值,则可得到重排后的数据向量,以xap表示如下:
(2.2)序列重排FFT函数求解:根据傅里叶变换的移位性质,以角标ka表示向量序号,则x’ka(n)的离散傅里叶变换X’ka(k)和x(n)的原始FFT频谱函数Xka(k)之间关系如下:
对上式的X′k(k)进行求和平均即为序列重排FFT函数:
解释上述过程如下:上面的分析结果表明:序列的序列重排FFT谱幅值为传统FFT谱幅值的平方,注意这里的平方关系是对所有谱线而言的,这意味着旁谱线相对于主谱线的比值也按照这种平方关系而衰减下去,从而使主谱线显得更为突出,因而序列重排后进行FFT的频率分析可以具有很好的抑制谱泄漏性能。对序号的分配进行解释:x(n)中的x(-N+1)实际上指的是长度为2N-1的x(n)中的第一个值,只是此时以x(0)为中心点,两侧序号分布是对称的。
(2.3)加窗序列频谱函数求解:对于长度为N的序列,若用窗函数Wf对单频复指数信号进行加窗FFT,以函数Fg为窗函数的频域形式,则得到的离散传统FFT的谱表达式为:
上式中,Δω=2π/N为频率分辨率单位,群延迟量为τ=(N-1)/2,根据FFT的LTI特性,对幅值为A、初相为θ、时间序列延迟为n0、角频率为ω*的复指数信号 再进行加窗FFT谱分析,则其离散谱表示为:
若对序列进行N个值和后N个值都执行加窗操作,则前、后窗函数均为Wf的双窗序列重排FFT的频谱函数表示为:
则根据序列重排FFT的线性特性,对信号进行同样的双窗序列重排FFT谱分析,其离散谱表示为:
/>
(2.4)构造相位差:以k*代表频率幅值谱中主频率谐波峰值的主谱线序号值,进而将传统相位谱值减去序列重排后的相位谱值/>则可得到如下的相位差/>表达式如下:
解释上述过程:上式表明,主谱线上的加窗FFT的相位谱和序列重排FFT相位谱的差值与频偏值dω=(ω*-k*Δω)成正比,其比例系数即为群延时τ。因此,只需通过相位之间的差值便可以实现频率和幅值的修正。
(2.5)频率和幅值修正:测出传统FFT和序列重排的两种谱分析在主谱线上的相位大小并作差,便可计算出频谱分析得到的频率与实际频率之间的偏移值,进而通过下式对信号真实频率的估计:
此处注意,若不只取主谱线的相位,也可以绘制相位差值的曲线。
分别取(式-13)和(式-15)的主谱线上的模值,有
将(式-15)进行进一步变换,即可得到频率幅值的修正值:
由(式-17)至(式-19)就可完成光强信号IP中3个谐波的频率、幅值估计,以频率值增序频排列,则经修正后的3个谐波的频率和幅值分别记为fp1、fp2、fp3和Ap1、Ap2、Ap3
在上面的计算中,窗序列Wf可以根据用户需求进行选取,常用的窗函数有汉宁窗、海明窗、矩形窗、三角窗、布莱克曼窗等,默认选择汉宁窗。
附图中图9所示的是频率和幅值修正所用到的曲线,在曲线中,需要对主频率进行定位,从而根据所定位到的相位值之差进行频率幅值修正。附图中图10所示的是一具体示例,在该例中,可以看到经过仿真软件模拟,FFT法和本方法的求解误差如下(由左至右分别是前表面、后表面、厚度变化信号):
其中,利用经典的FFT法求解的谐波频率和幅值误差如下:
FFT法求解的频率误差(Hz级别)=0.0040Hz,0.0080Hz,0.0053Hz;FFT法求解的幅值误差(-)=-0.1105,-0.0759,-0.0402。
而利用本方法求解得到的参数及误差如下:
序列重排FFT修正后的频率=7.1972Hz,2.3973Hz,9.6198Hz;真实频率与修正后频率的误差(10-3Hz级别)=0.0367×10-3Hz,0.0357×10-3Hz,-0.2639×10-3;振幅修正误差(-)=-0.0006,-0.0136,0.0430。
注意,上面结果中(-)代表无量纲。并且在附图中图9中也可以看出,相对于传统FFT,频谱泄漏效应得到的较好的抑制(因主要峰值附近的谱线高度较低)。上述误差分布和频谱分布结果充分验证了本方法的有效性。
对采样序列长度进行解释:在进行频域函数分析时,所采用的序列分析长度多为N(0,1,2,...,N-1),这是因为遵从传统的FFT计算的习惯而不致引起误解。对于序列重排傅里叶变换而言,所需要的总长度为2N-1(也即采样帧数,此处已经在步骤1中进行了强调),在进行序号标注时,2N-1个点的序号分别为(-N+1)~(N-1)的区间(边界含)。而经过步骤2中的步骤(2.1)所述的序列重排后,序列中的总点数由2N-1减少到N。在进行频率修正时,步骤2中是以信号中只有1个频率成分进行分析的,实则对不同频率组成的信号都是适用的,因本专利所针对的信号由3个谐波组成,因此主频率是3个、主谱线也是3个,逐一对3个主频率进行修正即可。而通过本专利所采用的补零方法以后(补零后的结果附图中图11所示),能够较好地对谐波信号的主频率进行探测,正如附图中图7所示。
至此,步骤2计算过程结束,转至步骤3。
步骤3:频谱折叠判定及相位求解:当(Ap2-Ap3)/Ap1大于折叠阈值EJS时,则判定不存在频谱折叠,若(Ap3-Ap2)/Ap1大于折叠阈值EJS则存在频谱折叠的问题。
当存在频谱折叠时,以FQS-1=fp1为厚度变化谐波的相位计算频率值、以FQS-2=fp2为前表面谐波的相位计算频率值、以FQS-3=fp3为后表面谐波的相位计算频率值;
当不存在频谱折叠时,FQS-1=fp1仍然为厚度变化谐波的相位计算频率值、以FQS-2=fp3为前表面谐波的相位计算频率值、以FQS-3=fp2为后表面谐波的相位计算频率值;
为各谐波的初始相位、Bp和Ap分别是各谐波对应的相位求解函数,进而干涉图中的初始相位的求解公式可以写为:
上式中U0表示干涉图中的光强,U0中包含着5N帧元素全为0的干涉图。各谐波对应的Bp和Ap分别通过下面式子确定:
上式中Re和Im分别表示求实部和求虚部,FQS-p、BP和AP中的右下角角标p=1、2、3时分别代表厚度变化信号、前表面信号和后表面信号。
求解出各谐波的包裹相位以后,对包裹相位进行解包裹、消倾斜以后便可根据波面高度和相位之间的关系进行各谐波面形的重建。
对上述过程进行解释,U0实际上为一矩阵,通过上面所描述的初始相位求解方法,可以快速地对初始相位的分布矩阵进行求解。U0(0)表示第一帧干涉图,也就是没有进行移相的干涉图,所求解的目标也就是该帧干涉图中的初始相位。应注意,在进行计算时,n=0直到7N-1,是将补零后的干涉图也计算在内的。从干涉图的角度进行补零,上面所提及的U0中包含着5N帧元素全为0的干涉图,也即包含着5N个与干涉图大小一致的全零矩阵。通过反正切求解出来的是被包裹在-π~π之间的,需要进行解包裹的操作,因被测面和参考面之间是存在着倾斜的,这种倾斜会体现在干涉信息和初始相位中,但是无法表现被侧面的面形信息,因此需要进行消倾斜操作,可理解为初始相位中包含着倾斜信息和面形信息。消倾斜过后,可根据面形高度与相位之间的线性关系,通过相位对面形高度进行重建。但是上述过程并非本专利的核心创新点,因此具体技术细节不再描述。
对上述过程进行解释:由于频谱折叠问题发生时,实际上是将超出频谱探测范围内(-fs/2~-fs/2,fs为采样频率)的谐波频率展现在-fs/2~-fs/2以内。对于多表面干涉问题而言,若实现对初始相位的准确计算,则只考虑FPS1~FPS3与fp1~fp3的对应问题即可,是一种高效的分析方法。但若同时要实现腔长的准确求解,则上述操作是不充分的,但本专利所解决的问题和创新点不集中在腔长计算中,因此在本专利中不再进行相关解释和分析。
附图中图12所示的是一测量示例,展示了模拟的三个叠加干涉波面以及通过本方法求解得到的3个重建波面及其包裹形式下的初始相位,并且结合附图中图13所示的该例对应的各表面重建误差,可以看到误差分布级别小于皮米级,这也证明了所提出的方法具有较高的测量精度。
至此,步骤3计算过程结束,测量完成。
实施例二:
在本实施例中,还需考虑到极端情况下的测量问题:
在步骤3中,若满足|(Ap2-Ap3)/Ap1|≤EJS则无法准确判断是否发生频谱折叠,将像素点移动到(xpp、ypp)后再次执行步骤1和步骤2;同时,折叠阈值可根据用户需求进行设定,一般应不小于0.2且不大于0.5。
上面对本发明实施例结合附图进行了说明,本发明全部流程图如附图中的图14所示,但本发明不限于上述实施例,还可以根据本发明的发明创造的目的做出多种变化,凡依据本发明技术方案的精神实质和原理下做的改变、修饰、替代、组合或简化,均应为等效的置换方式,只要符合本发明的发明目的,只要不背离本发明的技术原理和发明构思,都属于本发明的保护范围。

Claims (3)

1.一种基于频偏修正的多干涉谐波叠加信号频率及相位参数重建方法,其特征在于,包括以下步骤:
步骤1:通过波长调谐移相干涉仪采集2N-1帧干涉图,其中N应为不小于20的整数,选取干涉图中像素点(xp、yp)的光强记为Ip进行分析,设定近似熵阈值且对Ip进行近似熵计算,步骤1主要包括(1.1)向量序列分组、(1.2)向量距离量化计算、(1.3)量化向量相似度、(1.4)近似熵量化求解、(1.5)主频率搜索5个主要过程,具体计算步骤如下:
(1.1)向量序列分组:将光强Ip带入x(na)中,对于由NA个光强数据组成的时间序列x(na)=x(1),x(2),…,x(NA),近似熵的计算方法如下:
按序号组成一组维数为m的向量序列G,具体表示为其中向量序列G代表从第ia点开始的m个连续的x(na)中的元素值;
(1.2)向量距离量化计算:定义向量与/>的距离/>为:
上式中的max为求解最大值,||为计算绝对值;
(1.3)量化向量相似度:对于一组与/>计算/>小于等于误差半径ra的数量并将结果记为ci,进而计算εm,r
m和ra可以根据需求进行设定,总体误差变动范围用GA表示,数值上GA=2ra,但设定ra时应满足ra处于0.1FSD(x)和0.25FSD(x)之间的要求,FSD(x)为序列的x(na)的标准差;m一般设定为2;
(1.4)近似熵量化求解:增加维数到m+1,重复步骤(1.1)~(1.3),得到εm+1,r,并通过下式计算近似熵FApEn(m,r):
FApEn(m,r)=εm,rm+1,r (式-3)
(1.5)主频率搜索:对光强信号进行IP补零操作,补零长度为5N,补零后序列的总长度为7N-1且记为IP0,对IP0去均值后进行傅里叶变换并取单边幅值谱进行3个主要谐波峰值对应横坐标搜索,且将各个谐波峰值的横坐标记为主频率值;
至此,步骤1计算过程结束,如果近似熵大于预设熵值,则转至步骤2;
步骤2:将光强Ip带到x(n)中以后进行下一步处理,步骤2主要包括(2.1)构造重排序列、(2.2)序列重排FFT函数求解、(2.3)加窗序列频谱函数求解、(2.4)构造相位差、(2.5)频率和幅值修正五个具体步骤,流程如下:
以x(n)为单频复指数信号作为示例展示序列重排傅里叶变换的计算过程,若x(n)为一长度为N的序列,且其时域形式表示为:
上式中,信号的数字频率ω0表示为β倍频率间隔2π/N的形式,其中β可以是小数、θ为初始相位、N为总序列点数,j为虚数单位,则经过归一化以后,x(n)的原始FFT频谱函数可表示为:
上式中,k=0,1,…,N-1,进而可将上式展开,可以得到传统傅里叶变换以后长度为N的序列的频谱函数:
(2.1)构造重排序列:利用序列重排对包含着序列中的某个点的向量进行重新排序,进而进行循环抽取和排列,以x(0)为例,存在且只存在N个包含该点的N维向量:
上式中,x0到xN-1右下角的角标表示序列重排后的向量序号,将每个向量进行循环排序且把序列值x(0)移到序列中的第一位,则可得到另外的N个N维向量:
对准x(0)相加并取其平均值,则可得到重排后的数据向量,以xap表示如下:
(2.2)序列重排FFT函数求解:以角标ka表示向量序号,则x’ka(n)的离散傅里叶变换X’ka(k)和x(n)的原始FFT频谱函数Xka(k)之间关系如下:
对上式的X′k(k)进行求和平均即为序列重排FFT函数:
(2.3)加窗序列频谱函数求解:对于长度为N的序列,若用窗函数Wf对单频复指数信号进行加窗FFT,以函数Fh为窗函数的频域形式,则得到的离散传统FFT的谱表达式为:
上式中,Δω=2π/N为频率分辨率单位,群延迟量为τ=(N-1)/2,对幅值为A、初相为θ、时间序列延迟为n0、角频率为ω*的复指数信号 再进行加窗FFT谱分析,则其离散谱表示为:
对序列进行N个值和后N个值都执行加窗操作,则前、后窗函数均为Wf的双窗序列重排FFT的频谱函数表示为:
对信号进行同样的双窗序列重排FFT谱分析,其离散谱表示为:
(2.4)构造相位差:以k*代表频率幅值谱中主频率谐波峰值的主谱线序号值,进而将传统相位谱值减去序列重排后的相位谱值/>得到相位差/>表达式如下:
(2.5)频率和幅值修正:测出传统FFT和序列重排的两种谱分析在主谱线上的相位大小并作差,便可计算出频谱分析得到的频率与实际频率之间的偏移值,进而通过下式对信号真实频率的估计:
分别取(式-13)和(式-15)的主谱线上的模值,可得到:
将(式-15)进行进一步变换,即可得到频率幅值的修正值:
由(式-17)至(式-19)就可完成光强信号IP中3个谐波的频率、幅值估计,以频率值增序频排列,则经修正后的3个谐波的频率和幅值分别记为fp1、fp2、fp3和Ap1、Ap2、Ap3
至此,步骤2计算过程结束,转至步骤3;
步骤3:频谱折叠判定及相位求解:当(Ap2-Ap3)/Ap1大于折叠阈值EJS时,则判定不存在频谱折叠,若(Ap3-Ap2)/Ap1大于折叠阈值EJS则存在频谱折叠的问题;
当存在频谱折叠时,以FQS-1=fp1为厚度变化谐波的相位计算频率值、以FQS-2=fp2为前表面谐波的相位计算频率值、以FQS-3=fp3为后表面谐波的相位计算频率值;
当不存在频谱折叠时,FQS-1=fp1仍然为厚度变化谐波的相位计算频率值、以FQS-2=fp3为前表面谐波的相位计算频率值、以FQS-3=fp2为后表面谐波的相位计算频率值;
为各谐波的初始相位、Bp和Ap分别为各谐波对应的相位求解函数,进而干涉图中的初始相位的求解公式可以写为:
上式中U0表示干涉图中的光强,U0中包含着5N帧元素全为0的干涉图;各谐波对应的Bp和Ap分别通过下面式子确定:
上式中Re和Im分别表示求实部和求虚部,FQs-p、BP和AP中的右下角角标p=1、2、3时分别代表厚度变化信号、前表面信号和后表面信号;
求解出各谐波的包裹相位以后,对包裹相位进行解包裹、消倾斜以后便可根据波面高度和相位之间的关系进行各谐波面形的重建;
至此,步骤3计算过程结束,测量完成。
2.根据权利要求1中所述的一种基于频偏修正的多干涉谐波叠加信号频率及相位参数重建方法,其特征在于:步骤1中的近似熵判定的预设熵值可根据用户需求选择,默认选择0.3;步骤2中的窗序列Wf可以根据用户需求进行选取,默认选择汉宁窗。
3.根据权利要求1中所述的一种基于频偏修正的多干涉谐波叠加信号频率及相位参数重建方法,其特征在于:在步骤3中,若满足|(Ap2-Ap3)/Ap1|≤EJS则无法准确判断是否发生频谱折叠,将像素点移动到(xpp、ypp)后再次执行步骤1和步骤2;同时,折叠阈值可根据用户需求进行设定,一般应不小于0.2且不大于0.5。
CN202410133247.4A 2024-01-31 2024-01-31 基于频偏修正的多干涉谐波叠加信号频率及相位参数重建 Pending CN118089533A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202410133247.4A CN118089533A (zh) 2024-01-31 2024-01-31 基于频偏修正的多干涉谐波叠加信号频率及相位参数重建

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202410133247.4A CN118089533A (zh) 2024-01-31 2024-01-31 基于频偏修正的多干涉谐波叠加信号频率及相位参数重建

Publications (1)

Publication Number Publication Date
CN118089533A true CN118089533A (zh) 2024-05-28

Family

ID=91160524

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202410133247.4A Pending CN118089533A (zh) 2024-01-31 2024-01-31 基于频偏修正的多干涉谐波叠加信号频率及相位参数重建

Country Status (1)

Country Link
CN (1) CN118089533A (zh)

Similar Documents

Publication Publication Date Title
CN108872153B (zh) 基于非均匀傅里叶变换的平行平板光学均匀性的测量方法
EP0682771A1 (en) Method and apparatus for surface topography measurement by spatial-frequency analysis of interferograms
CN112097678B (zh) 基于频率盲估计的多表面面形测量方法
CN112923960B (zh) 用于校正非线性调谐效应的光纤参数测量装置
JP6207250B2 (ja) サンプルの性状を測定する方法及び装置
CN102221342A (zh) 一种时域多波长外差散斑干涉测量物体变形的方法
JPH05281048A (ja) 強度関係による空間波頭評価
CN102175332A (zh) 一种从含有移相误差的干涉图中恢复相位的方法
CN111397644A (zh) 一种用于光频域反射计的激光器非线性调谐效应补偿系统及补偿方法
CN110742584A (zh) 一种导管偏振敏感光学相干层析成像解调方法用偏振解算方法
CN102519597A (zh) 一种傅里叶变换光谱仪相位校正切趾方法
Huntley et al. Hyperspectral interferometry for single-shot absolute measurement of two-dimensional optical path distributions
CN105865371B (zh) 一种基于互相关计算的白光干涉显微轮廓复原方法
CN118089533A (zh) 基于频偏修正的多干涉谐波叠加信号频率及相位参数重建
CN111562088A (zh) 基于采样函数的平行平板光学参数的测量方法
JP6274555B2 (ja) 群遅延演算を用いたofdr方式光ファイバ計測方法及びそれを実施する装置
CN106441082B (zh) 一种相位恢复方法及装置
US20170138721A1 (en) Apparatus and method for processing the signal in master slave interferometry and apparatus and method for master slave optical coherence tomography with any number of sampled depths
CN114964033B (zh) 一种使用波长移相法测量超表面形貌分布的方法
Barajas et al. Towards an on-chip signal processing solution for the online calibration of SS-OCT systems
CN110500968B (zh) 基于稀疏傅里叶变换的数字莫尔干涉相位实时测量方法
CN109997010B (zh) 用于优化干涉仪的光学性能的方法及设备
Larkin Topics in multi-dimensional signal demodulation
Wei et al. Sub-pixel visualization of the envelope peak in a pulse train interferometer
Rajshekhar et al. Polynomial Wigner–Ville distribution-based method for direct phase derivative estimation from optical fringes

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination