一种基于混合基FFT的激光多普勒测速仪测速方法与流程

一种基于混合基fft的激光多普勒测速仪测速方法
技术领域
1.本发明涉及流体速度测量技术领域，更具体的说是涉及一种基于混合基fft的激光多普勒测速仪测速方法。


背景技术：

2.激光多普勒测速技术是利用激光器且基于光的多普勒效应测量并计算流体(固体)速度的一种技术，具有非接触测量、准确度高、空间分辨率高和动态响应快等独特优势，基于该技术设计的激光多普勒测速仪(laser doppler velocimetry,ldv)广泛应用于微粒测速、能源环保、临床医疗、生物医药、深潜勘探、卫星导航、工业生产等领域。
3.ldv主要由激光器、入射光学单元(分光系统)、收集光学单元(光接收系统)、多普勒信号处理系统、数据处理系统组成。其中多普勒信号处理系统负责对光接收系统收集的光电流信号进行预处理，提取多普勒信号后送至数据处理系统解算多普勒频率，最终计算得到速度。数据处理的过程通常由计算机完成，数据处理算法是ldv解算多普勒频率的关键。传统的数据处理算法有离散傅里叶变换(discrete fourier transform,dft)、基-2快速傅里叶变换(fast fourier transform,fft)与混合基fft的快速算法。这3种算法在ldv测速应用上都有各自的缺陷：dft将时域离散信号转为频域离散表示，计算时包含大量复数运算，在ldv实际测量过程中，大量复数运算会导致多普勒信号解算时间较长、测量实际速度产生一定的滞后，使现场测量的效率与实时性无法满足要求；基-2fft只能接收点数为2幂次的序列，造成了ldv测量时多普勒信号序列点数的局限性，无法满足解算任意点数信号序列的普适要求；混合基fft的快速算法是混合基fft的简化算法，在处理实际信号序列时仅将序列点数分解为两项乘积，虽然相比于dft在一定程度上提高了效率，但在序列点数较长以及序列点数为质数时运算效率与实时性仍旧不满足ldv现场测量的要求。
4.因此，如何提供一种能够有效化简复数运算，提高运算效率与实时性，有利于实现现场在线测量的测速方法是本领域技术人员亟需解决的问题。


技术实现要素：

