一种阶梯叠加式傅里叶变换微分方法与流程

本发明涉及波谱和成像信号处理领域，具体涉及一种阶梯叠加式傅里叶变换微分方法。

背景技术：

一个随时间t周期性变化的分析信号s(t)可通过傅里叶变换技术转换成频率ω的波谱s(ω),通用的数学表达式为：

波谱s(ω)经过逆傅里叶变换又可变回原有的分析信号s(t)：

上述傅里叶变换积分表达式乘以了一个复变量指数函数e^-iωt或e^iωt，其中

纯谐振分析信号在傅里叶变换后的三种基本峰形：吸收峰形，发散峰形和幅度峰形都是完全对称的。实际信号通常有衰减效应，傅里叶变换产生的基本峰形演变成洛伦兹峰形,或接近于洛伦兹峰形(如高斯峰形以及沃伊特峰形)，它们跟三种基本峰形一样仍然保持对称性。现有技术中实施计算过程比较复杂耗时，类似于波谱学的曲线拟合法,需结合一些先验知识才有可能对有限的信号数据和比较简单的图像进行叠加处理。一般来说提升傅里叶变换的波谱和成像技术的分析性能主要是如何有效地检测出弱小信号和分辨重迭谱峰。

技术实现要素：

有鉴于此，本发明的目的在于提供一种阶梯叠加式傅里叶变换微分方法，不但在分辨率和灵敏度都有很大的改进，而且大大节省了计算时间。

为实现上述目的，本发明采用如下技术方案：

一种阶梯叠加式傅里叶变换微分方法，包括以下步骤：

步骤s1:对分析信号在有限时间域进行采样，并通过离散化处理转换成谱图，通过其微分曲线鉴别谱峰顶点位置；

步骤s2:根据叠加式傅里叶变换原理对谱图中的各峰顶点以阶梯形式叠加；

步骤s3:对阶梯叠加后的谱图做微分，进一步获得阶梯叠加式傅里叶变换微分波谱或图像。

进一步的，所述步骤s1具体为:

设分析信号s(t)在一测量时间段离散化取n个样本，s(0)、s(1)、s(2)、···、s(n-1)；用离散化傅里叶变换得到n个频率谱数据s(0)、s(1)、s(2)、···、s(n-1)，傅里叶变换矩阵表达如下：

式中，傅里叶变换矩阵的w＝exp(-i2π/n)；

对傅里叶变换波谱s(ω)乘上一个对角矩阵diag[ak]：

式中如果有m个重迭在一块的峰，对角线元素ak取值为2n，n＝0，1，2，···，m；

阶梯叠加对角矩阵进一步写成：

通过微分学二阶导数的鉴别，定位出谱图曲线的峰顶点ω1、ω2、ω3和ω4。

进一步的，所述步骤s2具体为：

设谱图中有n个重迭峰，n个重迭峰分别f1、f2、f3、...、fn，首个峰f1定位在零阶梯，根据右叠加函数它的峰强度被乘以2，基线仍然是原有谱峰基线0；第二个峰f2的峰强度被乘以4，上升到第一个阶梯，依次将第三个峰f3乘以6，升到第二个阶梯...第n个峰fn乘以2n，升到第n-1个阶梯；如果是最后一个重迭峰，则阶梯应即刻返回为零阶梯；然后对阶梯叠加后的波谱作微分，这样第二个峰的叠加谱线在微分后基线回到第一阶梯位置，标注为基线1，同理标注出剩余基线，将基线连起来构成阶梯叠加式傅里叶变换微分谱图。

进一步的，所述步骤s3具体为：

设第k个重迭峰的峰顶点位于s(ωk)，阶梯叠加后变成2ks(ωk)，用相乘函数的微分法则：

式中δ(2k)/δω＝2k–2(k–1)＝2，等于与前面k–1阶的差值；对于后面所有同阶的峰值s(ω)，2k是一个常数，即δ(2k)＝0：

(5)和(6)表示除峰顶点外，通过微分运算波谱其余各点都被处理成谱图的本底。

