Fourier分析

Fourier分析（共7篇）

Fourier分析 篇1

青黛是一种常用中药,具有清热解毒、凉血散肿的功能[1],由青黛、生石膏、冰片、大蚯蚓等组成的青黛散在外治流行性腮腺炎方面疗效显著[2]。X衍射图谱分析方法在中药鉴定中具有广阔应用前景[3],但大部分都是对单味药的研究,用于复方鉴别的报道很少。本文采用X射线衍射的方法对复方青黛散及其组成药物进行鉴别分析,为制定该制剂的质量控制方法提供了新的思路。

1、仪器与试药

1.1 仪器与条件:

日本理学公司D/max-3c改进型X射线衍射仪,Cu辐射,石墨单色器,管压35kv,管流40m A,扫描速度16o/min,步长0.02o,2θ扫描范围5o~60o。Sartorius BP61S型电子分析天平。

1.2 药品:

青黛、生石膏、冰片、大蚯蚓药材购于西安藻露堂药店,经西安市药品检定所谢志民研究员鉴定为正品药材。

2、X-射线衍射法定性分析

2.1 供试品制备:

取青黛、生石膏、冰片、分别研末成细粉,过100目筛制成细粉。青黛散的制备按青黛、生石膏、冰片、大蚯蚓比例为5:5:2:2配制,混合均匀并研细后备用。

2.2 X-射线衍射图谱

青黛、生石膏、冰片及混合而成的青黛散的特征X-衍射图谱分别见于图1~4。其中标*、△、●分别为青黛、生石膏、冰片的特征衍射数据。特征数据的选取参考国际衍射数据中心(ICDD)提供的标准数据对药材衍射数据进行取舍,以排除杂质干扰。

3、图谱数据分析

青黛的X衍射图谱相应的衍射特征标记峰为:3.060/100*,2.514/15*,2.294/21*,2.104/20*,1.921/22*,1.884/23*,1.609/10。

冰片的X衍射图谱相应的衍射特征标记峰为:11.442/10,5.764/100●,5.007/37●。

生石膏的X衍射图谱相应的衍射特征标记峰为:4.304/87△,3.815/46△,3.074/100△,2.886/18△,2.692/18△,2.599/7△,2.462/5△,2.222/16△,2.085/10△,1.901/18△,1.816/17△,1.623/1 0△。

青黛散的X衍射图谱相应的衍射特征标记峰为:8.499/35,7.622/89△,5786/52●,4.962/20●,4.259/89△,3792/30△,3.329/36△,3.173/19,3050/100*,2.866/44△,2.778/19△,2679/39△,2.499/23△,2.444/21,2285/19*,2.213/22,2.075/35△,1.901/30*,1.877/31*。

4、讨论

本文获取了青黛和冰片中药材的X射线衍射图谱,经与标准X衍射图谱比较显示青黛所含晶体为Ca CO3;冰片虽为有机物,也有衍射特征较强的参数。通过分析青黛散的X衍射图谱,清楚地显示了青黛和冰片各自的特征衍射数据。显示利用X衍射用于分析矿物药是简便可行的,并且在分析矿物药组成的复方制剂时也有巨大潜力,对其进一步完善后有可能实现对药物的定量分析[4]。因此,X射线衍射分析法完全可以发展成为一种新的药物鉴定和质量控制的方法。

摘要：目的:用X射线衍射Fourier(傅立叶)图谱及其特征衍射峰值来鉴别青黛散及其组成的药物。方法:采用X衍射的方法,对青黛散及其各组分的X衍射Fourier图谱进行分析。结果:获得了青黛散中各组分及组方后的X射线衍射Fourier图谱的几何拓扑图形及其衍射特征峰值。结论:X-衍射Fourier图谱分析法可以用于青黛散类制剂的鉴定和识别。

关键词：青黛散,X-衍射,图谱分析,质量标准

参考文献

[1]华海清.青黛的古今应用[J].中国中药杂志.1992;17(2):117-120

[2]赵新洪,李永全,等.青黛散局部涂敷加静滴清开灵治疗流行腮腺炎疗效观察[J].中华临床医药杂志(广州),2003;(65):53-54

[3]朱志峰,王树春,等.中药蛤蚧的X衍射Fourier谱分析[J].中草药.2001;32(10):932-933

[4]陈海明,赵斌等.柏子养心丸的X衍射Fourier谱分析[J].中成药.2001;23(11):792-793

Fourier分析 篇2

随着便携式电子产品小型化、多功能化和环境保护的要求, 片上系统 (System on Chip, So C) 等高密度三维立体 (3D) 堆叠封装技术得到了迅速发展, 这种技术在减小芯片尺寸和实现紧密集成度等方面可以满足目前电子行业的要求[1,2]。

现代电子技术迅猛发展, 要求电子整机朝着短、小、轻、薄及高可靠、高速、高性能和低成本的方向发展[3]。热失效是电子设备的主要失效形式, 其失效率随温度增加将呈指数增长趋势[4]。在设计阶段, 热分析非常关键, 通过其可以发现电子设备的热缺陷[5,6]。

当前电子设备的集成度不断提高, 器件的特征尺寸也愈来愈小。当尺度 (空间尺度和时间尺度) 微细化后, 将会出现传热的尺度效应, 此时的传热情况与常规尺度下不同。对于多芯片组件, 更多的使用经典Fourier理论去分析解决问题, 但尺度微细化后的热量传播包含了不同于扩散的波动传播机制, 这是以经典Fourier定律为基础而建立起来的瞬态传导理论不能合理解释的。但目前关于电子封装的Non-Fourier效应研究较少, 相关的数据资料也比较有限。本文尝试以三维多芯片组件为模型, 分别用Fourier和Non-Fourier分析方法研究其热效应问题, 并对实验结果进行了分析对比。

1 多芯片组件导热模型

1.1 多芯片组件模型

本文对三维多芯片组件进行分析, 在直角坐标系中, 建立如图1 所示的三维多芯片组件模型。基板模型呈长方体状, 其几何尺寸为:xa× yb× zc= 70 × 70 × 5 mm3Z = 0 表面标有编号1~6 的位置布置有6 个芯片, 在此忽略芯片的实际尺寸, 只在基板表面有芯片的位置考虑有若干热源的情况。在热源位置内, 热源的热耗散均匀分布, 与坐标位置和时间均无关。

1.2 热传导问题的Fourier分析

在《热的分析理论》中[7], 法国科学家约瑟夫·傅里叶 (Joseph Fourier) 结合实验观察, 阐述了把温度梯度和热流联系起来的基本定律:对于均匀的各向同性固体, 傅里叶定律可表示为:

式中:热流密度q (r, t) 表示单位等温面上, 在温度降低方向上单位时间内的热流量 (单位:W·m2) , k为材料的导热系数 (单位:W/ (m·℃) ) , T (r, t) 为温度分布 (单位:℃) 。由于q (r, t) 指向温度降低方向, 式 (1) 中的负号就使热流量成为一个正的量。在直角坐标系中, 式 (1) 被改写为:

式中i, j, k分别为沿x, y, z方向的单位向量。

对一个很小的控制体V建立热量平衡方程, 可描述如下[8]:

对于上式中的每一项, 以其各自表达式表示, 可得式 (3) :

对式 (3) 进行转化, 可得静止均匀物体内含有热源的各向同性物质的热传导方程, 即:

1.3 热传导问题的Non-Fourier分析

随着工程技术的快速发展, 热传导领域的瞬态过程日益受到人们的重视[9,10]。对于稳定温度分布的无限大介质, 突然出现的热扰动会对介质的稳定状态造成破坏, 历经松弛过程之后, 介质将处于新的稳定状态。在温度场的重建过程中, 松弛时间以 τ0表示, 热量传播速度以Ch表示, 从热扰动的所在点算起, 热量传播深度为δh= Chτ0。由于热扩散率的量纲m2/s= (m/s) ·m, 因此可得δh= α Ch, 则热量传播速度的表达式 (以松弛时间和热扩散率表示) 如下:

修正经典的Fourier定律, 可得:

其中, 为在温度梯度截面上热流密度对时间的变化率, 为松弛时间间隔内截面上热流密度的改变量。

当热量传播速度为有限值时, 整理后的热传导微分方程如下:

若内热源Q与时间无关, 式 (7) 可简化为:

2 Non-Fourier差分方程的求解

相比于Fourier导热微分方程, Non-Fourier导热微分方程的求解更为困难。对于Non-Fourier导热微分方程, 仍然使用交替方向隐式法, 得到如下差分格式:

其中, τ0为松弛时间;ΔT为时间步长;θ 的下标j, i, i分别表示x, y, z三个方向上的网格节点号, θ 的上标n表示时间层;Δx, Δy, Δz分别为x, y, z三个方向上的空间步长。

使用中心差商近似替代Non-Fourier导热模型的边界条件, 以式 (9) 为例, 计算绝热侧面 (x = xa) 上的温度时, 差分格式如下:

通过同样的方法, 即得其他边界上相应的差分格式。利用初始条件, 联立边界条件上的差分近似式和差分方程式 (9) ~ (11) , 就可进行数值求解。

3 温度场计算及结果

首先建立Fourier和Non-Fourier导热微分方程相应的差分格式, 然后通过求解得到的关于网格节点温度的方程组, 即可求得所有节点温度值。

图2 分别是Fourier和Non-Fourier分析中, 顶面最高温度随时间变化的温度曲线图。在开始的100 s内, Fourier和Non-Fourier分析对应的顶面最高温度上升趋势相当。100 s以后, Fourier分析对应的温度上升速度逐步放缓, 并迅速进入到温度场稳定状态, 在600 s后稳定在大约70 ℃。Non-Fourier分析对应的温度曲线则比较陡峭, 在100 s之后温度上升的速度仍然比较快, 在1 200 s左右才开始逐渐进入温度场的稳定状态, 稳定之后温度大约为105 ℃。在整个过程中, Non-Fourier分析对应的顶面最高温度要高于Fourier分析, 两种分析相应的温度差也不断增大。Fourier分析可以看成Non-Fourier分析的特殊情况, 即松弛时间为零的情况, 在相同条件下, 其顶面最高温度的相应曲线落在图2 两条曲线之间。