5.有鉴于此，本发明提供了一种基于混合基fft的激光多普勒测速仪测速方法，可以有效化简复数运算，提高运算效率与实时性，有利于实现现场在线测量。
6.为了实现上述目的，本发明采用如下技术方案：
7.一种基于混合基fft的激光多普勒测速仪测速方法，包括：
8.按照预设采样间隔采集模拟光电信号，并进行预处理，得到离散多普勒信号序列；
9.将所述离散多普勒信号序列依次分解为较短的子序列；
10.构造dft迭代过程的三项变量因子矩阵；所述三项变量因子矩阵为：输入序列矩阵、系数旋转因子矩阵和dft旋转因子矩阵；
11.基于所述三项变量因子矩阵对所述子序列进行dft的逐层迭代计算，得到多普勒信号序列幅值频谱及其最大值，并求出频谱最大值的主频序号；
12.根据所述频谱最大值的主频序号依次计算多普勒频率和最终流体速度。
13.优选的，在上述一种基于混合基fft的激光多普勒测速仪测速方法中，第m层dft的表达式为：
[0014][0015]
其中，xm表示第m层的输出序列矩阵，表示第m层的dft旋转因子矩阵，x
m-1
表示第m层的输入序列矩阵，表示第m层的系数旋转因子矩阵。
[0016]
优选的，在上述一种基于混合基fft的激光多普勒测速仪测速方法中，所述将所述离散多普勒信号序列依次分解为较短的子序列，包括：
[0017]
将所述离散多普勒信号表示为n点序列x(n)形式，即x(0),x(1),
…
,x(n-1)；
[0018]
将n进行质因数分解，分解为n＝r1r2…rl
，其中，r1,r2,
…
,r
l
为l个从小到大依次排列的质数，用分解的l个质数作为混合基，将任意小于n的十进制非负整数表示为多基多进制形式。
[0019]
优选的，在上述一种基于混合基fft的激光多普勒测速仪测速方法中，第m层输入序列矩阵x
m-1
的构造过程为：
[0020]
设x
m-1
最初为n维向量，对n进行质因数分解后，所述离散多普勒信号序列由n维向量变为x
m-1
(k
l-1
,k
l-2
,
…
,k
l-m 1
,n
l-m
,
…
,n0)的l维矩阵，运用降维规则将l维矩阵降维为二维矩阵；
[0021]
对降维后的二维矩阵进行dft的迭代计算，输出二维矩阵xm；
[0022]
对输出的二维矩阵xm按照降维时顺序的逆序重新升维成新的l维矩阵xm(k
l-1
,k
l-2
,
…
,k
l-m
,n
l-m-1
,
…
,n0)，再按基顺序的逆序重新组合元素，构成n维向量，获得下一层的n维输入向量xm。
[0023]
优选的，在上述一种基于混合基fft的激光多普勒测速仪测速方法中，所述降维规则为：
[0024]
将l维矩阵x
m-1
的第m维元素x
m-1
(n
l-m
),n
l-m
＝0,1,
…
,r
m-1作为输入序列矩阵的第一列，rm表示第m个基；
[0025]
将剩下维数中的元素按基的逆序，即r
lrl-1
…rm 1rm-1
…
r1的顺序依次添加到矩阵的列；
[0026]
最终得到一个降维后的二维矩阵x
m-1
(n
l-m
,(k
l-1
…kl-m 1nl-m-1
…
n1n0))。
[0027]
优选的，在上述一种基于混合基fft的激光多普勒测速仪测速方法中，第m层的系数旋转因子矩阵由多个旋转因子子项的乘积组成，其表达式为：
[0028]
[0029]
其中，t＝1,2,
…
,m-1，表示第m层系数旋转因子子项的个数，为旋转因子子项；r
para
＝r
m 1rm 2
…rl
；k
para
为矩阵列数，0≤k
para
≤n/r
m-1；为k
para
的部分倒位序。
[0030]
优选的，在上述一种基于混合基fft的激光多普勒测速仪测速方法中，的取值过程为：
[0031]
将剩余基表示为原基r1,r2,
…
,rm,
…
,r
l
去掉rm形成的基为r1,r2,
…
,r
m-1
,r
m 1
,
…
,r
l
；
[0032]
确定第m层的输入为x
m-1
(k
l-1
,k
l-2
,
…
,k
l-m 1
,n
l-m
,
…
,n0)，则剩余基正位序排列为：
[0033]kpara
＝k
l-1
(r2…rm-1rm 1
…rl
)
…
k
l-m 1
(r
m 1
…rl
) n
l-m-1
(r
m 2
…rl
)
…
n0；
[0034]
剩余基倒位序排列为：
[0035]kl-1