本发明与现有技术相比具有以下有益效果：

1、本发明不但保持叠加式傅里叶变换法增强峰值一倍的效果，而且分辨率有效提高，同时借助于微分波谱的特性，提升了傅里叶变换的分析灵敏度和反卷积分辨能力。

2、本发明可以应用于用干涉仪和傅里叶变换技术的红外光谱仪器设备，磁共振成像仪器设备，对多维谱图和成像可分别对各维做同样的处理，从而快速实现倍增原有谱图和成像的辨识度。

附图说明

图1为本发明一实施例中常规傅里叶变换波谱出现重迭峰的示意图；

图2为图1中的重迭峰在实施叠加式傅里叶变换后产生的倍增分辨率和峰强度示意图；

图3为本发明的阶梯叠加原理示意图；

图4为图1中的重迭峰在实施本发明的阶梯叠加式傅里叶变换微分法后产生的高分辨高灵敏度示意图；

图5为本发明一实施例中采用本发明的阶梯叠加式傅里叶变换微分技术以及常规傅里叶变换技术获取的聚苯乙烯指纹区波数为1550-1650cm^-1和2970-3200cm^-1波段的红外光谱对比；

图6为本发明一实施例中用本发明的阶梯叠加式傅里叶变换微分技术以及常规傅里叶变换技术获取的乙基苯的苯环氢核磁共振谱图对比。

图7为本发明一实施例中用3特斯拉磁场获取的人造膜磁共振成像(像素512x512)；

图8为图7的原始成像图实施阶梯叠加和微分曲线的示意图；

图9为本发明一实施例中对图7的人造膜实施本发明的阶梯叠加式傅里叶变换微分技术后的成像；

图10为本发明一实施例中1.5特斯拉磁场获取的常规傅里叶变换磁共振头部医学成像；

图11为本发明一实施例中对图10的头部医学成像实施本发明的阶梯叠加式傅里叶变换微分技术后的成像。

具体实施方式

下面结合附图及实施例对本发明做进一步说明。

本发明提供一种阶梯叠加式傅里叶变换微分方法，包括以下步骤：

步骤s1:对分析信号在有限时间域进行采样，并通过离散化处理转换成谱图，通过其微分曲线鉴别谱峰顶点位置；

步骤s2:根据叠加式傅里叶变换原理对谱图中的各峰顶点以阶梯形式叠加；

步骤s3:对阶梯叠加后的谱图做微分，进一步获得阶梯叠加式傅里叶变换微分波谱或图像。

在本实施例中，参考图3，用四个重迭峰f1、f2、f3和f4(图下方的灰色曲线)示意阶梯叠加式微分方法的步骤；首个峰f1定位在零阶梯，它的左半峰归于零基线(标记为基线0)，右半峰依据右叠加函数，峰值被乘以2；由于峰是重迭的，第二个峰f2的左半个峰也被乘以了2定位在第一级台阶(如果这四个峰重迭的太拥挤，峰f3甚至f4的前面重迭部分都会被带到第一级台阶)，右半部分的峰强度被乘以4，同时带领第三个峰f3的前半部分抵达至第二级台阶；依此递增步骤将第三个峰f3的右半部分乘以6，升到第二级台阶，第四个峰f4的右半部分乘以8，升到第三级台阶。如果是最后一个重迭峰，则阶梯值应即刻返回基线0。然后对阶梯叠加后的波谱作微分，这样第二个峰的叠加谱线在微分后基线按图3中的箭头方向回到第一级台阶位置(基线1)；照此继续，第三个峰基线回到第二级台阶位置，标注为基线2；第四个峰基线回到第三级台阶位置，标注为基线3。然后把这些基线连起来就构成阶梯叠加式傅里叶变换微分谱图：即各谱峰顶点采取阶梯叠加微分，而基线回归于常规的微分谱。与单纯叠加式傅里叶变换技术相比依然保持了能把峰强度增加一倍的特点，而且采用微分大大提高了谱图的分辨率和灵敏度。