图3 为顶面温度场的相应最大温度差随时间变化的曲线图。短暂的时间过后, 两种分析的最大温度差都迅速稳定下来, 远快于整个温度场达到稳定状态所用的时间。相比于Fourier分析, Non-Fourier分析的最大温度差达到稳定所需的时间要略长一些。

图2 和图3 出现的这些差别, 是由于Fourier分析以热量传播速度无限大为前提, 同时热量被整个物体的全部质点不同程度地吸收;而Non-Fourier分析考虑了热量的传播速度, 进入某一截面的热流一部分以热传导的方式向物体内部传导, 而另一部分用在了该截面分子松弛过程所导致的热量改变。

应用有限元热分析软件Icepak建立三维多芯片组件模型, 如图4 所示。

由图4 可知, 软件仿真和通过计算所得的温度场分布大致相同, 但Non-Fourier模型更加接近。

分析在某些关键点和路径上的温度分布状况, 在基板上选取下述关键点和路径。

路径1:y=21, z=0, x∈ (0, 70)

路径2:y=35, z=0, x∈ (0, 70)

路径3:y=56, z=0, x∈ (0, 70)

路径4:x=42, y=42, z∈ (0, 5)

路径5:x=21, y=21, z∈ (0, 5)

路径6:x=14, y=56, z∈ (0, 5)

路径7:x= y, z = 0, x ∈ (0, 70) , y ∈ (0, 70)

路径8:x= 70 - y, z = 0, x ∈ (0, 70) , y ∈ (0, 70)

图5 分别为在Fourier和Non-Fourier分析两种情况下, 沿x方向路径1、2、3 的温度曲线图。由图5 可知, 除温度值不同, 及热源6 位置与其附近区域的温度值在曲线上略有不同外, 两种分析情况中的路径曲线走势基本是一致的。

图6 分别为在Fourier分析和Non-Fourier分析两种情况下, 沿z方向路径4、5、6 的温度曲线图。由图6 可知, 这三条路径上的温度从顶面到底面逐渐降低, 且温度变化不大, 变化规律基本一致。这是因为沿z方向的尺寸很小, 并且仅有底面对外传热。

图7 的路径7、8 是顶面上的两条对角线。路径8 的温度曲线左右对称, 越靠近路径8 的中点, 温度越高。整个温度场的分布关于路径7 对称。这是因为基板表面布置的热源在也是以路径7 为对称轴对称分布的, 而且发热量在对称轴附近是最大的。因而热源的分布状况对温度场的影响很大, 所以基板上芯片与其他热源的位置分布就显得尤为重要, 要尽量分散热源, 避免过于集中而导致局部温度过高。

路径贯穿热源和热源所在的位置, 分析中路径8 的温度曲线中部比Non-Fourier分析的曲线变化更为剧烈, Fourier分析中曲线的向上凸起程度不及Non-Fourier分析。同时, 与热源位置处的温度相比, 热源1 和6 之间相同区域的温度都比较低, 而Non-Fourier分析中没有出现这一情况。

由此可见, Non-Fourier分析中的热耦合现象明显强于Fourier分析。

对温度场进行分析可得, 使用Fourier分析和Non-Fourier分析得到的结果存在着很多差异。Non-Fourier导热模型的温度上升速度快于Fourier分析, 温度曲线比Fourier导热模型的更为陡峭, 进入稳定状态的时间远迟于Fourier分析且进入稳定状态后温度值较高, Non-Fourier模型的温度场最大温差要大于Fourier模型, 热耦合现象也较强。

4 结语

本文分别建立了Fourier分析和Non-Fourier分析导热微分方程的有限差分格式, 通过求解相应代数方程组, 得到两种分析情况下的温度场分布;分析比较了使用Fourier和Non-Fourier分析对基板温度场分布所产生的不同影响, 并在热分析软件Icepak中对模型进行了验证。由于对Non-Fourier分析的数学处理相对比较困难, 本文研究的是相对简单的基板模型, 没有考虑芯片与基板的连接、芯片的体积、以及基板上的布线结构等问题, 与多芯片组件结构的实际情况有一定的区别, 还有待进一步改进。

摘要：随着电子产品的尺寸和重量日趋变小, 经典Fourier理论已经不能很好地解决实际问题。为了使理论结果尽可能准确地反映实际情况, 需要用到Non-Fourier导热模型。首先, 建立了三维多芯片组件的Fourier和Non-Fourier导热模型;其次, 采用有限差分方法求解相应模型的传热方程, 得到了温度分布和温度响应;最后, 使用有限元热分析软件ICEPAK建立了相应的热分析模型并进行计算。实验结果表明, 与经典的Fourier模型相比, Non-Fourier模型更加接近实际温度, 且温度场进入稳态的时间较长, 热耦合的现象也更强。

关键词：多芯片组件,Fourier,Non-Fourier,有限差分法

参考文献

[1]LIU Shu-hua, CAO Li-qiang, LI Jun, et al.Study of mixed-signal crosstalk in 3-D package[C]//Proceedings of 2010 11th International Conference on Electronic Packaging Technology&High Density.Xi’an, China:ICEPT-HDP, 2010:559-562.

[2]RICKERT P, KRENIK W.Cell phone integration SiP, SoC, and PoP[J].IEEE Design&Test of Computers, 2006, 23 (3) :188-195.

[3]KUO W S, TZENG Y L, CHEN E, et al.POP package (Cavity BGA) warpage improvement and stress characteristic analyses[C]//Proceedings of 2007 International Microsystems, Packaging, Assembly and Circuit Technology.Orlando, FL:IMPACT, 2007:342-345.

[4]于慈远, 于湘珍, 杨为民.电子设备热分析/热设计/热测试技术初步研究[J].微电子学, 2000, 30 (5) :334-337.

[5]WEISS J, LANGHORST F.Thermal analysis of electronic packaging[J].Australian Electronics Engineering, 1988, 21 (8) :62-64.

[6]SWAGER A W.Thermal images provide reliability clues[J].END, 1990, (1) :47-58.

[7]王进春.热传导中的非傅里叶效应研究[J].中国电力教育, 2009 (3) :144-145.

[8]BERGLES A E.Heat transfer in electronic and microelectronic equipment[M].Troy, NY:Hemisphere Publishing Corporation, 1990.

[9]蒋方明, 刘登瀛.非傅里叶导热的最新研究进展[J].力学进展, 2002, 32 (1) :128-140.

Fourier分析 篇3

傅里叶变换轮廓术 (FTP) 具有单帧获取、全场分析和高分辨率等优点。随着计算机技术硬件和软件的发展, 以及图像获取设备分辨率的提高, 傅里叶变换轮廓术倍受人们的关注, 已成为三维传感中最重要和最活跃的研究领域之一。

利用光栅投影法原理进行三维物体轮廓测量, 关键是要选取合适的参数, 保证系统测量精度和测量范围, 并保证傅里叶频谱的完全分离。在此重点讨论了利用等效波长的概念, 对傅里叶变换轮廓术 (FTP) 中的参数进行选择和优化。从理论和实践两方面进行了深入研究。

1FTP中各参数的优化与选择理论分析

傅里叶变换轮廓术测量光路图如图1所示, G表示正弦光栅, P表示投影系统, C表示CCD摄像机, S表示监视器。在摄像机和投影仪光心的连线与参考平面平行以及摄像机光轴与投影仪光轴平行且L≫h (x, y) 的情况下, 被测物体的高度与相位差之间的映射[1,2,3,4]关系为:

$h(x,y)⧋\frac{L\text{Δ}\Phi (x,y)}{2\text{π}f{}_{0}D}=\frac{\text{Δ}\phi (x,y)}{2\text{π}f{}_0}\text{t}\text{g}\theta (1)$

式中, L为CCD摄像机入瞳到参考平面的距离;D是投影系统出瞳到摄像机入瞳的距离;tg θ=D/L;θ为投影仪光轴与摄像机光轴之间的夹角。由式 (1) 可知, 物体重建的精度取决于位相Δφ (x, y) 的测量精度, 光栅像的基频f0, 以及投影仪光轴与摄像机光轴之间的夹角θ。定义等效波长:λc=p0/tg θ, 用以表征系统的测量精度, 等效波长越长, 系统的测量精度越低。

$\text{d}h/\text{d}\phi =1/2\text{π}f{}_0\text{t}\text{g}\theta =\lambda {}_c/2\text{π}(2)$

用λc和 (dh/dφ) 表征系统的测量精度是一致的。定义局部空间频率为:

$\begin{array}{l}f{}_n=\frac{1}{2\text{π}}\cdot \frac{\partial }{\partial x}[2\text{π}nf{}_{0}x+n\phi (x,y)]\\ =nf{}_0+\frac{n}{2\text{π}}\cdot \frac{\partial \Phi (x,y)}{\partial x}(3)\end{array}$

引入物面倾面α (x, y) =tg-1 (∂h (x, y) /∂x) , 将式 (3) 代入后有:

$f{}_n=nf{}_0(1+\text{t}\text{g}\theta \cdot \text{t}\text{g}\alpha )(4)$

由式 (4) 可知, 当减小夹角θ时, tg θ亦减小, fn得到压缩, 而等效波长λc却增大, 所以频谱的压缩是以系统测量精度的降低为代价的。

压缩夹角θ, 使基频同零频和高频完全分离, 须满足条件:

$\begin{array}{l}(f{}_1){}_{\min }\ge f{}_b(5)\\ (f{}_2){}_{\min }\ge (f{}_1){}_{\max }(6)\end{array}$

式中, fb为零频最大值。解得各次谱完全分离时夹角θ所必须满足的条件为:

$\begin{array}{l}\text{t}\text{g}\theta \le \frac{f{}_0-f{}_b}{f{}_0}\cdot \frac{1}{\text{t}\text{g}\alpha {}_{\max }}(7)\\ \text{t}\text{g}\theta \le \frac{1}{3}\cdot \frac{1}{\text{t}\text{g}\alpha {}_{\max }}(8)\end{array}$