…
k
l-m 1
(r1r2…rm-2
) n
l-m-1
(r1r2…rm-1
)
…
n0(r1r2…rm-1rm 1
…rl-1
)；
[0036]
的值为剩余基倒位序排列中含k项的和，计算公式为：
[0037][0038]
优选的，在上述一种基于混合基fft的激光多普勒测速仪测速方法中，第m层dft旋转因子矩阵的表达式为：
[0039][0040]
优选的，在上述一种基于混合基fft的激光多普勒测速仪测速方法中，最终流体速度的计算公式为：
[0041]
v＝fd·
λ/2sinθ；
[0042]
其中，fd＝kd·fs
/n，fs为采样频率；λ为激光多普勒测速仪中激光的波长；θ为入射激光夹角；kd为频谱最大值的主频序号。
[0043]
优选的，在上述一种基于混合基fft的激光多普勒测速仪测速方法中，
[0044]
经由上述的技术方案可知，与现有技术相比，本发明公开提供了一种基于混合基fft的激光多普勒测速仪测速方法，在ldv测量时，将基于库利-图基思想的混合基方法与离散傅里叶变换(discrete fourier transform,dft)进行结合，对多普勒信号序列进行频谱分析并解算多普勒频率，最终实现ldv对流体速度的实时在线测量。该迭代的过程可以有效化简复数运算，提高运算效率与实时性，有利于实现现场在线测量，改善了ldv测量中dft的运算效率、打破了基-2快速傅里叶变换(fast fourier transform,fft)对输入信号序列长度的局限性，并使数据处理系统能够测量过程中高效处理任意长度的信号序列输入，具有普适性与实用性，也保证了ldv对流体速度的精确测量。
附图说明
[0045]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现
有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。
[0046]
图1为本发明提供的基于混合基fft的激光多普勒测速仪测速方法的流程示意图。
具体实施方式
[0047]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0048]
如图1所示，本发明实施例公开了一种基于混合基fft的激光多普勒测速仪测速方法，包括以下步骤：
[0049]
s1、按照预设采样间隔采集模拟光电信号，并进行预处理，得到离散多普勒信号序列；
[0050]
s2、将离散多普勒信号序列依次分解为较短的子序列；
[0051]
s3、构造dft迭代过程的三项变量因子矩阵；三项变量因子矩阵为：输入序列矩阵、系数旋转因子矩阵和dft旋转因子矩阵；
[0052]
s4、基于三项变量因子矩阵对子序列进行dft的逐层迭代计算，得到多普勒信号序列幅值频谱及其最大值，并求出频谱最大值的主频序号；
[0053]
s5、根据频谱最大值的主频序号依次计算多普勒频率和最终流体速度。
[0054]
下面对上述各步骤进行进一步描述。
[0055]
s1、设置采样频率fs，则采样间隔ts＝1/fs，信号处理系统依据采样间隔对光接收系统收集的模拟光电信号进行采样并进行预处理，数据处理系统接收预处理后的离散多普勒信号n点序列x(n):x(0),x(1),
…
,x(n-1)。
[0056]
s2、将n进行质因数分解，分解为n＝r1r2…rl
，其中，r1,r2,
…
,r
l
为l个从小到大依次排列的质数，称为基，可以用这l个质数作为混合基将数表示为多基多进制形式；也就是说，可以将任意小于n的十进制非负整数表示为多基多进制形式。
[0057]
s3、在构造迭代过程的三项变量因子矩阵前，需先确定dft的表达式，dft表达式如式(1)所示：
[0058][0059]
其中，n表示变换前序列的序号；k表示序列经dft后的主频序号；为旋转因子，是由数据处理系统构造的复数因子矩阵，用于处理多普勒信号序列；x(k)表示序列经dft后的序列频域结果，是复数。
[0060]
信号处理系统同时将十进制表示的序号n与k表示为由r1,r2,
…
,r
l
为混合基的多基多进制形式，将n表示为正位序、k表示为倒位序。其中ni＝r
l-i
,i＝0,1,
…
,l-1，kj＝r
j 1
,j＝0,1,
…
,l-1均为各个数位的系数。
[0061][0062]
信号处理系统将式(2)代入式(1)：
[0063][0064]
根据多项式乘法以及旋转因子的性质系统对式(3)中l2项幂项的旋转因子进行简化，将(l
2-l)/2项简化为0，剩余未简化幂项结果如表1所示。
[0065]
表1旋转因子剩余幂项
[0066][0067]
其中，“/”表示该幂项由于旋转因子的性质被消去，旋转因子的值为1。
[0068]
根据式(3)将dft分为l层并逐层求和及迭代，将表1的化简结果代入式(3)可得第一层dft公式。
[0069][0070]
第二层dft公式如下。
[0071][0072]
第三层dft公式如下：
[0073][0074]
以此类推，为统一变量，记x(n)为x0，第m层dft公式如下(1≤m≤l)。
[0075][0076]
根据式(4)～式(7)，构造dft迭代过程中的3项变量因子分别为：第m层的输入序列x
m-1
、系数旋转因子dft旋转因子具体如下：
[0077][0078]
将式(7)化为矩阵形式。
[0079][0080]
其中，x
m-1
与为同型矩阵，它们的乘法为矩阵各元素对应相乘。为同型矩阵，它们的乘法为矩阵各元素对应相乘。与的乘法为矩阵乘法运算，的列数等于的行数。
[0081]
下面，对输入序列矩阵、系数旋转因子矩阵和dft旋转因子矩阵三个变量因子的构造过程进行详细说明。
[0082]
1)第m层的输入序列x
m-1
矩阵的构造。
[0083]
x
m-1
是第m层的输入序列，为n维向量，当对n进行质因数分解后，根据dft同址运算的原则，x
m-1
变为x
m-1
(k
l-1
,k
l-2
,
…
,k
l-m 1
,n
l-m
,
…
,n0)的l维矩阵，运用如下降维规则将其变为二维矩阵：
[0084]
①
将l维矩阵x
m-1
的第m维元素x
m-1
(n
l-m
),n
l-m
＝0,1,
…
,r
m-1作为输入序列矩阵的第一列；
[0085]
②
将剩下维数中的元素按基的逆序即r
lrl-1
…rm 1rm-1
…
r1的顺序依次添加到矩阵的列；
[0086]
③
最终得到一个降维后的二维矩阵x
m-1
(n
l-m
,(k
l-1
…kl-m 1nl-m-1
…
n1n0))；
[0087]
用降维后矩阵按照式(9)进行运算，得到输出二维矩阵xm，对其按照降维时顺序的逆序重新升维成l维矩阵xm(k
l-1
,k
l-2
,
…
,k
l-m
,n
l-m-1
,
…
,n0)，再按基顺序的逆序重新排列向量元素即可获得下一层的n维输入向量xm。
[0088]
2)第m层的系数旋转因子矩阵的构造
[0089]
根据式(8)，系数旋转因子由许多旋转因子的乘积构成，如下式。
[0090][0091]
的取值过程如下：
[0092]
①
剩余基可以表示为原基r1,r2,
…
,rm,
…
,r
l
去掉rm后形成的基，形式如下r1,r2,
…
,r
m-1
,r
m 1
,
…
,r
l
。
[0093]
②
设矩阵列数为k
para
,0≤k
para
≤n/r
m-1。
[0094]
第m层的输入为x
m-1
(k
l-1
,k
l-2
,
…
,k
l-m 1
,n
l-m
,
…
,n0)，则剩余基正位序排
[0095]
列为：
[0096]kpara
＝k
l-1
(r2…rm-1rm 1
…rl
)
…
k
l-m 1
(r
m 1
…rl
) n
l-m-1
(r
m 2
…rl
)
…
n0ꢀꢀꢀ
(11)。
[0097]
剩余基倒位序排列为：
[0098]kl-1