优选的，在本实施例中，设一个时间域t的分析信号s(t)被截取到n个分立的信号[s(0)、s(1)、s(2)、···、s(k)、···、s(n-1)]；通过傅里叶变换得出一组n个频率域数据s(ω)＝[s(0)、s(1)、s(2)、···、s(k)、···、s(n-1)]，后者用逆傅里叶变换可以变回时间信号s(t)，数学上用矩阵方式简单地表达为：

傅里叶变换波谱和成像有几个非常好的特性：适宜于傅里叶变换的分析信号都是周期变化的比较单纯的信号，例如红外光谱仪用激光光源产生的干涉信号，核磁共振波谱仪由射频波激发的核自旋自由感应衰减信号，磁共振医学成像专注测量人体组织内水分子的氢核磁共振信号等等；傅里叶变换波谱的峰宽跟信号的频率没有关联，主要取决于分析信号的测量时间和衰减系数，所以在同等环境条件下，除了峰高外，所有峰都应具有统一的峰形，例如谐振衰减信号，它们的傅里叶变换波形都是洛伦兹形。

对傅里叶变换波谱s(ω)乘上一个对角矩阵diag[ak]：

式中如果有m个重迭在一块的峰，对角线元素ak取值为2n，n＝0，1，2，···，m；

阶梯叠加对角矩阵进一步写成：

通过微分学二阶导数的鉴别(参见图1)，很容易就定位出谱图曲线的峰顶点ω1、ω2、ω3和ω4；根据定位的各峰顶逐级实施阶梯叠加：2×1，2×2，2×3，和2×4等等。在最后一个重迭峰被阶梯叠加后应尽快返回到零(基线0)，避免它们的重迭信息被携带给后面的不相重迭的峰。由(3)产生的阶梯叠加谱进行微分后，高阶叠加峰的基线恰好回到相邻的低一级台阶。

设第k个重迭峰的峰顶点位于s(ωk)，阶梯叠加后变成2ks(ωk)，用相乘函数的微分法则：

式中δ(2k)/δω＝2k–2(k–1)＝2，等于与前面k–1阶的差值。对于后面所有同阶的峰值s(ω)，2k是一个常数，即δ(2k)＝0：

式(5)和(6)表示除峰顶点外，通过微分运算波谱其余各点都被处理成谱图的本底，因而至少使用三个阶梯叠加式傅里叶变换微分数据点就可以准确地定出一个波谱峰。通常离散化傅里叶变换计算让δω＝1，这三个点分别是前一阶梯的最后一个微分谱点2(k–1)δs(ωk-1)，峰顶点2s(ωk)+2kδs(ωk)，以及同阶梯紧跟着峰顶点的微分谱点2kδs(ωk+1)。

参考图4，本实施例中，显示的是四个模拟的重迭峰(底部的灰色谱图)用阶梯叠加式傅里叶变换方法被成功地完全分离开，每个峰的强度都增加了一倍(上部的黑色谱图)。把各个阶梯微分后的基线连接在一起就产生一个新的傅里叶变换波谱图，峰强度被增加一倍，灵敏度和分辨率可提高超过一倍。曲线峰顶点认定，阶梯叠加峰强度和微分计算对计算机来说都是比较简单的程序，因此在常规傅里叶变换波谱的基础上能快速有效地获取阶梯叠加式傅里叶变换微分谱图。图4中的微分谱本底还可以被抹平归于零基线，让谱图显得更整洁好看，但这些微分本底含有谱峰重迭度和藏匿在原始谱峰底下的噪音等有价值的信息，可利用于推算重迭峰的真实峰高；一般计算信噪比都是通过测量最接近谱峰的一段平滑基线的峰-峰值作为本底噪音，而不是真正藏在峰底的噪音，阶梯叠加式傅里叶变换微分技术提供了一种可能的途径测出真实的本底噪音。