必须对 (f0-fb) /f0与1/3进行比较, 以决定θ的选择。不考虑零频的影响, 假设各次频谱完全分离时的夹角为θc, 则θc由下式确定:

$\theta {}_c=\text{t}\text{g}{}^{-1}(1/3\text{t}\text{g}\alpha {}_{\max })(9)$

当θ<θc时, 可适当增大θ, 以提高测量精度;当θ>θc时, 可减小θ, 以满足频谱完全分离的条件。随着θ的减小, 等效波长增大, 在相同的位相测量误差条件下, 系统的整个测量精度下降。为保证系统测量的精度, 由等效波长的定义可知, 在减小夹角的同时, 提高投影光栅的空间频率f0, 可保证λc不变。由以上分析可得:压缩夹角, 提高投影光栅的空间频率, 可在频域中避免各次频谱的混叠, 并保证系统的测量精度。

2参数选择和优化的计算机模拟

为了证明结论的正确, 用计算机模拟了整个程序, 物体总长度为N=256像素, 投影光栅采用正弦光栅, 经过多次选择和处理后, 模拟的光学系统最佳参数如下:f0=1/16 (条纹/像素) , θ=25°, 滤波窗采用汉宁窗 (采用汉宁窗作为滤波窗, 可以更有效地抑制噪声, 减小误差) , 窗口大小选择60。

以下是模拟结果:图2为被测物体的高度分布图;图3为物体的条纹灰度图, 汉宁窗在x, y方向上的分布情况如图4所示。模拟物体折叠时的灰度图如图5所示, 展开后的相位图如图6所示。

在此, 对经参数选择和优化后的图像与未经参数选择和优化的图像进行比较。通过观察分析可以直观地看出, 图8与模拟图像原图2最为相近, 而图7相对误差很大, 噪声显著。

3人脸面具模型的实验验证

为了验证上述问题, 对一个人脸面具模型进行了实际测量。其投影图像由计算机产生, 并经液晶投影仪 (LCD ) 投影到被测物体表面。系统结构参数为:

$f{}_0=1/16(\text{条}\text{纹}/\text{像}\text{素})L{}_0=800\text{m}\text{m},d=350\text{m}\text{m}$

在此, 对经过参数选择和优化的物体三维相位网络图10与图11进行了比较, 可以明显看出: 经过参数选择和优化后的图像, 更有效地避免了零频分量和基频分量之间的频谱混叠, 提高了FTP 的测量精度, 图10更加精确清晰, 误差得到了显著降低。

为了更好地观察图像, 将图像进行渲染, 从不同角度得到如图12和图13所示的渲染图像。通过观察可以发现, 在清晰获得人脸三维信息的同时只有一些突变的区域比如人耳部分产生的错误。

为了进一步说明经参数选择和优化的优越性, 在标定的基础上得出以下的误差图:标定行程为10 cm, 步进电机携标定板每次移动10 mm, 这样在10 cm的行程中总共采集11组图片, 分别为50, 40, 30, 20, 10, 0, -10, -20, -30, -40, -50 (mm) 。选取10 mm处的图片进行分析, 可得到如图13所示的标准差图, 从图中发现, 实际高度为10 mm的平面, 在测量下其高度分布在9.97～10.005 mm之间, 其误差在0.03 mm以内。实验测得标准差为0.022 5 mm。

4结语

在参数选择和优化过程中, 需注意夹角越大测量精度越高, 但是夹角过大, 会使得视场变小, 视场变小会影响测量范围, 增大误差, 因此需压缩夹角。经过多次试验, 夹角在20°～30°效果最好, 精度最高;利用等效波长的概念, 压缩夹角, 提高投影光栅的空间频率, 这样既可以保证系统的测量精度, 又可以保证系统的视场范围。同时, 进行空间带宽积的合理选择, 以及滤波、相位展开等处理, 进一步提高了FTP的测量范围的同时, 系统的测量精度也得到了保证, 达到两全其美的效果。

Fourier分析 篇4

视频抖动的根本原因是由于承载相机的物理部件因为外部因素而震动,从而使得获取图像的传感器进入非理想的运动状态。相机不规则运动必然导致图像发生旋转、平移和尺度上的变化,影响人们的视觉心理。稳像算法的主要目的就是去除图像序列中不规则抖动,输出平滑、稳定的图像序列。

在空域稳像算法中[1,2,3,4]要保证较高的稳像精度,需要有丰富的纹理特征,且对图像形变、光照变化、局部遮挡和噪声等不敏感。文献[1]采用灰度投影算法,计算量小,但在相机调光和曝光过程中,局部最小均方误差(MSE)判别准则会失效,匹配精度将降低,更无法处理旋转、缩放等图像序列。文献[2]采用特征跟踪的方式稳定视频序列,无法处理跟踪目标被遮挡的情况。文献[3,4]提出采用角点特征稳像,对光照变化、缩放、旋转、局部遮挡均有较强的鲁棒性,但是对噪声敏感,需要进行多次高斯平滑和人工阈值设定,计算量大。

1994年Chen[5]等人提出了一种基于Fourier-Mellin变换的图像配准方法,在无需知道相机运动、焦距等物理参数情况下,通过计算两幅图像功率谱的反Fourier变换[5,6]所对应的峰值位置求取两幅图像间的相对平移,通过对幅度谱进行对数-极坐标变换,求取旋转角度和缩放因子。实验验证了该算法不仅精度较高,而且对光照变化、局部遮挡和噪声不敏感。2004年J R Martinez[7]等人将Fourier-Mellin变换用于视频稳像,并利用Cooley-Tukey的FFT算法的对称特性,减小计算量。文献[8]提出了一种单象限Fourier-Mellin域相位相关的近似运动估计算法(SQFMT),在利用了FFT变换的对称特性的基础上,采用下采样方法,在一定程度上提高了计算性能,精度上勉强得到保证,但距工程应用仍然较远。

为克服以上各种算法上的不足,本文提出一种改进的单相限实数Fourier-Mellin(SRFMT)算法,首先基于复数傅里叶变换(CFT)频谱共轭对称性,改进了传统实数傅里叶变换(RFT),并提出根据图像尺寸预先计算蝶形运算过程所要重复涉及的wi=e2πi/N(i=1...N)因子,进一步缩减计算量。同时分析分块运算对稳像算法鲁棒性和精确度的影响,并提出了优化的分块算法,利用统计直方图运动滤波处理块运动向量,使计算结果在保证很高精度的同时,却有了较强的鲁棒性。实验部分给出本文算法结果,证实了在频域内实现实时稳像的可能性。

1 Fourier-Mellin相位相关稳像算法与不足

稳定视频序列的核心是估计帧间旋转、平移和缩放等全局运动参数,帧间图像运动的基本模型如式(1):

式中:fn(x,y)和fn-1(x,y)为两连续帧,θ、(Tx,Ty)、s为旋转、平移和缩放全局运动参数,式(1)的幅度谱经对数-极坐标变换后,即Fourier-Mellin变换得式(2):

式中:(ρ,θ)为对数极坐标,log s与θ0为所求的尺度因子与旋转角度。

相位相关算法获取需要进行四次FFT正变换和两次反变换,求解过程均采用在笛卡尔坐标和对数-极坐标系下归一化功率谱的脉冲函数,通过计算脉冲峰值的位置来确定运动补偿中所需的全局运动参数。在计算过程中,常利用频谱共轭对称性,截取频谱的上半平面减少计算开支,对整个计算性能有一定提升。当前Fourier-Mellin算法有以下不足:1)Fourier-Mellin运算对象为复数,而图像像素均为实数表示,因此存在巨大计算冗余。2)稳像算法有自身的特点,可以通过局部运动参数来估计全局运动向量。因遮挡,局部运动等因素导致的个别不正确局部运动参数估计,并不完全影响全局参数估计置信度。3)在尺度、旋转和平移校正中,不宜采用相同计算方法对待,导致计算量过大而计算精度和算法鲁棒性难以保证。因此,结合当前Fourier-Mellin相位相关稳像算法的不足与稳像算法本身特点,本文将采用第三节与第四节的方法进行改进,保证亚像素稳像精度,提升实时性和算法鲁棒性。

2 改进的相位校正稳像算法

2.1 实数傅里叶变换

若f(x)为实数,那么其频谱F|(ω)|是关于垂直轴几何对称。同样对于二维实函数F|(ω)|,其频谱F|(ωx,ωy)|是关于原点对称。图1为二维实函数傅里叶变换的一个简单实例。

根据共轭对称性,二维实函数f(x,y)在频域的二维复数序列F(u,v)有如下性质:

Re{F(u,v)}与Im{F(u,v)}分别为F(u,v)的实部和虚部。算法主要性能取决于FFT运算量的大小,如果选取当前帧之前一帧作为参考帧,那么要完成旋转、尺度和平移校正,需分别对参考帧和当前帧做2次傅里叶前向变换和2次傅里叶反向变换以及一次对数极坐标变换,那么完成一帧图像的稳定共4次前向FFT、2次反向FFT和1次对数极坐标变换。本文将实数傅里叶变换应用于稳像算法,并加以改进。传统实数傅里叶变换基本思想是:傅里叶运算对象针对的是复数,而图像是以实数表示,将含有N点实数的图像按奇偶位置顺序分离重新组合成一个N/2的复数序列,然后进行FFT运算,最后利用共轭对称性与线性特性对频谱奇偶分解。在传统实数傅里叶变换上,本算法基于频谱的共轭对称性,在对实数FFT变换后的频域N/2点错格奇偶分解[8,9]时取单相限频谱,完成空间幅度谱到对数极坐标空间谱的映射。同时为减少计算量,将当前帧在完成一次傅里叶前向计算后,划分一定空间存储以便下一帧运动估计需要,当前帧计算的全局运动补偿参数与前一帧校正参数之和为当前帧校正。

除此以外,考虑Cooley-Tukey算法的蝶形运算中wi=e2πi/N因子,由于在FFT计算过程中为常数值,且只取决于长度N,而用于FFT运算的向量长度仅与图像的分辨率有关,因此在FFT运算开始前完成wi=e2πi/N因子的计算,以后每次FFT运算过程中无需再次重复计算该因子,只需以查表方式完成整个傅里叶变换。