…
k
l-m 1
(r1r2…rm-2
) n
l-m-1
(r1r2…rm-1
)
…
n0(r1r2…rm-1rm 1
…rl-1
)
ꢀꢀꢀ
(12)。
[0099]
③
的值为剩余基倒位序排列中含k项的和。
[0100]
即即为k
para
的部分倒位序。
[0101]
3)第m层的dft旋转因子矩阵的构造：
[0102]
根据式(8)可知，可直接构造为如下对称矩阵。
[0103][0104]
s4、构造完迭代过程中的3项变量因子后，根据式(9)进行l层迭代运算后得到x(k)序列，最终由信号处理系统输出多普勒信号序列经混合基fft的频谱，计算多普勒信号序列幅值频谱的最大值max{|x(k)|}，得到频谱最大值的主频序号kd。
[0105]
s5、根据频率计算公式fd＝kd·fs
/n计算多普勒频率fd，并根据流速换算公式计算
出最终流体速度v＝fd·
λ/2sinθ。其中λ为ldv中激光的波长，θ为入射激光夹角，均为设定值，fs表示采样频率。
[0106]
在其他实施例中，还包括：将多普勒频率fd以及对应速度v实时显示在屏幕上，并绘制图像，达到现场实时测量的目的。
[0107]
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的装置而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
[0108]
对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。