波谱学中常用的微分谱法具有高灵敏度和对重迭峰有去卷积分辨能力，本发明的阶梯叠加式傅里叶变换微分技术继续保持了这两个优点，又弥补了通常微分波谱峰形不规则和峰值计算繁杂的不足。理论上，洛伦兹峰形和峰高都可以用半峰高处的峰宽w1/2来表达，设两个傅里叶变换谱峰间距为d，则它们的分辨率rs可简单的写成：

rs＝d/(2w1/2)(方程式7)

以图1接近于洛伦兹峰形的重迭谱带为例，表1列举了适用半定量的快速构建阶梯叠加式傅里叶变换微分波谱和图像的最佳范畴。

表1

傅里叶变换谱峰的定位需要有充分的采样数据支持，半定量分析最好有10倍以上奈奎斯特采样标准(nyquistcriterion)的数据点。对于分辨率≥0.5的重迭谱带，借助于通常微分谱法的一阶和二阶导数就可以准确地定位出谱图曲线的凸驻点，即峰顶点；当谱峰之间的分辨率小于0.5但大于0.35时，可借助于通常二阶和三阶导数找出一阶导数的斜率变化点(凸驻点和凹驻点)定位次峰的峰顶点。

当然可以不局限于表一列举的参数，例如通过分析谱峰的不对称性，重迭程度，抑制本底噪声，采样频率大于20倍奈奎斯特采样标准等等手段获取更好的定量结果，这都需要更多的分析和计算时间为代价。

实施例一

本实施例提供一种基于阶梯叠加式傅里叶变换微分技术快速获取高分辨率高灵敏度红外光谱的方法。在本实施例中，采用nicoletprotégé460商品型傅里叶变换红外光谱仪，配置632.8纳米(6.328×10^-5厘米)波长的氦氖激光红外光源。干涉仪光程双方向移动阶数3295点，设置709个波数读取点，以波数间距3.85cm^-1和16cm^-1分辨率获取的聚苯乙烯傅里叶变换红外光谱。图5下层灰色谱线是聚苯乙烯红外光谱指纹区的原始谱图，在波数2970-3200cm^-1区间只显示四个特征峰,分别为2854cm^-1，2924cm^-1，3028cm^-1和3062cm^-1。在波数1550cm^-1-1650cm^-1区间只出现一个特征峰，波数1601cm^-1。通过采用本发明的阶梯(右)叠加式傅里叶变换微分技术后，如图5上层黑色谱线所示，在波数2970-3200cm^-1区间显示出八个特征峰，分别为2850.3cm^-1、2908.2cm^-1、2942.9cm^-1、3004.6cm^-1、3027.7cm^-1、3062.5cm^-1、3085.6cm^-1和3108.7cm^-1。在波数1550-1650cm^-1区间轻易就分辨出两个特征峰，1581.4cm^-1和1600.7cm^-1。图5中的波数2850.3cm^-1、2908.2cm^-1和2942.9cm^-1峰旁边伴随的微分隆包表明这三个峰(或许还有更多的峰)重迭的比较严重。在分辨率为16cm^-1的原始数据基础上采用本发明中的阶梯叠加式傅里叶变换微分技术得到红外光谱完全可以与高分辨率4cm^-1聚苯乙烯傅里叶变换红外光谱图相媲美。

实施例二

本实施例提供一种基于阶梯叠加式傅里叶变换微分技术快速获取高分辨率高灵敏度核磁共振谱的方法。图6对比从美国通用qe300核磁共振谱仪获取的化合物乙基苯的300兆赫兹氢核(¹h)核磁共振谱的前后效果。主要仪器和工作参数为：磁场强度7特斯拉，驻留时间250微秒，扫描带宽4000赫兹，偏置频率1850赫兹，采样时间0.512秒，采样点间距1.95赫兹，通过对核磁共振双通道自由感应衰减信号各添加一倍的零，把原始数据点从2048个扩充到8192个数据点。图6底部灰色实线显示的是乙基苯的苯环氢300赫兹常规傅里叶变换核磁共振谱，波谱范围从7到7.4ppm，主要分布在两个谱带：苯基邻位&对位氢谱带(7.05到7.15ppm)和苯基间位氢谱带(7.16到7.26ppm)。邻位&对位氢谱带可勉强辨别出拥挤在一堆的五个峰；间位氢谱带只出现三个峰。图6上部的黑色实线谱是用本发明的阶梯(右)叠加式傅里叶变换微分技术后产生的阶梯叠加式傅里叶变换微分核磁共振谱，本发明的不但彻底分辨开邻位&对位氢谱带的五个重迭峰，而且把干扰严重的溶剂峰从间位氢谱带里完全分离出来。阶梯叠加式傅里叶变换微分技术从300兆赫兹核磁共振谱仪获取的谱图，无论分辨率还是灵敏度都丝毫不逊色于从500兆赫兹核磁共振波谱仪获取的乙基苯谱图。