2.2 双线性插值旋转校正

在Fourier-Mellin算法中,极径ρ分辨率的高低将直接影响运算时间,因此,在不影响估算精度的前提下,不宜采用过细的采样精度。同理,θ的采样精度通常选0.01°,并根据实际应用场合限定在一定的角度范围。文献[8]提出了一种下采样单相限Fourier-Mellin算法提取ρ,θ,有一定的可取之处,速度上有明显优势,但是精度有限,不免会影响平移校正的精度。当图像旋转校正后,将会出现退化现象,使稳像后图像质量降低。因此,结合文献[8]的优越性,采用以下改进措施:利用双线性插值对经旋转变换和尺度变换后的图像做插值处理,提高图像质量,方法简单,效率较高。

2.3 分块平移校正

稳像平移校正主要利用局部块运动向量来确定全局运动向量,优化的分块方案,可有效增强算法的鲁棒性,提高计算精度与速度。图2(a)为视频序列中被截单帧512 pixels×512 pixels图像,被分割为9块,在计算过程中采用内框计算与外框计算相比,精度上没有任何差异,而内框计算时间比外框少很多,那么适当的选取图像块大小有益于减少计算时间,图像块大小选取采用如下关系式进行约束:

WBlk、HBlk为选取块的宽度与高度,Woffset、Hoffset分别为平移运动偏移量,与物理系统有关,可以通过事先实验测量获取,对于无法测量的手持设备,可以根据实际情况设定一个阈值。图像分割块数根据块的大小而确定,通常采用3×3的分割方式。

对于图像发生局部遮挡,出现如图2(a)、2(b)、2(c)、2(d)四种形式。会导致被遮挡块运动估计结果失效。图像在未被遮挡的情况下,统计直方图集中度高,否则会出现如图3中第83帧统计直方图现象,被遮挡块运动估计参数将远离未被遮挡的块。当被遮挡像素面积一定,块选取愈小,遮挡对运动估计影响越小。同时,被遮挡的图像块数愈多,相应块估计也会受到影响,从而影响全局运动参数估计,最终将影响算法鲁棒性。

3 统计直方图运动滤波

利用统计直方图处理块运动向量,获得用于运动补偿的全局运动向量。图3为提取的实验视频序列中第82～86连续5帧的块运动向量统计直方图,横轴为运动偏移量,纵轴为该偏移量所对应的图像块数。当块运动估计结果在直方图的集中度,满足如下条件认为估计结果有效:

其中:bmvRange[n]=bmvmax[n]-bmvmin[n],m取过去的1,…,5帧,n为当前帧;

其中:bmvThrmax=round(.06×2×BMVmax),bmvThrmin=round(.035×2×BMVmax),bmvmax[n]为第n帧块运动向量中的最大值,BMVmax为先前帧中的最大值;bmvmin[n]为第n帧块运动向量中的最小值。利用bmvDynamicThr[n]来表示块运动向量的动态阈值。采用式(7)通过动态阈值求取集中度。

在满足集中度条件下,如果某一运动估计值大于N/2(N为图像块数),认为该运动估计值为全局运动估计值。否则通过加窗滑动处理原始统计直方图,窗口大小等于bmvDynamicThr[n]。将窗内非零估计值与窗口中心值相加如此循环,如果最大值满足是其它所有值之和的2倍,就认定该运动估计值为全局补偿运动参数。此方法主要用于去除伪运动估计和窗口内有相等块数的估计值对全局估计的影响。图3经过直方图滤波以后,结果如图4所示。

4 实验与分析

实验所采用硬件平台位Intel P4 3.0 GHz,内存2 G,Windows操作系统下VC2005编译环境。测试灰度视频序列共132帧,分辨率为256 pixels×320 pixels。采用改进后的单相限实数傅里叶变换与统计直方图滤波,测试算法的时间性能和鲁棒性。

4.1 改进的Fourier-Mellin算法时间性能测试实验

将分辨率为256 pixels×320 pixels的图像序列插值处理生成512 pixels×512 pixels的图像序列后,在稳像算法中,取当前帧的前一帧作为参考帧。在相位校正过程中,取下采样因子为2,预先计算wi=e2πi/N,此时N=512×512/4,并放入预先开辟的存储空间,待每次FFT运算时以查表方式访问;利用改进后的单相限实数傅里叶变换,如图2(a)经极坐标变换后为图5,图6是图5经实数Fourier-Mellin变换后的频谱图,在计算过程中仅需取单相限频谱,通过空间频谱映射,完成FFT运算,并对角度补偿后的图像进行双线性插值以改变图像的PSNR性能。在平移校正过程中,通过预先对物理系统的震动测试,本图像序列是通过算法方式获取,可以认为是预先通过对物理系统测试后获得,其偏移量在-40～40 pixels内,因此在将图像分为9个大小相等的块后,无需下采样,以各块中心为基准,选取每块大小为128 pixels×128 pixels,此时对块进行FFT运算的过程与相位校正中使用方法一样,仅wi=e2πi/N中N=128×128/2。并利用统计直方图对块运动向量滤波,获得用于运动补偿的全局运动向量。测试结果如表1,表2。

采用改进后的单相限实数傅里叶变换,旋转补偿的时间相对于CFT将近提升了70%,平移补偿的计算时间也提升了70%,稳定一帧图片的时间约为39.6 ms,满足实时性要求。相对于文献[9],本文算法性能提升了将近42%。同时可以看出,平移补偿的计算与块的大小选择有着密切关系,较小的块分割必然减少计算时间,但分割方式受物理部件的约束。

4.2 改进的Fourier-Mellin算法鲁棒性测试实验

利用人工方式对测试的132帧图像序列的40帧以后的10图像加入局部面积为80 pixels×80 pixels随机遮挡,用来测mulated0试扰算的法能补力偿。抖以动平视移频估缩计放中、的旋水转平和分平量移估过计程为中例抵,抗采外用部文干献ntalaccumotions-10[8]的算法,下采样因子为4,同时未采用直方图滤波算法,orizo-20本文算法在平移过程中不使用下采样,并对平移估计采用YH-30统计直方图方法,依据集中度对块运动向量进行滤波处理,两者相比,估计结果如图7。

带黑色方块的曲线为真实偏移量,下圆点黑色曲线为本文算法估计的水平平移运动量,而黑色的带上三角符号的曲线为采用文献[8]估计的平移运动量。可以看出本文改进后的算法与真实偏移量有很高的一致性,相对于改进前算法,抵抗外部随机干扰的能力进一步增强。本文以相对真实偏移量的误差来说明抵抗外部干扰的能力,本文算法相对于真实偏移量的均方误差为0.008;改进前的算法相对于真实偏移量的均方误差为0.663。从给估计结果带来的误差方面来说,本文算法的鲁棒性得到很大增强。

图8为任意抽取原始图像序列,上面为被人工遮挡随机图像抽取序列,下面为稳像后的效果。

5 结论

在保证亚像素稳像精度的前提下,提出了一种改进的单相限实数Fourier-Mellin变换(SRFMT)相位校正稳像算法,稳像的实时性得到了较大提高,较改进前提高了约42%,单帧稳像时间在40 ms以内,满足实时性要求;同时采用优化的分块方案,运用统计直方图算法,利用集中度原则,判定块运动向量,去除伪运动估计,解决了因图像局部遮挡而引起的错误估计问题,使算法的鲁棒性得到进一步提高,方法简单有效。本文算法性能的提升,一方面通过改进Fourier-Mellin算法本身,采用实数Fourier-Mellin算法取代传统的复数Fourier-Mellin算法,利用共轭对称性,采用单相限频谱实现,另一方面,通过改进计算方法,避免重复的单帧傅里叶运算以及提前蝶形计算过程中涉及的wi=e2πi/N因子,而计算本身对硬件没有特殊要求,便于实际运用实现。

参考文献

[1]孙辉,张葆,刘晶红,等.航空光电成像电子稳像技术[J].光学精密工程,2007,15(8):1280-1286.SUN Hui,ZHANG Bao,LIU Jing-hong,et al.Electronic image stabilization for aerial E2O imaging system[J].Optics and Precision Engineering,2007,15(8):1280-1286.

[2]王斌,赵跃进,尹德森.基于电子稳像的特征跟踪算法[J].红外与激光工程,2008,37(4):607-610.WANG Bin,ZHAO Yao-jin,YIN De-sen.Digital image tracking algorithm based on electronic image stabilization[J].Infrared and Laser Engineering,2008,37(4):607-610.

[3]CAI J,WALKER R.Robust video stabilization algorithm using feature point selection and delta optical flow[J].IET Computer Vision(S1751-9632),2009,3(4):176-188.

[4]BATTIATO S,GALLO G,PUGLISI G,et al.SIFT Features Tracking for Video Stabilization[C]//14th International Conference on Image Analysis and Processing,Modena,Sept10-14,2007:825-830.

[5]Chen Q,Defrise M,Deconinck F.Symmetric phase only matched filtering of Fourier Mellin transform IBS for image registration and recognition[J].IEEE Translated Pattern Analysis and Machine Intelligence(S0612-8828),1994,16(12):l156-1168.

[6]高莹莹,杨建峰,马晓龙.基于Fourier–Mellin算法的干涉图像配准[J].光学精密工程,2007,15(9):1415-1420.GAO Ying-ying,YANG Jian-feng,MA Xiao-long.Interference image registration based on Fourier–Mellin algorithm[J].Optics and Precision Engineering,2007,15(9):1415-1420.

[7]Martinez de-Dios J R,Ollero A.A real-time image stabilization system based on Fourier-Mellin transform[J].Lecture Notes in Computer Science(S0302-9743),2004,32(11):3762383.

[8]王志民,徐晓刚.基于单象限Fourier-Mellin域相位相关的快速稳像算法[J].中国图象图形学报,2009,14(9):1791-1796.WANG Zhi-min,XU Xiao-gang.Fast Electronic Image Stabilization Based on Phase Correlation in Single2quadran t Fourier-Mellin Domain[J].Journal of Image and Graphics,2009,14(9):1791-1796.

[9]李仕,张葆,孙辉.航空光电成像模糊的实时恢复[J].光学精密工程,2007,15(8):1287-1292.LI Shi,ZHANG Bao,SUN Hui.Real-time restoration for motion blurred in aerial E-O image system[J].Optics and Precision Engineering,2007,15(9):1287-1292.

Fourier分析 篇5

Fourier变换 (Fourier Transform, FT) 作为最主要的信号分析工具主要用于处理频率不随时间变化的平稳信号, 在时频平面时间轴与频率轴相互垂直, 即Fourier变换是从时间域旋转π/2到频率域。但是Fourier变换通常无法表述信号的时频域性质, 不能表示某种频率分量发生在哪个时间, 而这种性质恰恰是非平稳信号最关键的性质, 这对非平稳信号十分重要。分数阶Fourier变换 (Fractional Fourier Transform, FrFT) 是Fourier变换的一种推广形式, 揭示了信号从时间域到频率域变化过程中所呈现的特征, 即从时间域和频率域同时表示信号旋转π/2分数倍时的特征, 从而克服了传统Fourier变换不能反映非平稳信号的统计量随时间变化的缺陷[1,2,3,4]。

分数阶Fourier变换是一种线性变换, 与经典的Fourier变换有着天然的联系, 又提供了Fourier变换所不具备的某些特点, 而且它与小波变换、Wigner-Ville分布 (WVD) 都有密切的关系[3], 因此近年来分数阶Fourier变换受到了研究者的广泛关注, 相应的离散变换算法也相继提出[4]。此外, 有关分数阶Fourier变换应用的研究也越来越多, 如[5,6,7,8,9,10]等文献就介绍了分数阶Fourier变换在信号处理和通信系统方面的应用。

线性调频 (Chirp) 信号是一种典型的非平稳信号, 它的瞬时频率随时间呈线性变化, 常见于雷达、声纳和移动通信等系统中。对Chirp信号的研究, 是非平稳随机信号处理理论及方法的基础。研究证明, 分数阶Fourier变换就是一种很适合处理Chirp信号的变换。近年来, 将分数阶Fourier变换用于线性调频信号 (包括时不变、时变幅度线性调频信号) 的检测和参数估计引起了越来越多的关注[11,12]。本文介绍了FrFT的定义、性质和简单应用, 分析了高斯脉冲Chirp信号的FrFT以及最优化分析角度的选取问题。

1 分数阶Fourier变换的基础

分析和处理平稳信号最常用和最主要的方法是Fourier变换 (FT) 。Fourier变换建立了信号整个时域与整个频域的对应关系, 其中:

$\begin{array}{l}X(\omega )=1/\sqrt{2\text{π}}{\displaystyle {\int }_{-\infty }^{\infty }x}(t)\text{e}{}^{-\text{j}\omega t}\text{d}t(1)\\ x(t)=1/\sqrt{2\text{π}}{\displaystyle {\int }_{-\infty }^{\infty }X}(\omega )\text{e}{}^{\text{j}\omega t}\text{d}t(2)\end{array}$

时域和频域构成了分析一个信号的两种表达方式, Fourier变换在整体时域上将信号分解为不同的频率分量, 但它没有将时域和频域组合成一个域, 不能提供时间和频率的联合信息。谱X (ω) 只是显示任一频率ω包含在信号x (t) 内的总强度, 无法表述信号的时频域性质, 不能表示某种频率分量发生在哪个时间, 而这种性质恰恰是非平稳信号中的最根本和最关键的性质。如果x (t) 是由几个非平稳分量组成的, 那么时间上的任何变化都会改变X (ω) 。此时, 传统的Fourier变换就不能满足信号分析的要求了。

时频分析的基本思想就是设计时间和频率的联合函数, 用它描述信号在不同时间和频率上的能量密度或强度, 从而克服传统Fourier变换不能反映非平稳信号的统计量随时间变化的缺陷。分数阶Fourier变换 (FrFT) 是借用时频面的概念, 以时间和频率分别为横轴和纵轴, 旋转一定的角度进行的一种线性变换。

传统的Fourier变换X (ω) 就是x (t) 旋转π/2, 即x (t) 由时间轴t变到频率轴ω的表示形式。令:

$\alpha =p\text{π}/2(3)$

并且定义线性算子:

$R{}^{\alpha }=R{}^{p\text{π}/2}(4)$

记作Fp。F2相当于t轴连续两次逆时针旋转π/2, 得到x (-t) ;F3相当于t轴连续三次逆时针旋转π/2, 得到指向-ω轴的函数;F4表示t轴连续四次逆时针旋转π/2, 得到原函数。因此线性算子Rα有以下数学性质:

(1) 零旋转:

$R{}^0=x(t)(5)$

(2) 与Fourier变换等价:

$R{}^{\text{π}/2}=F[x(t)](6)$

(3) 旋转相加性:

$R{}^{\alpha }R{}^{\beta }=R{}^{\alpha +\beta }(7)$

(4) 2π旋转 (恒等变换) :

$R{}^{2\text{π}}=x(t)(8)$

如果角度α以π/2的非整数倍进行旋转, 则得到函数x (t) 的广义Fourier变换, 记作:

$\{R{}^{\alpha }x\}(u)=\{F{}^{p}x\}(u)=X{}_p(u)(9)$

这就是x (t) 的分数阶Fourier变换:

Xp (u) ={Fpx} (u) =∫∞-∞x (t) Kp (t, u) dt (10)

式中:变换核Kp (t, u) 定义为:

${\rm K}{}_p(t,u)=\{\begin{array}{l}\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2\text{π}}}\text{e}{}^{\text{j}\frac{t{}^2+u{}^2}{2}\cot \alpha -\text{j}ut\text{c}\text{s}\text{c}\alpha },\alpha \ne n\text{π}\\ \delta (t-u),\alpha =2n\text{π}\\ \delta (t+u),\alpha =(2n+1)\text{π}\end{array}(11)$

式中:n为整数。

从以上的定义可以看出, 分数阶Fourier变换是经典Fourier变换的广义形式, 它包含了信号的时间域和频率域表示。旋转角度为π/2时, 即阶数为1的分数阶Fourier变换就是传统的Fourier变换;不旋转或旋转的角度为2π的整数倍则为信号本身;当旋转角度不在以上两个位置即p为分数时, 它同时从时间域和频率域给出了信号的特征。在讨论分数阶Fourier变换时, 由于角度以2π为模, 所以只需要考虑$0\le \left|p\right|\le 2$的旋转阶数。α角度的分数阶Fourier反变换对应着-α角度的分数阶Fourier变换, 即:

$x(t)={\displaystyle {\int }_{-\infty }^{\infty }X}{}_p(u){\rm K}{}_{-p}(t,u)\text{d}u(12)$

Fourier变换在线性系统分析、光学系统、信息处理系统等方面起着核心作用, 并应用于众多的工程技术领域。作为广义Fourier变换, FrFT具有比Fourier变换更普遍的特性和更广的应用场合, 尤其是普通Fourier变换技术不能解决问题的场合, FrFT更显其优越性。在FrFT研究的启示下, 许多学者推广了分数阶的概念, 得到分数阶卷积和分数阶相关, 还把分数阶的概念应用到Hadamard变换、Hartley变换等, 得到相应的分数阶变换。而且, 分数阶Fourier变换已应用到图象处理的优化图象恢复方面。此外, 将分数阶Fourier变换算子的分数幂推广至复数幂也成为目前研究的一个热点。

在20世纪90年代中期, 分数阶Fourier变换被引入了信号处理领域。在信号处理中, FrFT有很多应用, 其中两个典型的应用是信号滤波和信号分离[5]。实验已经证实分数阶Fourier变换滤波的效果明显优于Fourier变换;特别是对于线性调频Chirp信号, 分数阶Fourier变换能够获得最佳的能量积聚效果, 是分数阶Fourier变换最合适的应用领域之一。

2 高斯脉冲Chirp信号的分数阶Fourier变换

线性调频Chirp信号是一种特殊的非平稳信号, 它的瞬时频率随时间呈线性变化, 广泛地出现在通信、雷达、声纳和地震勘探等系统中。在工程实践中, 高斯调制Chirp信号有着广泛的应用, 而FrFT对于Chirp信号良好的检测和分析效果自然让我们想到能否将此分析用具引入到高斯调制的Chirp信号分析中。本文对此进行了深入研究, 给出了相关解析结果。

高斯信号的标准形式为:

$f(t)=\frac{1}{\sqrt{2\text{π}\delta {}^2}}\text{e}{}^{-\frac{(t-t{}_0){}^2}{2\delta {}^2}}(13)$

这是一个服从正态分布 (t0, δ) 的函数, 它的Fourier变换和分数阶Fourier变换都具有高斯信号的形式[12]。利用已有的这些性质, 可以研究高斯脉冲Chirp信号的分数阶Fourier变换。

令高斯脉冲Chirp信号为:

$x(t)=\text{e}{}^{-\frac{(t-t{}_0){}^2}{2\delta {}^2}}\text{e}{}^{\text{j}(at{}^2+bt+c)}(14)$

将高斯脉冲Chirp信号写成幅度和相位的函数, 得到:

$x(t)=\left|x(t)\right|\text{e}{}^{\text{j}\theta (\omega )}(15)$

其中:

$\begin{array}{l}\left|x(t)\right|=\text{a}\text{b}\text{s}(\text{e}{}^{-\frac{(t-t{}_0){}^2}{2\delta {}^2}}\text{e}{}^{\text{j}(at{}^2+bt+c)})=\text{e}{}^{-\frac{(t-t{}_0){}^2}{2\delta {}^2}}\\ \theta (\omega )=at{}^2+bt+c\end{array}(16)$

式 (14) 所表示的高斯脉冲Chirp信号的分数阶Fourier变换为:

$\begin{array}{l}X{}_p(u)=\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2\text{π}}}\text{e}{}^{\text{j}\frac{u{}^2}{2}\cot \alpha }{\displaystyle \int \text{e}}{}^{-\frac{(t-t{}_0){}^2}{2\delta {}^2}}\text{e}{}^{\text{j}(at{}^2+bt+c)}\text{e}{}^{\text{j}\frac{t{}^2}{2}\cot \alpha }\text{e}{}^{-\text{j}ut\text{c}\text{s}\text{c}\alpha }\text{d}t\\ =\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2\text{π}}}\text{e}{}^{\text{j}c}\text{e}{}^{\text{j}\frac{u{}^2}{2}\cot \alpha }\text{e}{}^{-\frac{t{}_0^2}{2\delta {}^2}}{\displaystyle \int \text{e}}{}^{-[(\frac{1}{2\delta {}^2}-\text{j}a-\text{j}\frac{\cot \alpha }{2})t{}^2-(\frac{t{}_0}{\delta {}^2}+\text{j}b-\text{j}u\text{c}\text{s}\text{c}\alpha )t]}\text{d}t\\ ={\rm P}{\displaystyle \int \text{e}}{}^{-({\rm M}{}^{2}t{}^2-2{\rm M}{\rm N}t)}\text{d}t(17)\end{array}$

其中:

$\begin{array}{l}{\rm P}=\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2\text{π}}}\text{e}{}^{\text{j}c}\text{e}{}^{\text{j}\frac{u{}^2}{2}\cot \alpha }\text{e}{}^{-\frac{t{}_0^2}{2\delta {}^2}}\\ {\rm M}=\sqrt{\frac{1}{2\delta {}^2}-\text{j}a-\text{j}\frac{\cot \alpha }{2}}\\ {\rm N}=\frac{\frac{t{}_0}{\delta {}^2}+\text{j}b-\text{j}u\text{c}\text{s}\text{c}\alpha }{\sqrt{\frac{1}{2\delta {}^2}-\text{j}a-\text{j}\frac{\cot \alpha }{2}}}\end{array}(18)$

所以式 (17) 化为:

$\begin{array}{l}X{}_p(u)={\rm P}{\displaystyle \int \text{e}}{}^{-({\rm M}{}^{2}t{}^2-2{\rm M}{\rm N}t+{\rm N}{}^2-{\rm N}{}^2)}\text{d}t\\ ={\rm P}\text{e}{}^{{\rm N}{}^2}{\displaystyle \int \text{e}}{}^{-({\rm M}t-{\rm N}){}^2}\text{d}t(19)\end{array}$

令Mt-N=Q, 则$t=\frac{Q+{\rm N}}{{\rm M}},\text{d}t=\frac{\text{d}Q}{{\rm M}}$, 因此得到:

$X{}_p(u)=\frac{{\rm P}\text{e}{}^{{\rm N}{}^2}}{{\rm M}}{\displaystyle \int \text{e}}{}^{-Q{}^2}\text{d}Q(20)$

然后取积分限为 (-∞, ∞) , 由于${\displaystyle {\int }_{-\infty }^{\infty }\text{e}}{}^{-Q{}^2}\text{d}Q=\sqrt{\text{π}}$, 所以可得:

$X{}_p(u)=\frac{\sqrt{\text{π}}{\rm P}\text{e}{}^{{\rm N}{}^2}}{{\rm M}}(21)$

再将P, M, N的值代入式 (21) 得:

$\begin{array}{l}X{}_p(u)=\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2}}\text{e}{}^{\text{j}c}\text{e}{}^{\text{j}\frac{u{}^2}{2}\cot \alpha }\text{e}{}^{-\frac{t{}_0^2}{2\delta {}^2}}\cdot \\ \frac{1}{\sqrt{\frac{1}{2\delta {}^2}-\text{j}a-\text{j}\frac{\cot \alpha }{2}}}\text{e}{}^{\frac{\left(\frac{t{}_0}{\delta {}^2}+\text{j}b-\text{j}u\text{c}\text{s}\text{c}\alpha \right){}^2}{\frac{2}{\delta {}^2}-\text{j}4a-\text{j}2\cot \alpha }}(22)\end{array}$

将$\frac{(t{}_0/\delta {}^2+\text{j}b-\text{j}u\text{c}\text{s}\text{c}\alpha ){}^2}{2/\delta {}^2-\text{j}4a-\text{j}2\cot \alpha }$的分子分母同时乘以2/δ2+j4a+j2cot α, 则式 (22) 化为:

$X{}_p(u)=\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2}}\text{e}{}^{\text{j}c}\text{e}{}^{\text{j}\frac{u{}^2}{2}\cot \alpha }\text{e}{}^{-\frac{t{}_0^2}{2\delta {}^2}}\frac{1}{\sqrt{\frac{1}{2\delta {}^2}-\text{j}a-\text{j}\frac{\cot \alpha }{2}}}\cdot$

$\text{e}{}^{\frac{\left(\frac{t{}_0}{\delta {}^2}+\text{j}b-\text{j}u\text{c}\text{s}\text{c}\alpha \right){}^2\left(\frac{2}{\delta {}^2}+\text{j}4a+\text{j}2\cot \alpha \right)}{\frac{4}{\delta {}^4}+(4a+2\cot \alpha ){}^2}}(23)$

若t0=0, 可以得到:

$X{}_p(u)=\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2}}\frac{1}{\sqrt{\frac{1}{2\delta {}^2}-\text{j}a-\text{j}\frac{\cot \alpha }{2}}}\cdot$

$\text{e}{}^{\text{j}c}\text{e}{}^{\text{j}\frac{u{}^2}{2}\cot \alpha }\text{e}{}^{\frac{-\text{j}(4a+2\cot \alpha )(u\text{c}\text{s}\text{c}\alpha -b){}^2}{\frac{4}{\delta {}^4}+(4a+2\cot \alpha ){}^2}}\text{e}{}^{\frac{-\frac{2}{\delta {}^2}(u\text{c}\text{s}\text{c}\alpha -b){}^2}{\frac{4}{\delta {}^4}+(4a+2\cot \alpha ){}^2}}(24)$

将Chirp信号的分数阶Fourier变换写成幅度和相位的函数, 得到:

$X{}_p(u)=\left|X{}_p(u)\right|\text{e}{}^{\text{j}\rho (\omega )}(25)$

因为$\left|\text{e}{}^{\text{j}A}\right|=1$, 在计算Xp (u) 的幅值时可以忽略所有模为1的部分, 这样就能得到:

$\left|X{}_p(u)\right|=$

$\text{a}\text{b}\text{s}\left(\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2}}\frac{1}{\sqrt{\frac{1}{2\delta {}^2}-\text{j}a-\text{j}\frac{\cot \alpha }{2}}}\text{e}{}^{\frac{-(u-b\text{s}\text{i}\text{n}\alpha ){}^2}{2\sin {}^2\alpha [\frac{1}{\delta {}^2}+\delta {}^2(2a+\cot \alpha ){}^2]}}\right)(26)$

令:$\sigma {}^2=\sin {}^2\alpha [\frac{1}{\delta {}^2}+\delta {}^2(2a+\cot \alpha ){}^2](27)$

得:

$\begin{array}{l}\left|X{}_p(u)\right|=\text{a}\text{b}\text{s}[\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2}}\\ \frac{1}{\sqrt{\frac{1}{2\delta {}^2}-\text{j}a-\text{j}\frac{\cot \alpha }{2}}}\text{e}{}^{\frac{-(u-b\text{s}\text{i}\text{n}\alpha ){}^2}{2\sigma {}^2}}](28)\end{array}$

由式 (24) 和式 (25) 得:

$\text{e}{}^{\text{j}\rho (\omega )}=\text{e}{}^{\text{j}\left[c+\frac{u{}^2}{2}\cot \alpha -\frac{(4a+2\cot \alpha )(u\text{c}\text{s}\text{c}\alpha -b){}^2}{\frac{4}{\delta {}^4}+(4a+2\cot \alpha ){}^2}\right]}(29)$

因此Xp (u) 的相位:

$\rho (\omega )=c+\frac{u{}^2}{2}\cot \alpha -\frac{(4a+2\cot \alpha )(u\text{c}\text{s}\text{c}\alpha -b){}^2}{\frac{4}{\delta {}^4}+(4a+2\cot \alpha ){}^2}$ (30)

由式 (28) 可得$\left|X{}_p(u)\right|$具有高斯函数的形式, 服从正态分布 (bsin α, σ2) 。当方差σ2最小时, 高斯脉冲Chirp信号的分数阶Fourier变换Xp (u) 的能量最大, 因此这时的α就是最优化角度αopt。从式 (27) 可知, σ2最小时, 有2a+cot αopt=0, 即:

$\alpha {}_{\text{o}\text{p}\text{t}}=-\arctan \frac{1}{2a}(31)$

同样, 由式 (3) 可知, 此时Xp (u) 的最优化阶数popt为:

$p{}_{\text{o}\text{p}\text{t}}=-\frac{2}{\text{π}}\arctan \frac{1}{2a}(32)$

由式 (31) 和式 (32) 对比文献[12]可得, Chirp信号和高斯脉冲调制Chirp信号的αopt, popt是一样的, 这说明加上高斯窗的Chirp信号 (即高斯脉冲Chirp信号) 的最优化参数并没有改变。这就为高斯Chirp信号的FrFT最优化分析奠定了坚实的理论基础, 可以充分利用FrFT对于LFM信号最优的能量积聚性质, 将其应用于高斯Chirp信号的分析中, 在低信噪比的环境中这一点更具有实际意义。

3 结 语

分数阶Fourier变换是Fourier变换的广义形式。作为一种新的时频分析工具, 分数阶Fourier变换既与经典的Fourier变换有着天然的联系, 又提供了Fourier变换所不具备的某些特点, 而且它与小波变换、Wigner-Ville分布 (WVD) 都有密切的关系[12,13], 是一种很适合处理线性调频Chirp信号等非平稳信号的变换。

在Chirp信号的检测和参数估计方面, 分数阶Fourier变换具有重要的应用价值。尤其是它所具有的线性特性, 可以方便地用于多分量线性调频信号的检测和参数估计而不受交叉项的干扰, 大量被用于恒定幅度的线性调频信号, 并取得了良好的效果。事实上, 时变幅度线性调频信号在工程实际中有着更为广泛的应用背景, 在工程实践中, 会遇到很多经过幅度调制的线性调频信号, 其中的高斯窗调制就是典型的一种。