实施例三

本实施例中，傅里叶变换磁共振成像虽然基于二维或三维空间核磁共振原理，但它与上述的傅里叶变换核磁共振波谱有几个显著不同之处。

(1)图像是按照测量时间顺序来定像素的空间位置，名义上的时间域分析信号；

(2)磁共振成像测量同一种原子核(如水分子的氢核)在三维梯度磁场产生的核自旋共振频率，这些频率信号以梯度零为中心排列存储的数据称为k-空间，名义上属于频率域空间；

(3)所以磁共振成像是从k-空间用逆傅里叶变换转换而成，但逆傅里叶变换与(正)傅里叶变换有很多共性，如处理时间域信号的添零，插值及截趾等常用技术依然对k-空间有效；

(4)磁共振成像被测身体部位的每个像素都是‘峰顶点’(身体部位之外的像素才是本底)，所以每一维的像素就相当于肩并肩的重迭‘峰’。

为了实施本发明提出的阶梯叠加方法，磁共振成像的每个像素都需要在各自身旁额外提供一个像素点以便形成台阶。可以利用在本发明实施例二中简单的添零技术或者更复杂的总变量约束数据外延法(totalvariationconstraineddataextrapolation)将k-空间规模扩充一倍或一倍以上。

图7是从西门子verio3ttim核磁共振成像仪获取的人造膜成像，像素为512x512。主要仪器和工作参数为：磁场3特斯拉，自旋回波脉冲序列，层厚4.0毫米，驻留时间15.6微秒，脉冲重复时间600毫秒，回波时间6毫秒，像素带宽250赫兹。从西门子rhp软件获取的该图像显示有重影。在其k-空间添零扩充为1024x1024矩阵并校正重影后，即可实施阶梯叠加式傅里叶变换微分技术。

图8显示其中一条y轴1024个像素的阶梯(右)叠加曲线(灰色细实线，以左侧纵坐标作图)，和微分像素曲线(黑色粗实线,以右侧纵坐标作图)，图8的初始台阶不是0阶，而是直接从22台阶开始，这样一个简单的操作节省了定位像素边界在微分时造就的‘负’峰及抹平的操作时间，因为数字成像的计算量远远大于波谱作图。图9是用本发明的技术制作的人造膜成像，忠于原始数据，消除了重影，真实地显示出人造膜的平滑表面。

图10是用1.5特斯拉磁场获取的头部医学磁共振成像，主要工作参数：层厚6.0毫米，驻留时间7.2微秒，脉冲重复时间480毫秒，回波时间11.7毫秒，像素带宽93赫兹。图11显示使用了本发明的阶梯叠加式傅里叶变换微分技术后，把图10原始图像的辨识度提高到接近于3特斯拉磁场成像水平。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图1个流程或多个流程和/或方框图1个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图1个流程或多个流程和/或方框图1个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图1个流程或多个流程和/或方框图1个方框或多个方框中指定的功能的步骤。

以上所述，仅是本发明的较佳实施例而已，并非是对本发明作其它形式的限制，任何熟悉本专业的技术人员可能利用上述揭示的技术内容加以变更或改型为等同变化的等效实施例。但是凡是未脱离本发明技术方案内容,依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与改型，仍属于本发明技术方案的保护范围。