从分析结果看, 经过高斯窗调制的Chirp信号保持了原本的特性, 在FrFT这个分析工具下能够获得优异的性能。

Fourier分析 篇6

chirp类信号是一种典型的非平稳时变信号,它的瞬时频率随时间呈线性变化,可以看成是频域时变信号的一阶近似。因此,FRFT很适合处理chirp信号的变换,将其用于chirp信号(包括时不变、时变幅度线性调频信号)的检测和参数估计是近年来兴起的热点方向[1,2,3]。针对FRFT数值计算的快速算法中,当chirp信号的最高频率超过采样频率的一半时,即处于欠采样,仍然可以精确地估计出调频率值的问题[4],对FRFT进行研究,提出了采用分数阶Fourier域的采样定理来解决。很多文献对采样定理从不同方面进行了研究[5,6]:文献[6]提出了原始信号与采样信号的FRFT关系,但运算过程复杂,且结算量大,没有完整的说明量纲归一化的步骤,缺乏仿真实例的验证;文献[5]中简化了运算过程,降低了运算复杂度,还采用了仿真实例进行验证,但是没有讨论加入噪声时对欠采样信号还原的影响,以及没有进一步在三维空间中对chirp信号的FRFT进行搜索仿真。

本文尝试利用分数阶Fourier域的采样定理,将chirp信号通过FRFT,取不同调频数值,即具体分析在其对应的分数阶Fourier域欠采样时采样信号实现能量聚焦的情况。

1 信号的FRFT定义

FRFT实质是一种统一的时频变换,同时反映了信号在时频域的信息,与通常的二次型时频分布的区别是它的单一变量表示时频信息,且没有交叉项这个难点,与传统的Fourier变换相比,更适合处理如chirp类非平稳信号。下面从线性积分变换方面得到FRFT的定义

$\begin{array}{l}f{}_p(u)\equiv {\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\rm K}}{}_{{\rm P}}(u,t)f(t)\text{d}t=\\ \{\begin{array}{l}\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2\text{π}}}{\displaystyle {\int }_{-\infty }^{+\infty }\text{e}}\text{x}\text{p}\left(\text{j}\frac{u{}^2+t{}^2}{2}\cot \alpha -\frac{\text{j}ut}{\sin \alpha }\right)f(t)\text{d}t‚\\ \alpha \ne n\text{π}\\ f(t)‚\alpha =2n\text{π}\\ f(-t)‚\alpha =(2n\pm 1)\text{π}\end{array}(1)\end{array}$

式中:fp(u)是原函数f(t)在变换域的表示;而KP(u,t)≡Aαexp[jπ(u2cotα-2utcscα+t2cotα)],称为FRFT的核函数;$A{}_{\alpha }\equiv \sqrt{1-\text{j}\text{c}\text{o}\text{t}\alpha }$;α≡pπ/2,p≠2n,n为整数。

2 chirp信号采样

下面进行分数阶Fourier域带通信号的采样定理推导[5]。

设chirp信号(看作原始信号)为g(t),g(t)的分数阶Fourier变换域上的带限信号满足:G(u)的带宽为Φ=Φh-Φ1,则当u>Φh或u<Φ1,0<Φ1<Φh时,G(u)=0。其中G(u)为g(t)经α阶分数阶Fourier变换域的谱信号。将g(t)与采样冲激串δ(t-nTs)(采样周期是Ts)相乘,得到理想采样信号gs(t)为

$g{}_s(t)=g(t)\times {\displaystyle \sum _{-\infty }^{+\infty }\delta }(t-n{\rm T}{}_s)(2)$

再利用分数阶Fourier卷积定理得到α阶FRFT后的采样谱信号Gs(u)与g(t)经α阶分数阶Fourier变换域的谱信号G(u)的关系

$\begin{array}{l}G{}_s(u)=\frac{1}{{\rm T}{}_s}\times \text{e}{}^{\text{j}\frac{\cot \alpha }{2}u{}^2}\times \\ \left[G(u)\text{e}{}^{-\text{j}\frac{\cot \alpha }{2}u{}^2}\times {\displaystyle \sum _{-\infty }^{+\infty }\delta }\left(u-n\frac{2\text{π}\sin \alpha }{{\rm T}{}_s}\right)\right](3)\end{array}$

通过图1可以看出在分数阶Fourier域的chirp信号谱函数G(t)与其采样信号谱函数Gs(u)之间的变换流程。

从式(3)和图1中给出关系:信号g(t)在时域采样, 那么Gs(u)就看作以周期$\frac{2\text{π}\sin \alpha }{{\rm T}{}_s}$重复,且可被线性调频信号$\text{e}{}^{-\text{j}\frac{\cot \alpha }{2}u{}^2}$加权。那么当$\frac{2\text{π}\sin \alpha }{{\rm T}{}_s}\ge 2\Phi {}_{\text{h}}$(即采样频率ωs≥2Φh/sinα,Φh为信号的频域最高频率)时就不会发生混叠,这样便能设计合适的分数阶Fourier变换域理想低通滤波器来取出零次谐波,从而不失真地恢复出原信号g(t)。传递函数Hα如

${\rm H}{}_{\alpha }=\{\begin{array}{l}{\rm T}{}_s,|u|<\Phi {}_{\text{c}\text{u}\text{t}}\\ 0,\Phi {}_k<\Phi {}_{\text{c}\text{u}\text{t}}<\omega {}_s\sin \alpha -\Phi {}_k\end{array}(4)$

可以得到

2ωh>2Φh/sinα (5)

即gh>uh/sinα(gh为信号的频域最高频率,uh为信号的$\frac{2\alpha }{\text{π}}$阶分数阶Fourier域最高频率)。

满足用低于2Φh的频率来进行采样,再根据图2进行推导:通过时移和频移使其时频分布关于原点对称,即

$\begin{array}{l}u{}_{\text{h}}=u{}_{1,\text{h}}=g{}_{\text{h}}\frac{\cos (\beta -\alpha {}_1)}{\sin \beta },\\ \beta -\alpha {}_1<\text{π}/2\text{且}0<\alpha {}_1,\beta <\text{π}(6)\end{array}$

$\begin{array}{l}u{}_{\text{h}}=-u{}_{2,\text{h}}=g{}_{\text{h}}\frac{\cos (\beta -\alpha {}_2)}{\sin \beta },\\ \beta -\alpha {}_1<\text{π}/2\text{且}0<\alpha {}_2,\beta <\text{π}(7)\end{array}$

式(6)和(7)中:α1,α2是分数阶Fourier阶数;β是chirp信号的Wigner分布与时间轴的逆时针夹角。

联合(6)、(7)式得到

$u{}_{\text{h}}=g{}_{\text{h}}\frac{\cos (\beta -\alpha )}{\sin \beta },0<\alpha ,\beta <\text{π}(8)$

将式(8)代入式(5)则有$\frac{\cos (\beta -\alpha )}{\sin \beta \sin \alpha }<1$,最后得到cotαcotβ<0, 0<α,β<π。

显然,只要α与β不在同一象限,则式(5)成立,可知chirp信号可用低于2gh的采样频率进行采样而不会造成分数阶Fourier谱信号的混叠,相应的p阶分数Fourier域重构低通滤波的截止频率Φcut的取值区间为[2πuh,2πgssinα-2πuh]。

3 信号仿真

下面列举具体的chirp信号来进行探讨验证,设chirp信号模型为g(t)=ejπwt2,该信号的采样率为fs=240 Hz,持续时间Td=2 s,得到采样信号为

$\text{g}{}_{\text{s}}(\text{t})=\text{e}{}^{j\text{ω}\text{t}{}^2}{\displaystyle \sum _{-\infty }^{+\infty }\text{δ}}(\text{t}-\frac{\text{n}}{2}),-1\le \text{t}\le 1(9)$

令分数阶Fourier域的变换角度为αr(cotαr=-120,0<αr<π),利用采样定理,根据式(8)和$f{}_s=\frac{2f{}_{\text{h}}}{\sin \beta {}_r}$,可得到

$\begin{array}{l}f{}_s=\frac{2u{}_{\text{h}}}{\sin \alpha {}_r}=2f{}_{\text{h}}\frac{\cos (\alpha {}_r-\beta )}{\sin \beta \sin \alpha {}_r}=2f{}_{\text{h}}(1+\cot \alpha {}_r\cot \beta )=\\ \lambda {}_c{\rm T}{}_{\text{d}}\left(1+\frac{\cot \alpha {}_r}{2\text{π}\lambda {}_c}\right)(10)\end{array}$

将数值cotαr=-120代入式(10),求出λc=140 Hz/s,即为采样信号gs(t)的αr角度分数阶Fourier谱不混叠的临界调频率。为了验证不同调频率λ(分小于λc、等于λc、大于λc时3种情况)对采样的影响,下面对该chirp信号g(t)在用离散尺度法实现归一化的前提下,进行分数阶Fourier变换域上的采样。从图3a,3b,3c可以看出,在λ=100 Hz/s,120 Hz/s,140 Hz/s时,采样信号gs(t)在其分数阶Fourier域的谱信号Gs(t)没有发生混叠,能很好地分离出各次谐波,且在λ=120 Hz/s时,是gs(t)的最优阶次的FRFT,在该变换域上形成了能量聚焦;但是当λ=150 Hz/s,即大于临界调频率λc时,Gs(t)的各次谐波已经发生了较为明显的信号混叠。

在与前面相比,为了更形象地表现出欠采样对Gs(t)的影响,采用对chirp信号g(t)的p阶次(最优阶次)变换信号的能量聚焦的三维空间图进行仿真,取傅里叶变换阶次p=[1,2],在搜索步长为0.01的范围内进行二维搜索的FRFT的仿真,可设加入的高斯白噪声后信噪比等于0 dB,为了验证加入噪声后,欠采样信号的还原情况,设持续时间、采样频率、变换阶次不变的情况下,通过分别取调频率ω=100 Hz/s,120 Hz/s,140 Hz/s,150 Hz/s这几个分别能代表过采样、二倍采样、临界采样值以及欠采样时的典型值,对比图4中各个图像可以看出:

1) 在图4a,4b,4c中看出,当调频率均不大于λc时,Gs(t)仍能在u域的能量聚焦。但是,在图4d中,当λ=150 Hz/s时,超过了最大值λc=140 Hz/s,在参数平面上,混叠信号出现了,还原信号已经失真,不能满足在欠采样时恢复原始信号,这个现象和前面推导的结果是吻合的。

2) 高斯噪声对信号的误差影响相比与在欠采样取近似阶次p时对信号还原的影响还小,而且当λ>120 Hz/s时,即信号处于欠采样状态时,在归一化后,相对误差还是很小,不影响还原原始信号。

4 结语

chirp信号在分数阶Fourier域欠采样情况下,通过仿真结果可以看到采样信号在最优分数阶域变换域的能量聚集很显著,参数平面上,混叠信号一直没有出现,直到采样频率与信号频率相等时,才出现几个信号大致相等的幅值峰值(信号混叠),造成还原信号失真。因此,可以验证下面的推导是可行的:当满足cotαcotβ<0, 0<α,β<π时,就能得到$u{}_{\text{h}}=g{}_{\text{h}}\frac{\cos (\beta -\alpha )}{\sin \beta }$。这样chirp信号能在欠采样的情况下作低通滤波,并能很好地恢复原信号。

摘要：针对chirp信号在分数阶Fourier域上被低于二倍信号频率采样(欠采样)后,仍能恢复原始信号的问题,依据分数阶Fourier域上带限信号的采样定理,分析了分数阶Fourier域上欠采样时应满足的条件,并对其进行了仿真验证。结果表明:对chirp信号进行欠采样后能够得到采样信号的能量聚焦。取不同的调频率,欠采样的采样信号仍能较好地估计出调频率,并能恢复原始信号。

关键词：分数阶Fourier域,分数阶Fourier变换,带限信号,调频率,能量聚焦

参考文献

[1]齐林,陶然,周思永.基于分数阶Fourier变换的多分量LFM信号的检测和参数估计[J].中国科学:E辑,2003,33(8):749-759.

[2]SUN Hongbo,LIU Guosui.Application of the fractional Fourier transformto moving target detection in airborne SAR[J].IEEE Transaction on Aer-ospace and Electronic Systems,2002,38(4):142-146.

[3]董永强,陶然,周思永.基于分数阶傅立叶变换的SAR运动目标检测与成像[J].兵工学报.1999,20(2):132-136.

[4]赵兴浩,邓兵,陶然.分数阶傅里叶变换数值计算中的量纲归一化[J].北京理工大学学报,2005,25(4):360-364.

[5]陶然,邓兵,杨曦.分数阶Fourier域采样及分辨率分析[J].自然科学进展,2007,17(5):655-661.

Fourier分析 篇7

Chirp信号即线性调频 (linear frequency moderate, LFM) 信号, 其瞬时频率随着时间线性变化, 具有抗截获破译能力、抗快速移动造成的衰落和多径分辨能力强等优点, 已广泛应用于无线通信、雷达和声纳等众多领域。Chirp扩频调制[1,2,3]是IEEE802.15.4a定义的两种物理层标准之一, 目前应用较多的调制方式是二进制正交键控 (Binary Orthogonal Keying, BOK) 。BOK调制方式单一, 数据率低, 且不能实现多址。在1962年Winkler提出用Chirp信号进行抗干扰扩频通信时, 即提出采用不同的调频率来携带信息, 即高阶的多进制Chirp-rate调制 (MCr SK) , 可以大大提高信息率。目前MCr SK的解调算法主要是分数阶Fourier变换 (Fr FT) [4,5], 获得良好的性能。但是Fr FT解调算法在某一调频率信号达到最佳匹配滤波时, 其他调频率的信号会成为频带内的噪声, 而没有完全滤除, 即Fr FT解调算法对于不同调频率的chirp信号不具有正交性。本文将双正交Fourier变换 (Biorthogonal Fourier Transform, BFT) [6,7]引入多MCr SK系统, 提出基于BFT的MCr SK解调算法, 解决了多进制调制的正交性问题, 实现了不同的调频率携带不同的信息的数字化通信。

1 BOK调制和BFT变换

在时域Chirp-rate信号可以表示为:

式中为载频, r为调频率。当r为正时, 称为p-Chirp信号;当r为负时, 称为down-Chirp信号。BOK就是利用upChirp信号和down-Chirp信号的准正交性分别来传输二进制信号, 其自相关函数为:

其归一化的自相关和互相关函数如图1所示。

理想的BOK调制信号的自相关函数应该为sinc函数, 而互相关函数应该为零。但是如图1所示, 因为up-Chirp信号和down-Chirp信号不是完全的正交, 匹配输出出现多峰现象, 其自相关函数近似sinc函数, 而互相关函数则是有一定幅值的震荡波形, 这些性质为解调带来了不良影响, 限制了BOK调制的性能提高。

BFT解决了不同调频率Chirp信号间的正交问题, 可以通过不同的调频率的信号同时传递信息而不会相互混淆。BFT把信号采用双正交Chirp基函数展开, 得到信号的调频率密度谱, 它是一种类似于频谱的密度分布。通过BFT, 就像单频信号的频谱是冲激函数一样, 单调频率Chirp信号在斜率谱就会表现为相应的冲激函数, 而且不同调频率的Chirp信号彼此正交。BFT定义为:

式中为时域信号, 为调频率密度谱。BFT可以把Chirp信号压缩成冲激函数, 其在调频率谱的相应位置出现单一的峰值, 而多Chirp函数则是单个Chirp函数的组合, 不存在交叉项, 也不用搜索, 因此适合多调频率的Chirp组合信号分析。

实际中得到的都是信号的离散值, 因此有离散BFT, 通常记作:

式中k和l分别为时间和调频率的离散量, fs为采样频率, N为总的采样点数。

图2分析包含3分量的Chirp信号的BFT的性能, 同时给出了其时频Wigner-Ville分布 (WVD) 如图2 (a) , 以及其RadonWigner变换 (RWT) 如图2 (b) 、分数阶Fourier变换 (FRFT) 如图2 (c) 和BFT如图3 (d) 。首先WVD属于Cohen类分布, 其严重不足是存在交叉项, 会产生虚假信号, 在Chirp信号分量较多时, 其WVD非常混乱。其次Radon变换需要对角度进行搜索, 即不同调频率对应不同角度, 这样某一个角度下, 只有一个信号分量能够得到最佳的估计, 其它信号则不能有效估计。FRFT虽然没有交叉项, 但是同样需要搜索, 一般是先检测最强分量的信号, 然后对其进行遮蔽处理, 再进行第二强分量信号的估计, 一直重复这个过程, 最后得到所有信号的信息, 无疑这个过程是复杂且对噪声敏感的。当Chirp信号调频率变大时, FRFT变得不敏感, 估计精度下降。而BFT通过变换直接精确给出三个分量的调频率和其幅值, 不用搜索, 也没有交叉项, 简单方便。

2 MCr SK调制

MCr SK调制把不同信息调制在Chirp信号的调频率上, 解调时可以实现两个信号的完全分离。MCr SK可表示为:

式中ω为载频, T为码元长度, g为码元窗函数, θ为相位, 其中由r携带所传递的信息。调制时把要传输的数字信号先经过串并转换得到一组多进制的数值, 然后由这组值得到相应斜率的锯齿波电压函数, 再由锯齿波控制压控振荡器, 就得到MCr SK信号, 最后调制到载波上。信号的时频图如图3所示。

根据离散BFT的性质, 其调频率的分辨率为:

式中T为码元长度, 可见其的分辨率只与码长度有关。设每个码元的采样点数为M, 则

在基带调制时, 设每个码元的Chirp信号都从零频开始, 则最大的频率为:

根据采样定理要求, 有:

带入式 (6) 、式 (7) 和式 (8) 至式 (9) 可得:

上式表明了MCr SK的容量, 即在满足最小分辨率和最大采样频率的基础上, MCr SK的最大容量为每码元采样点数的1/4。

3 性能仿真

在BFT解调算法的框架下, MCr SK各个调制分量实现了正交分解, 因此基于BFT的MCr SK的解调属于正交信号相干检测, 图4给出二进制的MCr SK信号的误比特率与信噪比的关系曲线。图4在相同情况下, 二进制的MCr SK的性能和BFSK相近。

4 结语

由MCr SK的调制方式可以看出, MCr SK可以与MFSK、MPSK和MASK三种调制一样, 作为基本的调制解调体制。MCr SK可以充分发挥Chirp信号的一些特点, 具有频带较宽、传输容量大、恒包络和旁瓣衰减迅速等特点, 可广泛用于抗干扰和扩频通信等方面。由于对MCr SK和BFT的研究较少, 对于适用于工程应用的MCr SK调制解调特性, 还要进一步研究。

参考文献

[1]Sarni Y, Sadoun R, Belouchrani A.On the application of Chirp modulation in spread spectrum communication systems[J].IEEE International Symposium on Signal Processing and its Application (ISSPA) .2001 (2) :501-504.

[2]贺鹏飞, 吕英华, 张洪欣, 等.基于Chirp-BOK调制的超宽带无线通信系统研究[J].南京邮电大学学报:自然科学版, 2006, 26 (2) :21-25.

[3]王晓炜, 李昕, 费敏锐.Chirp扩频通信系统抗噪声性能研究[J].通信技术, 2009, 42 (3) :36-37, 50.

[4]邓兵, 陶然, 平殿发.基于分数阶Fourier变换的chirp-rate调制解调方法研究[J].电子学报, 2008, 36 (6) :1078-1083.

[5]肇启明, 张欣宇, 张乃通.多进制chirp-rate键控调制及分数傅里叶变换法解调[J].通信学报, 2010, 31 (6) :1-6.

[6]王本庆, 李兴国.LFM信号调频率的双正交Fourier变换分析算法[J].电子与信息学报, 2009, 31 (7) :1620-1623.