快速傅里叶变换算法-盘古文库

快速傅里叶变换算法（精选10篇）

快速傅里叶变换算法第1篇

电力系统频率测量的实质是信号观测模型的动态参数辨识问题,即利用真实物理信号输入,通过一定的信号处理和数值分析过程,实现对预定模型参数的较好估计。精确的频率测量是电力系统运行与多种控制调节的基础,并从根本上保证了其他电气量计算的精度。假设系统中仅含有基波分量,其额定频率为f0,由于系统的真实频率f未知,因此只能根据f0进行采样。若用 Δf表示频差,则真实频率f可表示为f = f0+ Δf。所有的测频原理都是基于如何确定 Δf的值[1]。

目前,国内外学者提出了很多检测频率的理论和算法,可主要分为以下几个方面: 基于正弦信号模型的检测算法、周期法、随机模型算法、基于周期信号模型的算法、现代谱分析和基于HHT变换的频率检测算法。此外,二次型商法、虚拟转子法、正交信号法、人工智能技术( 人工神经网络、模糊逻辑、专家系统和遗传算法) 等算法也逐渐在电力系统频率测量领域得到应用或初步探索[2]。上述理论和算法在频率测量的精度、计算量、响应速度、抑制谐波和噪声抗干扰能力以及在实现难易程度等各方面各有优劣。

上述算法的主要矛盾是实时性、准确性与快速性不统一。有些算法虽然易于实现但其实时性与准确度低,如传统的周期法; 有些算法为了提高测量的准确度而牺牲了计算量,实时性不佳,如解析法; 有些算法虽然具有较强的抗干扰能力但计算量过大,从而影响了测量的实时性,如误差最小化类算法、人工神经网络法; 有些算法虽然能够满足测量的快速性,但抗干扰差,计算量偏大,如正交去调制法。DFT( FFT) 类算法虽然具有一定的时滞,但其改进算法却提高了其测量范围、精度和稳定性,并且不受谐波的干扰[3]。

1 频率偏差算法基础理论

传统频率测量算法是使用傅里叶变换,但是该算法有一个比较明显的缺陷,就是由于傅里叶变换后,频域分辨率与信号时间成反比关系,即用来分析的信号时长越长,频域分辨率也越高,例如当采样信号为时间长度为1s,则频域分辨率为1Hz,若信号时间长度为10s,则频域分辨率为0. 1Hz。而由于本文系统是在线分析的系统,具有较高的实时性要求,采样时长约为0. 1s,则频域分辨率为10Hz,显然,传统的傅里叶变换很难保证系统的测量精度和实时性要求,需要对傅里叶算法进行改进。

文献[1]中提到了一种被广泛采用的递推傅式算法,但该算法需要计算两个相邻周期的相位来计算准确的相位差值,因此需要根据频率偏差对算法进行迭代修正。文献[1]和[4]中提到了采用一种改进的FFT算法测频,具有很好的测量精度和实时性,本文即采用此算法作为基础算法,并在此基础上提出了一种改进算法。算法原理如下。

以理想正弦电压信号( 其中Um为电压有效值,θ0为初始相位)为输入信号,由于实际角频率 ω 在额定角频率 ω0上下波动,在 ω 未知的情况下,以额定角频率进行近似计算,近似傅里叶变换得到[4]:

联立式( 1) 和式( 2) ,可以得到:

由tan( πΔf T0+ θ0) cot( πΔf T0+ θ0) = 1 ,化简可得基于傅里叶变换的频率测量算法的基本公式:

由于余弦函数为偶函数,因此,式( 4) 中 Δf的符号完全取决于角度 πΔ2f T0的符号,由式( 1 )可得:

若- 12. 5Hz < Δf < 12. 5Hz ,对于50Hz工频电压、电流来说,已经满足测量需求,此时,2πΔf T0的取值范围为( - π/2,π/2) ,并满足下列条件:

由此可知,a的符号将完全取决于cos( πΔf T0+θ0) ,且b的符号则取决于sin( πΔf T0+ θ0) ,同理,a' 的符号完全取决于cos( 3πΔf T0+ θ0) ,b' 的符号取决于sin( 3πΔf T0+ θ0) 。因此可以依据a 、b的符号确定 πΔf T0+ θ0在四象限中的哪一象限,记为P象限,同理,可以依据a' 、b' 的符号确定3πΔf T0+θ0在四象限中的哪一象限,记为Q象限。由于( 3πΔf T0+ θ0) - ( πΔf T0+ θ0) = 2πΔf T0,且2πΔf T0∈ ( - π /2,π /2) ,因此P和Q为同一象限或者相邻象限[1]。

当 πΔf T0+ θ0和3πΔf T0+ θ0位于同一象限,也即P = Q时,可以通过三角函数的单调增减性判断角度增量2πΔf T0的符号,从而得到 Δf的符号。

①若P = Q = 1 或P = Q = 2 ,余弦函数为单调递减函数,此时若a > a' ,则2πΔf T0> 0 ,反之则为负。

②若P = Q = 3 或P = Q = 4 ,余弦函数为单调递增函数,此时若a < a' ,则2πΔf T0> 0 ,反之则为负。

当 πΔf T0+ θ0和3πΔf T0+ θ0位于不同象限,也即P ≠ Q时,由于P和Q为相邻象限,其相互关系可以分为P、Q为顺时针排列或逆时针排列。

①当P、Q为逆时针排列,也即Q = P + 1 或Q = P - 3 时,角度增量2πΔf T0为正。

②当P、Q为顺时针排列,也即Q = P - 1 或Q = P + 3 时,角度增量2πΔf T0为正。

由上述方法判断得到 Δf的符号,并将其带入式( 4) ,即可得到待测信号的精确频率。由于式( 4) 不需要计算相邻两个周期的信号相位得到精确的相位差值,因此无需根据频率偏差值对算法进行迭代修正,即可得到系统的频率f的精确值。

2 频率偏差算法改进

上述算法虽然精度高,实时性好,但算法较为复杂,特别是在判断频率偏差 Δf的符号时,逻辑复杂,不易编程实现,为此,本文提出了一种改进算法,大大简化了对频率偏差符号判断的逻辑。原理如下:

由式( 1) 可得:

由式(7)且可以得到tan(πΔf T0+θ0)与b/ a同号。

由式( 3) ,可得:

因此,当cos(2πΔf T0)-a'/a>0时,sin(2πΔf T0)与tan(πΔf T0+θ0)同号,即与b/ a同号。当cos(2πΔf T0)-a'/a<0时,sin(2πΔf T0)与b/ a异号。由2πΔf T0∈(-π/2,π/2),则Δf与sin(2πΔf T0)同号,由此即可判断频率偏差Δf的符号。

该算法与原算法相比,其判断逻辑大大简化了。使用程序实现起来也非常方便,例如使用MATLAB程序实现如下:

程序中p即为由式( 4 ) - ( 8 ) 计算得到的cos( 2πΔf T0) 的值,m即为 Δf的符号。

使用该算法,程序大大简化,仅仅需要4行程序,而原算法的实现需要40余行程序。使用MAT-LAB仿真,得到原算法和改进算法的时间占用如图1所示,改进算法的时间占用较原算法的时间占用减少了一半多,这对于算法实时性的提高具有实际意义。

该方法的计算周期仅为两个工频周期,即0. 04s,具有很好的实时性,另一方面,该算法在50Hz附近又具有较高的测量精度,满足本系统的设计要求。虽然该方法是针对仅含基波分量的模型所设计的,但因风电信号中基波分量仍是占据绝对的主导,经过MATLAB的理论仿真和试验实测,此时该算法仍然适用,能够保证足够的测量精度和很好的实时性要求。

3 频率偏差算法的MATLAB仿真

使用MATLAB对上文所述基于傅里叶变换的频率改进算法进行编程测试。选取含有基波以及2 ~ 5 次谐波分量的电压信号:

对算法进行仿真分析,在49. 5Hz ~ 50. 5Hz频率范围内均匀选取测试信号。所得结果如表1所示。

由表1可知,该频率测量算法在50Hz附近具有极高的测量精度,绝对误差不超过0.0005Hz,相对误差更是低至十万分之一以下。

4 结束语

传统频率测量算法是使用傅里叶变换,本文提出了一种新型的频率偏差符号判断算法,相对原有算法,本算法大大简化了判断逻辑,减少了计算时间,提高了算法的实时性,精度也较好。

本文所设计的算法是基于固定采样率的非同步采样,这种测量方式会产生频谱泄露,虽然本文所设计算法已经具有了足够高的频率测量精度,但仍有改进的余地。

摘要：频率测量是电能质量测量的重要基础,传统频率测量算法是使用傅里叶变换,很难保证系统的测量精度和实时性要求,需要对傅里叶算法进行改进。文中介绍了一种基于傅里叶变换的频率偏差的精确测量算法,并提出了一种新型的频率偏差符号判断算法。该算法大大简化了判断逻辑,减少了计算时间,提高了算法的实时性,经过MATLAB仿真,精度也较好。

关键词：频率偏差,傅里叶变换,实时性

参考文献

[1]束洪春.电力工程信号处理应用[M].北京:科学出版社,2009.

[2]欧立权.电力系统频率测量方法及应用的研究.[D].长沙:湖南大学,2007.

[3]李璟.基于DSP的电力系统频率检测的研究[D].福州:福州大学,2003.

快速傅里叶变换算法第2篇

2、将后面的正弦函数展开：

于是得到：

那么如何计算an,bn,a0这些参数成为能否展开成为正余弦函数的关键。

上面的这些积分为0被称之为正余弦函数的正交性。这些证明很简单，可惜当初学习正余弦函数的时候可能遇到过，但是却不知道这些东西能干什么用。下面的处理手段凸显了大师的风范：

如果我们队原函数进行如下积分，得到很神奇的东西：

后面的积分很明显是0，于是我们求出了a0的值。

那么如何求出an,如果让原函数乘以cos(nx)再进行积分。

利用三角函数的正交性，可以得到：

再用sin(nx)乘，再进行积分就会得到bn,于是乎得到了一个任意函数展开成为正余弦函数的通用表达式，同时为什么会出现A0/2而不是直接的A0的原因也很明朗：就是让整个表达式更具有通用性，体现一种简洁的美。

通过了以上的证明过程，应该很容易记住傅里叶变换的公式。

到此为止，作为一个工程人员不用再去考虑了，可是作为每一个数学家他们想的很多，他们需要知道右侧的展开式为什么收敛于原函数，这个好难，有个叫Dirichlet的家伙证明出如下结论：

有兴趣的可以继续找书看，可惜我有兴趣没时间····

至此以2π为周期的傅里叶变换证明完毕，只不过我们经常遇到的周期函数我想应该不会这么凑巧是2π，于是乎任意的一个周期函数如何知道其傅里叶变换呢，数学向来都是一个很具有条理性的东西，任意周期的函数的傅里叶变换肯定也是建立在2π周期函数的基础之上的。

也就是说如何让一个以2l为周期的函数变成一个以2π为周期的函数，于是乎可以使用z=2π*x/(2l),这样就z就是一个以2π为周期的函数了，于是乎得到如下公式：

傅里叶函数看起来其实还是比较复杂的，有没有一种更简单的表达形式来表示呢。既然提出这个问题，肯定是有的，我个人猜想肯定是复变函数大师在挖掘复变函数的时候，用复变函数去套用经典的傅里叶变换，偶然间发现的······

一个基本的欧拉公式eiθ=cosθ +i*sinθ,这个很容易可以从复数的几何意义上得知，我们通过取两个互为相反数的θ可以得到两个式子，进而可以得到cos 和 sin 的复数

表达形式：



fT(t)c0



[cne

n1

jnt

cne

jnt

]

即：fT(t)



n

cne

快速傅里叶变换算法第3篇

通过数字方法对图像进行加密，解密过程用光学装置实现。计算机模拟结果表明，该加密方法解密图像质量好，系统安全性良好。

关键词：加密；傅里叶变换； Gyrator变换；相位密钥

中图分类号： O 438 文献标志码： A doi： 10.3969/j.issn.1005-5630.2015.01.016

Abstract：A novel optical encryption is proposed based on Fourier transform and Gyrator transform. The original image is firstly multiplied with the first random phase function. Then the Fourier transform is performed and the information of the frequency domain minus the second random phase function. A complex function is obtained and the Gyrator transform is performed to get the encrypted image. The transform angle of the Gyrator transform serves as a key and the second random phase function serves as a phase key. The two keys ensure the security of the encryption system. The encryption process is performed digitally and the decryption process can be implemented in the digital and optoelectronic way. The optical decryption scheme is designed. Computer simulations indicate that the decrypted image is of good quality and the encryption system has a high security.

Keywords：encryption； Fourier transform； Gyrator transform； phase key

引言

近年来，光学图像加密技术因其高处理速度、高并行度，可以提供相位、振幅、偏振态、波长、空间频率等多种加密自由度而引起了人们的极大兴趣。Refregier等提出双随机相位编码方法[1]，在空域和频域对加密对象进行编码，并将输出平面得到的复振幅图像作为加密结果。双随机相位编码技术是光学信息加密领域中的经典技术，具有较高的安全性和鲁棒性，但也存在一些缺点，如系统对精度要求高，容偏能力低，加密结果为复振幅，不便于记录和传输等。研究人员深入研究了双随机相位编码方法并提出了许多改进的或者是新的加密方案。2004年，Situ等根据菲涅耳夫琅和费衍射，提出了菲涅耳域双随机相位编码方法[2]，利用两次菲涅耳变换和两个随机相位模板对原图像进行双随机相位编码加密。该系统的光学加密解密装置和4f系统类似，但是不需要透镜，并且增加了密钥空间。分数傅里叶变换和Gyrator变换都是传统傅里叶变换的推广形式。基于分数傅里叶变换，研究人员提出了相应的双随机相位编码方法[3]。其加密解密过程并没有增加光学元件，但增加了密钥空间，大大提高了系统的安全性。1993年，Lohmann深入研究光学分数傅里叶变换并给出了两种实现分数傅里叶变换的光学装置[4]。2007年，Rodrigo等给出了Gyrator变换的积分定义形式[5]，并设计了一种Gyrator变换的光学实现装置[6]。Gyrator变换的变换阶数可用作密钥，并且其变换阶数的改变也是通过绕光轴旋转透镜实现的，易于调节并减少了光轴校准带来的误差。

本文基于Gyrator变换提出了一种新的加密方法。原图像先经第一个随机相位函数扰乱后，再将其频域信息归一化后直接减去第二个随机相位函数，所得的频域信息经Gyrator变换后得到加密图像。该加密方法中，第二个随机相位函数作为密钥，Gyrator变换的变换角度作为解密密钥。该加密方法、加密过程直观，有较高的安全性。

1 理论分析

1.1 Gyrator变换

加密过程用数字的方法实现，解密过程可以用数字的方法实现也可以用光电混合系统实现，实现装置如图2所示。图中，SLM1和SLM2为空间光调制器，GT为Gyrator变换的光学装置，半透半反棱镜、全反射镜、傅里叶透镜组成逆傅里叶变换的光学装置。

将SLM1加载为全息图E′（x，y），置于Gyrator变换的输入平面，经离轴的参考光照射后做变换角度为-α的Gyrator变换得到恢复的g（u，v）。将SLM1加载为随机相位函数R（u，v），R（u，v）与恢复的g（u，v）进行相关叠加。将叠加结果置于逆傅立叶变换的输入平面，最后，在逆傅里叶变换的输出平面用CCD接收得到解密图像。

2 仿真实验及性能分析

用MATLAB 7验证方法的正确性和有效性。原图像“Babara”像素为512×512，灰度值为0到256，如图3（a）所示。加密过程中，原图像空域信息被扰乱后，其频域信息减去随机相位函数R（u，v），得到复函数g（u，v），图3（b）为R（u，v）。取Gyrator变换的变换角度α为1.2，则得到的加密图像如图3（c）所示。解密过程中，各个密钥都正确时，恢复所得的g（u，v）的实数部分如图3（d）所示，R（u，v）与恢复的g（u，v）的叠加结果的实数部分如图3（e）所示，在逆傅里叶的输出平面，用CCD接收到的解密图像如图3（f）所示，可以看出，原始信息完全恢复。

nlc202309021303

为研究各个密钥对系统抗攻击性能的影响，图3（g）给出了相位密钥R（u，v）错误时的解密结果，MSE为7.541×104，解密图像为噪声分布，无法得到原图像信息。图3（h）为Gyrator变换的变换角度α错误时的解密结果，MSE为7.612×104，无法得到原图像信息。通过直接观察及对MSE值分析可得：如果密钥R（u，v）和Gyrator变换的变换角度α有一个不正确，将无法恢复原图像，说明该方法具有较高的安全性。

为了通过原图像和解密图像的MSE值与Gyrator变换的变换角度α的偏离之间关系，研究该加密方法对α的敏感性。为此，改变α，使其变化区间为[1.1，1.3]，步长为0.005。原图像和解密图像的MSE值随α的变化曲线如图4所示。可以看出，稍微变换一点，就无法得到正确的解密图像，说明该加密系统对α的变化很敏感，具有较高的安全性。

3 结论

本文基于傅里叶变换和Gyrator变换提出了一种图像加密方法。原图像的频域信息减去一个随机相位函数后经Gyrator变换得到加密图像，密钥为随机相位函数及Gyrator变换的变换角度，并设计了解密过程的光学装置。用计算机模拟得到了加密图像，并验证了所有参数正确时，才可以完全地恢复原图像。对于非法用户，由于没有正确的密钥，无法得到正确的解密图像。研究表明，解密图像对Gyrator变换的变换角度的敏感性很高，变换角度变化一点，就无法得到正确的解密图像。该加密系统具有良好的安全性。

参考文献：

[1] REFREGIER P，JAVIDI B.Optical image encryption based on input plane and Fourier plane random encoding[J].Optics Letters，1995，20（7）：767-769.

[2] SITU G H，ZHANG J J.Double random-phase encoding in the Fresnel domain[J].Optics Letters，2004，29（14）：1584-1586.

[3] UNNIKRISHNAN G，JOSEPH J，SINGH K.Optical encryption by double random phase encoding in the fractional Fourier domain[J].Optics Letters，2000，25（12）：887-889.

[4] LOHMANN A W.Image rotation，Wigner rotation，and the fractional Fourier transform[J].Journal of the Optical Society of America A，1993，10（10）：2181-2186.

[5] RODRIGO J A，ALIEVA T，CALVO M L.Gyrator transform：properties and applications[J].Optics Express，2007，15（5）：2190-2203.

[6] RODRIGO J A，ALIEVA T，CALVO M L.Experimental implementation of the gyrator transform[J].Journal of the Optical Society of America A，2007，24（10）：3135-3139.

（编辑：刘铁英）

快速傅里叶变换算法第4篇

在各种场合的扩音工作中,音箱出现正反馈是啸叫的主要成因。啸叫的出现,轻则破坏现场气氛,影响演出效果,严重的可能损坏扩音设备,所以,如何解决现场的啸叫,是扩音工作的一大难题。数字信号处理技术主要研究用数字序列来测试信号,并用数学公式和运算来对这些数字序列进行分析和处理,包括有数字滤波,波形分析,幅值分析和频谱分析等。FFT[1](Fast Fourier Transformation),即为快速傅里叶变换,是离散傅里叶变换的快速算法,它是根据离散傅里叶变换的奇、偶、虚、实等特性,对离散傅里叶变换的算法进行改进获得的。本文拟采用FFT算法检测音频信号的频谱,从而结合滤波网络实现对啸叫频率的检测和抑制。

1 系统方案设计

总体系统设计方案如图1所示。将麦克风音频信号送入拾音电路,而后经通带为200Hz~10KHz的带通滤波器滤除干扰信号。当啸叫产生时,单片机采集部分信号从而经FFT算法检测啸叫发生频率,然后经低通滤波器滤除易产生啸叫的频段;若无啸叫发生,则直接将信号送入功率放大器放大,最终送喇叭输出[2]。

1.1 拾音电路模块

采用TI公司生产的双运放LM358通过外接电阻电路构成一个能放大多倍的拾音电路,运放的输出端接入压控放大器VCA822。LM358是一款内部包括两个独立的、高增益、内部频率补偿的双运算放大器芯片。使用LM358搭建拾音电路设计容易、制作简单且成本低,其放大倍数可以满足设计要求。

1.2 可控增益放大器的选择

VCA822是V/V线性连续可变增益放大器,双列14脚封装,可程控输出信号的幅度。其可调增益带宽为150MHz,带宽更大使其峰峰值在0~6V之间连续可调。

1.3 啸叫检测与抑制方案选择

将部分信号经单片机采集,从而通过编程实现FFT快速傅里叶变换算法得到功放模块输出信号的频谱。这样即可直观地测量啸叫频率并使用滤波器滤除啸叫频率达到抑制啸叫的目的。而且,数字滤波的算法和技术较之于模拟滤波更成熟计算结果也更加精确。

2 系统理论分析与计算

2.1 啸叫频率的分析与计算

功率放大器的输出波形经过AD采样获得U(t),而后通过以下程序实现傅里叶变换[3]:

S(t) .real= U(t) = ADC(ch) * 3.3 / 4096;

S(t).imge = 0;

FFT(S(t));

Vol(t) = S(t).real;

可求得啸叫频率:

f= T(Vol(t).max);

3.2啸叫功率的分析与计算

啸叫平均功率由式(1)计算可得:

4 系统硬件电路设计

4.1 拾音电路设计

拾音电路如图2所示,使用Multisim仿真电路麦克风可以看成一个电流源,所以在仿真时可以用电流源I1替代,使用双运放LM358进行设计,一个运放用来把电流变换成电压,第二个运放可调节增益,提高功放的输出电路。单电源设计,给运放一个偏置电压,从而供电由双电源变成单电源。

4.2 频率200Hz~10KHz带通滤波器[4]设计

使用Filter Pro滤波电路设计软件设计4阶巴特沃思带通滤波器,然后使用Multisim软件进行仿真验证,带通滤波器的性能指标符合本系统要求。滤波器电路如图3所示:

4.3 噪声抑制电路设计

在一定环境条件下,啸声的频率范围一般在10KHz以上,所以设计一个4阶低通滤波器,滤除易出现啸叫的高频信号,在一定程度上达到啸声抑制功能。滤波器电路如图4所示。

综上,该系统各个模块的设计方案已论证完毕,下面给出整体系统的原理图(如图5所示)。

5 系统软件程序设计

软件部分可以实现设置键盘和液晶显示的功能。键盘可以设置增益、换页以及启动程序运行。程序设计框图如图5所示,设计思路:首先通过ADC采集电压值,经过FFT变换,求出基波频率,即为输入的频率。当产生啸叫时,频率单一稳定,此时的基波频率即为啸叫频率。通过对电压积分求出有效值,负载确定时可以求出负载功率,即为输出功率。

6 系统测试方案与部分测试结果

6.1 测试方案

硬件部分测试:各模块上电之后,输入有效信号,检测输出信号是否正确,以及按键液晶是否正常工作。

硬件软件联调:通过按键改变输入有效信号,测量输出信号,液晶显示各项数据指标。

6.2 测试方法和测试仪器

测试方法:正常室温环境下,接入输入信号,然后接通电源,在电路正常工作的前提下,测量输出信号。

测试仪器:数字示波器Agilent(DSO-X,2002A,70MHz)、F120型数字合成函数/任意波信号发生器/计数器、学生电源(MODEL, MPS-3005L-3)、乐儿飞M-01USB喇叭、QG-500MIC全向麦克风、VICTOR,VC890+数字万用表。

6.3 测试结果和分析

表2中给出了输出信号的频谱。从数据中可以看出,经带通滤波器滤波,输出信号的有效频带约在500Hz~10KHz的范围内,在通带内信号有较好的频率特性。

测试结论:测试结果显示,该系统可程控设置功率放大器输出功率在50m W~5W的范围内;功率放大器的频率响应范围为:500Hz~10k Hz;当啸叫发生时,低通滤波模块能有效地滤除啸叫发生的高频段信号,从而实现了实时检测和抑制啸叫的功能。

7 结束语

分数阶傅里叶变换的数值实现第5篇

分数阶傅里叶变换的数值实现

信号及其傅里叶变换可以分别反映信号在时频两域内的信息.傅里叶变换是一种常用的数学工具,在数学、物理及工程技术领域都得到了十分广泛的应用.介绍了一种崭新的信号分析工具--分数阶傅里叶变换,并用经典的傅里叶变换的观点对分数阶傅里叶变换进行了解释.对于分数阶傅里叶变换的实现,因一般情况下分数阶傅里叶变换给不出解析表达式,故分数阶傅里叶变换的.数值算法的研究是十分重要的.给出了分数阶傅里叶变换的较准确的数值计算方法.利用此方法对被线性调频函数污染混叠的高斯信号进行了滤波分离.

作者：陈明杰作者单位：重庆工商大学,计算机与信息工程学院,重庆,400067 刊名：重庆大学学报(自然科学版) ISTIC EI PKU英文刊名：JOURNAL OF CHONGQING UNIVERSITY(NATURAL SCIENCE EDITION) 年，卷(期)： 26(5) 分类号：O174.22 关键词：分数阶傅里叶变换傅里叶变换滤波干扰

基于快速傅里叶变换的形状检索第6篇

关键词：图像检索,形状,傅里叶变换

0 引言

本文旨在研究基于形状区域特征的图像检索,提出了一种基于二维快速傅里叶变换的的形状检索算法。第一节介绍了的相关知识,第二节介绍了基于快速傅里叶变换的频域特征提取的具体实现技术,第三节介绍了特征匹配的方法,最后是实验的结果与分析。

1 傅里叶变换

本文算法是在快速傅里叶变换的理论基础上提出的,所以本文首先对傅里叶变换进行了研究。

1.1 傅里叶变换的基本概念

傅里叶变换在数字图像处理中有着广泛的应用,其在数学中的定义非常严格,定义如下:设f(x)为x的函数,若f(x)满足狄里赫拉条件:(1)具有有限个间断点;(2)具有有限个极点;(3)绝对可积。

但图像必须在空间和灰度上都离散化才能被计算机处理,要对数字图像进行傅里叶变换,必须引入离散傅里叶变换(DTF,Discrete Fourier Transform)的概念。

1.1.1 一维离散傅里叶变换

对一个连续f(x)等间隔采样,可得到一个离散序列。设共采了N个样,则这个离散序列可表示为{f(0),f(1),f(2),…,f(N-1)},则其离散傅里叶变换F(u)为:

离散傅里叶反变换为:

根据欧拉公式e±ix=cosx±isinx,式(1)可写成:

所以,一维离散傅里叶函数F(u)也可写成复数函数形式:

其中R(u)和I(u)分别为F(u)的实部和虚部。

上式称为傅里叶变换的幅度或频率谱。

上式称为变换的相角或相位谱。

上式称为功率谱。

1.1.2 二维离散傅里叶变换

在数字图像中,像素点是二维离散的,设一个图像尺寸为M×N(M为图像的高,N为图像的长)的函数为f(x,y),则图像可表示成为一个二维的离散点序列,如式(8)表示:

对于图像的像素点f(x,y),其二维离散傅里叶变换为:

其中:u=0,1,…,M-1;v=0,1,…,N-1

其傅里叶反变换为:

其中:x=0,1,…,M-1;y=0,1,…,N-1

考虑到在数字图像处理中,图像取样一般是方阵,即M=N,则二维离散傅里叶变换公式为:

其中:u,v=0,1,…,N-1

其中:x,y=0,1,…,N-1

与在一维时的情况类似,二维傅里叶变换的频谱、相位角、功率谱如下:

1.2 傅里叶变换的性质

由于我们主要研究的是二维数字图像,而数字化图像都是离散化了的图像,因此下面对二维离散傅里叶变换的一些重要性质进行简述。

(1)线性。离散傅里叶变换是一个线性变换。

(2)分离性。一个二维离散傅里叶变换也可以用二次一维离散傅里叶变换来实现。

(3)对称性。如果离散函数f(x,y)的离散傅里叶变换为F(u,v),那么:I[F(x,y)]=MN·f(-u,-v)(16)

(4)周期性。离散傅里叶变换和傅里叶反变换均是以N为周期。

(5)共轭性。如果离散函数f(x,y)的傅里叶变换为F(u,v),F*(-u,-v)为f(-x,-y)离散傅里叶变换的共轭函数,则它的傅里叶变换具有共轭对称性。其共轭性说明离散函数f(x,y)经过正交变换后得到的F(u,v)是以原点为中心对称的。利用该特性,只要求出半个周期内的值就可以得到整个周期的值。

(6)平移特性。离散傅里叶变换对的平移特性可写成:

式(17)表示将f(x,y)与一个指数相乘就相当于把其变换后的频域中心移动到新的位置。式(18)表明将F(u,v)与一个指数相乘就相当于把其反变换后的空域中心移动到新的位置。

(7)旋转特性。如果空间域离散函数旋转角度θ0,则在变换域中该离散函数的离散傅里叶变换函数也将旋转同样的角度。

(8)尺度变换特性。如果函数f(x,y)的傅里叶变换为F(u,v),给定两个标量a和b,则:

(9)平均值。一幅图像的灰度平均值可由DFT在原点处的值求得,即:

(10)卷积定理。设f(x,y)和g(x,y)分别用尺寸A×B和C×D,周期为M和N的离散数组表示。为了防止卷积后产生重迭误差,需要选择和。

则二维离散傅里叶卷积定理如下:

其中,Fe(u,v)为离散函数fe(x,y)的离散傅里叶变换,Ge(u,v)为离散函数ge(x,y)的离散傅里叶变换。

(11)相关定理。对于离散函数f(x,y)和g(x,y)的相关运算莓,同样必须先对f(x,y)和g(x,y)进行扩展处理,然后再定义相关运算莓如下:

则相关定理如下:

其中Fe(u,v)为函数fe(x,y)的傅里叶变换,Ge(u,v)为函数ge(x,y)的傅里叶变换;Ge*(u,v)为Ge(u,v)的共轭,ge*(x,y)为ge(x,y)的共轭。

1.3 快速傅里叶变换

进行离散傅里叶变换需要的计算量太大,运算时间太长,在某种程度上限制了它的使用。按照式(1),计算一个长度为N的一维离散傅里叶变换,对u的每一个值需要做N次复数乘法和(N-1)次复数加法。那么对N个u,则需要N2次复数乘法和N(N-1)≈N2次复数加法。很显然,当N很大时,计算量是相当可观的。

库里(Cooley)和图基(Tukey)在1965年首先提出一种快速傅里叶变换(FFT)算法,采用该算法进行离散傅里叶变换,复数乘法和加法次数正比于Nlog2N,这在N很大时计算量会大大减少,如表1。

可见,采用FFT可以减少运算量,图像越大减少越多。对于长为1024的离散序列,用普通的离散傅里叶变换往往要计算几十分钟,而采用FFT,一般只要几十秒。快速傅里叶变换是傅里叶变换的一种改进算法,它分析了离散傅里叶变换中重复的计算量,并尽最大的可能使之减少,从而达到快速计算的目的。

在这里,只对一维情况下的快速傅里叶变换进行讨论,因为根据傅里叶变换的分离性可知,二维傅里叶变换可由连续两次一维傅里叶变换得到。

设N为2的整数幂,即:N=2n,并令M为正整数,满足N=2M,则式(1)可写成:

由式(27)可知W2M2um=WMux,所以式(28)可写成:

现在定义:

式(29)可转化为:F(u)=1/2[G(u)+H(u)Wu2M](32)

同样,因为和,式(30)通过式(32)得到:

仔细分析式(30)至式(33)可知,一个N点变换可以通过把原始表达式分成两部分来计算,如式(32)和式(33)所示。计算F(u)的前半部分要对式(30)和式(31)给出的两个N/2点变换进行计算。G(u)和H(u)的计算结果被代入式(32)中得到F(u),u=1,2,…,(N/2-1),另外一半可直接从式(33)得到,而无需另外的变换计算。因此可以将求N个点的离散傅里叶变换转换成求N/2点的离散傅里叶变换,当N为2的整数幂时,式中的G(u)和H(u)可以再被细分成更短的序列,进一步加快运算速度。

1.4 傅里叶变换在图像检索领域的应用

通过上述研究可知,图像的频率是表征图像中灰度变化剧烈程度的指标,是灰度在平面空间上的梯度。如:大面积的沙漠、海水、天空在图像中是一片灰度变化缓慢的区域,对应的频率值很低;而对于地表属性变换剧烈的边缘区域在图像中则是一片灰度变化剧烈的区域,对应的频率值较高。从物理效果看,傅立叶变换是将图像从空间域转换到频率域,其逆变换是将图像从频率域转换到空间域。换句话说,傅立叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数,傅立叶逆变换是将图像的频率分布函数变换为灰度分布函数。

作者认为傅里叶变换具有以下优点:第一,傅里叶变换是在图像的频域对图像进行特征分析,对噪声不敏感,而基于空域的图像分析方法是很难控制噪声对图像带来的影响的;第二,一幅图像由空域转换到频域,使高度相关的空间样值变为相关性较弱的变换系数,减少了空间样值之间的冗余度。正是这些优点,使得傅里叶变换被广范地应用于图像处理和图像分析领域。一维傅里叶变换应用于物体形状的轮廓描述上,它首先要确定图像中物体的轮廓,然后将轮廓映射到复平面上后,再对其进行傅里叶变换。而二维傅里叶变换不需要确定物体的边界,可直接对图像进行变换,对图像分别从x和y方向进行傅里叶变换,应用于对物体形状的区域描述上。

2 基于快速傅里叶变换的特征提取

本文在快速傅里叶变换的理论基础上提出了一种基于区域形状特征的提取方法。该方法首先对图像进行二维的快速傅里叶变换,然后对变换的结果进行分块统计,由每块的频域统计值构成图像的形状特征向量,提取步骤如图1所示。

2.1 图像归一化处理

考虑到图像资源尺寸的不统一,使得傅里叶变换的结果矩阵大小也不统一,将造成后期无法对图像统一提取特征这一问题。作者对图像进行了大小的归一化处理。又考虑到离散快速傅里叶变换仅适用于大小为2n的图像,可以将图像归一化后的尺寸选择为64×64、128×128、256×256、512×512。

归一化的过程实际是对图像进行缩放,缩会使多点合并成为一个点,放会使图像产生空白点。作者采用的具体方法是:对多点合并一个点时,先计算周围8领域的平均值,然后再利用这个平均值作为合并点的值;对产生的空白点需要进行填充处理时,先计算周围8邻域像素点的平均值,然后利用该平均值填充。通过归一化处理的图像,从视觉上看是比较接近的,为后期特征的提取带来了方便。但作者也发现它有可能对图像检索产生一些负面影响,一是如果图像较大,图像压缩则可能丢失一定的图像信息;如果图像较小,这样的拉伸会使图像信息产生一些冗余的附加信息。二是由于原有图像的宽度和高度如果不是1:1的比例,这样的缩放会使图像的形状产生变形和扭曲。

在对图像进行归一化处理时,若原图像的比例不是1:1,则造成了图像的扭曲变形,为图像检索带来了负面影响。同时图像归一化的过小,图像可能会丢失掉一些重要的细节,影响图像的检索效果,而归一化的过大,又增加了傅里叶变换的计算量,影响图像的检索效率。作者权衡以上两个因素,将图像归一化为256×256个像素点。

2.2 对图像进行频域转换

傅里叶变换是将空间域转换为频率域的有效方法,利用傅里叶变换的这种特性就可以将图像的空间信息转换为频率信息进行处理。要提取图像的形状区域信息,则需要对图像进行二维离散傅里叶变换,1.1节中已经详细地介绍了二维图像傅里叶变换的公式,这里就不再重复。作者考虑到二维傅里叶变换需要的计算量依赖于图像的大小,计算处理时间随着图像的增大迅速增长。因此,作者采用了傅里叶变换的快速算法对图像进行二维傅里叶变换,这一做法使傅里叶变换的计算量降至Nlog2N,大大缩减了计算时间,其方法已在1.3节做了详细分析,这里不再重复。快速傅里叶变换虽然可以使计算量下降,但其局限性也是很显然的,它的N要求必须是N=2n,所以作者在2.1节对图像做的归一化处理时就考虑到了这个问题。

由1节可知傅立叶变换是基于复数的,其结果为复数形式。因此,要直接表示一幅图像的二维FFT变换的结果就必须用到两幅图像:一幅表示实部,一幅表示虚部。在频谱图像中,灰色表示的频域值为0,黑色表示的频域值为负值,白色表示的频域值为正值。由图2可以看到4个角上的黑色更黑,白色更白,表示其幅度更大,其实4个角上的系数表示的是图像的低频组成部分,而中心则是图像的高频组成部分,低频部分反映了图像的整体轮廓,高频部分反映了图像的细节。尽管这样表示的频谱图已经能够反映出原始图像的能量分布情况,但是不够清晰、直观,因此,另一种表示方法是引入变换结果的模作为值在频谱图中表示出来,用灰度的明暗来代表模的大小,如图3所示。

研究表明,二维快速傅里叶变换是一种正交变换,变换后信息不会丢失,具有能量保持性,并且能够对能量进行重新分配与集中,使高度相关的空间样值转变为了相关性较弱的变换系数,从而减少空间样值之间的冗余度,降低了对噪声的敏感度。同时,对不同的图像进行二维快速傅里叶变换后得到的频谱图是不同的,且二维快速傅里叶变换的结果直接反映了图像的形状区域信息,因此可以通过对变换后的图像频谱结果进行分析来提取图像的区域形状特征。

2.3 频域特征的提取

对一个N×N的灰度图像进行二维快速傅里叶变换后,结果为一个N×N的复数矩阵,它反映了图像中物体的区域形状信息。通过分析发现,直接从傅里叶系数计算频谱(即模)作为描述符,是不实用的,因为得到的特征不具备旋转不变性。对一个形状描述符来说,是否具有旋转不变是非常重要的,因为近似的形状可能具有不同的方向。作者通过对图2和图3(a)所示结果进行分析后发现,虽然由傅里叶系数实部、傅里叶虚部和系数的模绘制出的频谱图外观上虽然不尽相同,但它们的能量分布却是完全一致的,体现了相同的图像特征;通过实验,作者分析发现对变换的结果仅提取实部作为描述符就可以有效地解决旋转不变性问题,因此作者选择傅里叶系数的实部来描述形状的区域特征,具体做法如下:

由于对图像进行傅里叶变换前已经将图像归一化为256×256,因此这里得到的变换结果为一个256×256的矩阵,如果利用每一个频谱特征构成图像的特征向量,则一幅图像的数据量就达到216,大大增加了图像检索中特征匹配过程的计算量,降低检索效率。在这里,作者把矩阵均匀分割为16×16的块,如图4所示,共提取出256个特征值组成图像的区域形状特征向量。由于每个块中包含了256个傅里叶系数,我们需要对它们进行计算得到该块的特征值,作者选择对每块内的256个系数统计求和再开方,这样做可以获得性能更一致的全方位的能量响应。方法如下:

假设傅里叶系数为:f=A+Bi

其中,Fi为第i个块,Aj为块中第j个傅里叶系数的实部,j的计算公式如下:

其中,r为第i个块所在的行,c为第i个块所在列,m为第i个块中第j个系数所在行,n为第i个块中第j个系数所在列。

由于傅里叶系数实部的值比较大,为了特征匹配后的结果易于理解和分析,作者对特征Fi进行了归一化处理,即

这样就得到了图像频谱的256个特征分量SFi(i=1,2,3,…,N),由它们构成图像的频谱特征向量,即图像的区域形状特征向量,在图像中的特征分量从左到右分别计为1-N。

3 特征匹配

图像检索时,作者用查询图像的频谱特征向量直接和图像库中的每一幅图像的频谱特征向量运用欧式距离法进行相似性比较,得到相似图像,并按相似性进行排序。

式(37)中S(q,d)表示总得相似性距离,SFid表示图像库中被检索的图像的第i个特征分量,SFiq表示待检索图像的第i个特征分量。

4 实验结果与分析

为了更好地评价本章算法的检索性能,本文选择了MPEG-7 Shape B标准形状测试集。测试集中有1400幅图像,被分成70类,每类20幅图像。图像库内的图像均为灰度图像,且每幅图像目标唯一、背景单一,其中目标为白色,背景为黑色。

对性能的评价,本文采用“查准率”对本章算法输出的结果与人们期望结果的一致性进行比较,查准率的公式见式(38)。实验中,作者将P(G)(即输出近似图像的数目)分别选择为10和20。具体实验步骤如下:以MPEG-7Shape B形状库为基础,从中选择10类作为本实验的测试形状库,其中每类包含了20幅图像,这些图像的分类是基于语义的分类,具体分类如图5所示,其中每一类中所包含的图像在视觉上是相似的,一些不同类之间也有相似性,如camel和elephant。对每一类图像随机抽取10幅分别作为关键图像,用查准率公式计算出P10和P20,然后对每类计算平均查准率P軈10和P軈20。

表2为本章算法(记作算法A)和采用傅里叶系数的模作为特征的检索算法(记作算法B)的平均检索准确率。

从表2中的结果可以看出,针对上述各测试图像类,本章算法的平均检索准确率均高于采用傅里叶系数的模作为特征值的算法。由于本章算法在频域提取图像的区域形状特征,故对噪声不敏感,具有平移、尺度不变性;只选用傅里叶系数的实部建立图像的特征向量对图像进行检索,获得了较好的旋转不变性,检索效果如图6,第一副为示例图像。

参考文献

[1]章毓晋.图像处理和分析(图像工程上册).北京:清华大学出版社,2004.

[2]张洁.数字图像边缘检测技术的研究.硕士学位论文,合肥工业大学,2009.

[3]Bober M..MPEG-7visual shape descriptor.IEEE Trans.on Cirrcuits and Systems for Video Technology.2001,11(6):716-719.

[4]冈萨雷斯.数字图像处理.北京:电子工业出版社,2005.

[5]Latecki L.J.Shape Data for the MPEG-7Core Experiment CE-Shape-1,http://www.cis.temple.edu/-latecki/TestData/mpeg7shapeB.tar.gz,2002.

利用FFTW实现快速傅里叶变换第7篇

快速傅里叶变换 (Fast Fourier transform简称FFT) 作为时域和频域转换的分析工具, 广泛应用在数字通信、音频信号处理、图像处理、谱估计、雷达信号处理、地震勘探等各领域中。在FFT的软件实现中, 比较有名的库有fftpack, mkl及FFTW等。

一、FFT及FFTW

FFT是离散傅里叶变换的快速算法, 也可用于计算离散傅里叶变换的逆变换。快速傅里叶变换有广泛的应用, 如数字信号处理、计算大整数乘法、求解偏微分方程, 等等。

对于复数序列{x (n) } (n=0, 1, 2, ..., N-1) , 则该复数序列的离散傅里叶变换为[1]:

直接按照上述变换公式的计算复杂度是O (n2) ;快速傅里叶变换 (FFT) 可以计算出与直接计算相同的结果, 但只需要计算复杂度O (nlog (n) ) 。而实现FFT的算法有AMD Core Math Library (ACML) , intel mkl, cwplib, dfftpack, kissfft, mixfft, monnier, jmffte, numutils, FFTW3等, FFTW以其在各个平台上的优秀表现, 应用相当的广泛。

FFTW是Fastest Fourier Transform in the West (西方快速Fourier变换) 的缩写。它是计算离散Fourier变换 (DFT) 的快速C程序的一个完整集合, 它可计算一维或多维、实和复数据以及任意规模的DFT。FFTW还包含对共享和分布式存储系统的并行变换。FFTW由麻省理工学院计算机科学实验室超级计算技术组开发。

二、FFTW的工作原理

FFTW在各种平台上都能实现相当高的性能, 主是由于以下两个FFTW工作特点[2]。

1. 将长度为N的序列分解成长度为N1和N2的短序列进行计算。这些短序列的计算由一系列经过高度优化的C代码片段来完成, 这些代码片段称为solvers。solvers由planner来管理与组合起来完成计算, 其中一种组合称为一个plan。

2. 对于同一个长度为N, 其分解方式有多种, 不同的组合路径对于运算的速度影响不同, 在不同的机器与平台下其表现也会有所不同, 因此需要确定哪种组合的效果最佳。FFTW通过对各种组合的评估, 选择出其中执行速度最快的一个plan。将此plan记录下来, 以后每次做同样长度的变换时, 只要使用相当的plan即可, 这样使FFTW能适应各种运行环境下的最优。

FFTW的核心工作模块主要由planner, solver和plan来组成。其中solver是解决变换的一种独立的方法;plan是由solver生成的包含计算变换的路径和所有输入输出数组的指针的一个数据结构, 可以供后面的算法使用;planner则是管理solver生成plan, 并通过一定的算法找出执行速度最快的plan, 并将其返回给用户。具体的执行情况如下:

首先planner对一系列的solver进行初始化, 然后对给出一个指定长度变换, planner依次调用每个solver来生成plan, 当此solver可以解决的话就返回一个plan指针, 不能解决时则生成一个空指针。Planner对生成的各个plan进行执行并测量, 选择执行速度最快的plan返回给用户。而用户则调用返回的plan来执行傅里叶变换。可以将此包含执行路径及组合算法的plan保存下来, 重复使用。

由于FFTW的planner是一些独立的解决特定问题的plan的模块。故同样的planner可以用于复数的DFT, 实数的DFT及其它的变换。

三、FFTW的使用方法

在使用FFTW之前应该配置好其开发环境, 本文以VC6.0为例来讲述其环境的配置方法。其步骤如下:

1.下载FFTW for windwos版本 (ftp://ftp.fftw.org/pub/fftw/fftw-3.2-dll.zip) 。

2. 解压fftw-3.2-dll.zip到一指定目录, 如c:fftw3。在此目录下有编译所需要的动态链接库, 库的定义文件等, 但缺少静态链接的lib文件, 可利用def文件和VC的命令来生成。

3. 使用VC下的命令lib来生成VC静态链接时需要的lib文件。如下命令将生成libfftw3-3.lib文件。

4. 在VC的包含文件路径中加入fftw3头文件的路径, 如c:fftw3。

5. 在新建的工程中添加上生成的libfftw3-3.lib。

6. 然后在工程就可以使用FFTW库。

FFTW的使用过程要使用到的数据结构为fftw_complex。它的定义如下:

typedef double fftw_complex;

其是一个复数, 其数据的存储格式为:部分存储实数部分, 部分存储虚数部分。此数据结构由fftw_malloc和fftw_free函数分配与释放。其函数原型为:

在使用过程中使用到的函数有fftw_plan_dft_* (*号代表各种方式的变换) 和fftw_execute函数。其中fftw_plan_dft_*可实现生成一维及多维的及时变换和复杂变换的plan。fftw_execute则根据生成的plan来执行变换。其函数原型为:

void fftw_execute (const fftw_plan plan) ;

FFTW函数的使用步骤如下:

首先, 使用fttw_malloc分配变换中要使用的输入输出数组。

第二步, 利用它生成一个plan, 用来包含FFTW计算FFT是要包含的所有的数据。

当输入输出数组使用同一个数组时, 将执行原位运算。

第三步, 利用先前生成的plan作为参数, 调用fftw_execute执行傅里叶变换。当需要多次执行时, 可以重复调用fftw_execute来完成任务。

最后, 当不使用时, 可以调用fftw_destroy_plan, fftw_free函数将相应的plan和输入输出数组销毁。

具体的使用实例如下所示:

四、FFTW的性能评测

在2.2 GHz 32位Dual Core AMD Opteron Processor 275处理器上, 分别对AMD Core Math Library (ACML) , intel mkl, cw-plib, dfftpack, kissfft, mixfft, monnier, jmffte, numutils, FFTW3对双精度一维的不同长度的DFT进行测试。其测试的速度如图1所示[3]。

从上图可以看出其性能明显要优于其它的FFT算法的执行速度。在实践中使用FFTW将会显著缩短算法的执行时间, 提高整个系统的实时性能。

本人在做地震实时相关处理项目中使用FFTW做相关运算的处理, 取得了显著的性能提升。由此可见在实践中使用FFTW实现FFT是切实可行的, 并且能够高效地执行。

参考文献

[1]刘益成, 孙祥娥.数字信号处理.北京:电子工业出版社, 2004.

[2]王刚.遥感图像快速傅立叶变换新算法的研究与实现.中科院遥感应用研究所, 2003.

快速傅里叶变换算法第8篇

通过监测胎心率的变化来评价胎儿在宫内的安危是胎心监护中一种常用的胎儿监护方法。在围产阶段对母体内的胎儿进行监护, 不仅可以了解胎儿的发育状况, 还可以减少因缺氧、缺血而导致的胎儿出生缺陷, 甚至胎儿死亡的情况。因此, 通过对胎儿心率的变化情况进行监测对提高生育质量有着重要的意义。早至19世纪初就有产科医师通过胎心听诊估计胎儿在宫内的状况。随着超声多普勒技术的发展, 产时电子胎心监护 (electronicfetal monitoring, EFM) 已成为目前最常用的胎儿监护方法。而超声多普勒测量方法是目前最常用的胎儿心率测量方法[1]。

现阶段, 对超声多普勒原理获得的胎心音信号瞬时心率的检测一般都采用自相关算法。自相关算法提取胎儿心率, 是在时域上通过计算自相关函数2个峰值之间的时间获得胎心率值[2,3]。当胎儿心音信号受到外界噪声的干扰, 幅值变得微弱时, 很难在时域上定位出相应的峰值。因此, 噪声对瞬时心率的结果有很大的干扰性。

胎心音信号作为一种时变、非平稳且由多种成分组成的复杂信号, 传统的分析方法 (如傅里叶变换) 难以刻画不同时刻的频率成分, 无法对其进行全面的分析[4]。时频分析作为非平稳信号分析的有力工具, 该方法的主要思想是将一维时域信号变换到二维的时频平面, 通过时域和频域的联合分布信息, 从而清楚描述信号频率随时间变换的关系。时频分析方法与传统方法相比, 在非平稳信号的处理中具有突出的优越性[5]。短时傅里叶变换 (STFT) 是最早和最常用的一种用于非平稳信号的时频分析方法, 相对于其他算法, 短时傅里叶变换具有计算简单、能实时计算瞬时心率的优势。鉴于此, 本文提出了一种基于时频域的短时傅里叶变换和归一化互相关匹配的胎儿瞬时心率检测算法。首先, 将一维时域的胎心音信号通过短时傅里叶变换映射到二维时频平面;其次, 在二维时频平面上利用胎心音信号中S1音和S2音的先验信息提取特征模板;最后, 利用归一化互相关匹配算法计算瞬时心率。实验表明, 提出的算法抗干扰能力强、准确率高、易实时实现。

1 短时傅里叶变换简介

短时傅里叶变换的基本思路是把信号分成许多小时间间隔, 再用傅里叶变换分析每个时间间隔, 这样就可得到信号的一组局部频谱。STFT定义为[6]:

式中, h (t) 为时间中心和频带中心均为零的分析窗。信号x (u) 与短时窗信号h (u-t) 相乘可以有效地抑制分析时刻u=t的领域外的信号, 所以STFT是信号x (u) 在时刻t的邻域内的局部频谱, x (u) 的频谱特性则可看作是所有局部频谱的加权形式。在用STFT分析心音信号时, 要考虑STFT在时间分辨率和频率分辨率上的折中性。

2 归一化互相关匹配原理

图像相关匹配是目标跟踪[7]、物体识别、景象匹配等领域中的基本方法。归一化互相关算法是经典相关匹配方法中鲁棒性相对较高的一种方法。在将一维时域的胎心音信号映射到二维时频平面后, 本文利用归一化互相关匹配算法计算瞬时心率。文中利用心音信号中S1音和S2音的先验信息从时频图中提取特征生成二维模板。

2.1 心音信号中S1音和S2音的先验信息

正常心脏有4个心音[8]:第一 (S1) 、第二 (S2) 、第三 (S3) 、第四 (S4) 心音, 多数情况下只能听到第一和第二心音。从出现第一心音到出现第二心音的期间定义为心室收缩期;从出现第二心音到下一心动周期中出现第一心音的期间定义为心室舒张期。在一个心动周期中, 心音的主要成分为:第一心音、收缩期、第二心音、舒张期, 通过它们可以完整地表述出心音的时间特征[7]。对于一般胎儿来说, 正常情况下心脏的收缩期短于舒张期。正常胎儿心率为120~160次/min, 心动周期约为0.5 s, 其中收缩期大约占0.2 s, 舒张期大约占0.3 s。也就是说, 一般胎心音信号中, S1音到S2音在时域上的间隔大约为0.2 s。

因此, 提取特征模板时用到的S1音和S2音的先验信息有2点: (1) S1音到S2音在时域上的间隔大约为0.2 s, 特征模板在时域上的宽度应大于或等于0.2 s; (2) 整个心动周期约为0.5 s, 特征模板在时域上的宽度应小于0.5 s。

2.2 归一化互相关匹配算法

文中利用归一化互相关算法计算二维时频平面中特征模板M与二维时频图S的相关系数, 并绘制互相关曲线。归一化互相关系数定义[10]为:

式中, c (M) 为模板图像的均值, c (Sxy) 为在以 (x, y) 为起点、大小为UV区域的均值, 模板尺寸为UV。

将模板M和所匹配的图像区域按行展开成一维矢量的形式, 各自去均值后按行展开的矢量分别记为:模板矢量HM={HM (i) , i=0, , U-1}, 待匹配区域矢量HSxy={HSxy (i) , i=0, , U-1}。其中, 模板行矢量为HM (i) ={hm (i, j) -c (M) , j=0, , V-1}, 待匹配图像区域行矢量为HSxy (i) ={h Sxy (i, j) -c (Sxy) , j=0, V-1}, 那么, 式 (2) 可以改写成向量内积的形式:

式中, Re (HM (i) 、HSxy (i) ) 为HM (i) 和HSxy (i) 的内积, ||HM||2和||HSxy||2分别为HM和HSxy的2-范数。

同样, 将模板图像和待匹配区域按列展开, 也将得到类似的表达形式。式 (2) 中的归一化互相关运算是由按行 (列) 展开的各个行 (列) 矢量之间的相关运算组成的。如果2幅图像的行 (列) 与行 (列) 之间有较好的相关性, 那么这2幅图像之间也将具有较好的相关性[9]。

2.3 瞬时心率的计算

互相关曲线的峰值是二维时频图中与特征模板相似的片段和特征模板相乘得到的结果。在本算法中只要计算出互相关曲线中2个波峰之间的时间间隔, 就能求出该时刻胎心音的瞬时心率。它们之间的关系为:

3 仿真实验及结果分析

基于短时傅里叶变换的胎心音瞬时心率检测算法在MATLAB7.11.0上编程实现, 在Windows XP环境中运行。实验中用到的信号是由超声多普勒胎心仪采集得到的多普勒胎心音信号。心音的采样频率为11 025 Hz, 采用14位的A/D转换, 通过串口通信将采样数据传入上位机。

实验原始胎心音数据如图1所示, 图2是二维时频图, 图3是二维时频图对应的等高曲线。从图1和图2中可以看出, 一维时域信号不能表征心音信号的全部特征;二维时频图在时域和频域上很好地表征了信号, 并且在200~400 Hz的特征频段处特征的分辨率很高 (通过大量的数据仿真结果也证明了这一点) 。因此, 可以通过截取出200~400 Hz的特征频段代替时频图作为目标图像, 并从目标图像中利用S1和S2的先验信息提取特征模板 (也可利用先验信息提前生成一个模板) 。

图4是由特征模板与目标图像匹配出来的归一化互相关曲线, 曲线的峰值用“*”标记。从曲线中可以看出, 互相关系数在0.8以上, 表明归一化互相关算法的可靠性高。图5 (a) 是多普勒胎心音信号, 图5 (b) 和图5 (c) 分别是基于时频分析的胎心音瞬时心率检测算法和自相关算法计算出来的瞬时心率图。从图5 (b) 可知, 自相关算法提取胎儿心率是在时域上通过计算自相关函数2个峰值之间的时间获得胎心率值, 由于胎儿心音信号易受外界噪声的干扰, 幅值波动很大。当胎心音信号的幅值突然变小时, 很难在时域上定位相应的峰值, 因此, 瞬时心率的结果对噪声有很大的依赖性。采用基于时频两域联合分析的方法可以降低瞬时心率的结果对噪声的依赖性。如图5 (a) 所示 (箭头所指的数据段) , 胎心信号在第21~23 s和第33~36 s时幅值非常小, 接近于0, 采用自相关算法难以定位, 因此, 瞬时心率会出现异常。如图5 (b) 所示, 在21~23 s和33~36 s的瞬时心率值小于80。由图5 (c) 可见, 基于时频分析的胎儿瞬时心率检测算法克服了胎儿心音信号幅值变小时难以定位的弱点, 提高了胎心率计算的准确性, 特别是当胎心音信号干扰严重或胎心率超出正常范围时, 此算法更具有明显的优势。

以基于短时傅里叶变换的胎心音瞬时心率检测算法为核心的胎儿监护仪在经过4个月的临床实验中, 分别对不同孕期的90名孕妇进行了临床胎心监护, 通过与自相关算法相比较, 监护统计结果见表1。

例数

从表1中可以看出, 基于短时傅里叶变换的胎心音瞬时心率算法计算胎儿瞬时心率, 克服了胎心音信号的幅值突然变小时难以在时域上定位相应的峰值, 造成胎心率异常的现象, 比传统自相关算法的准确率高5%。

4 结论

本文提出了一种基于短时傅里叶变换的胎心音瞬时心率检测算法。通过将一维非平稳的胎心音信号变换到二维时频平面, 再在二维时频平面上利用归一化互相关匹配计算瞬时心率。该算法克服了胎心音易受噪声干扰的弱点, 实现了胎心音瞬时心率的计算, 而且还提高了胎心率计算的准确性, 为胎心音波形分析和胎心率监护提供了一种新的方法。

摘要：目的:针对母亲的活动和胎儿在子宫内的活动都会导致胎心音信号的幅度出现较大波动, 传统自相关算法难以准确测量胎儿瞬时心率的问题, 提出一种基于短时傅里叶变换的胎心音瞬时心率检测新方法。方法:对胎心音进行短时傅里叶变换, 从二维时频图中提取出胎心音特征模板并与目标图像进行归一化互相关匹配, 绘制出互相关曲线, 根据互相关曲线中的相邻峰值点时间间隔, 计算出胎儿的瞬时心率。结果:实验结果表明, 该算法的识别准确率比传统算法高5%。结论:该算法克服了胎心音易受噪声干扰的弱点, 实现了胎心音瞬时心率的计算, 而且还提高了胎心率计算的准确性, 为胎心音波形分析和胎心率监护提供了一种新的方法。

关键词：短时傅里叶变换,胎心音信号,瞬时心率,二维时频图,归一化互相关匹配

参考文献

[1]张翠英, 张尚军, 高上凯.一种用于测量胎儿心率的改进自相关算法[J].生物医学工程学杂志, 2001, 18 (3) :434-437.

[2]赵继印, 刘海英, 马洪顺, 等.基于coin5小波的多普勒胎心音信号提取算法的研究[J].中国生物医学工程学报, 2006, 25 (5) :538-541.

[3]杨晓峰, 张欣, 王金浦, 等.基于小波变换的多普勒胎儿心率检测研究[J].西安交通大学学报, 2007, 41 (8) :917-921.

[4]宋智, 李焱淼.时频分析在心电信号分析中的应用[J].大连交通大学学报, 2010, 31 (2) :56-59.

[5]王衍文, 王海滨, 程敬之.心音信号的时频分析[J].生物医学工程学杂志, 1998, 15 (1) :41-46.

[6]周静, 杨永明.心音信号的时-频分析[J].重庆大学学报, 2004, 27 (5) :159-162.

[7]Liu Gang, Su Xiu-qin, Hu Xiao-dong, et al.Real-time image tracking based on adaptive threshold in high speed TV[J].Acta Photonica Sinica, 2005, 34 (8) :1 262-1 265.

[8]韦哲.心音时-频分析方法及虚拟动态心音分析系统的研究[D].兰州:兰州理工大学, 2011.

快速傅里叶变换算法第9篇

傅里叶变换在地球物理数据降噪中应用较多, 主要原因是因为磁法、电法和地震等物探数据连续性好于地球化学数据[4,5]。而傅里叶变换在地球化学数据降噪处理上的应用尚处于尝试研究的阶段, 在文献[6,7]中,傅里叶变换用于了土壤乙烷、水系沉积物的降噪处理,但还没有出现岩石土壤地球化学降噪处理方面的应用。

1快速傅里叶变换

根据地球化学数据有序、离散的采样特点,可适用于傅里叶变换中的离散傅里叶变换( discrete fourier transform) 。土壤地球化学的每一组数据x( n) 都是由各个元素具体的含量值组成的,离散傅里叶变换表示为[8]:

在实际处理中,面对动辄成千上万的地球化学数据,计算量十分庞大,所以在计算机实现傅里叶变换的过程中为了避免繁琐的计算过程采用快速傅里叶变换算法( Fast Fourier Transformation) ,极大地提高了处理效率。计算方法如下。

从地球化学角度看,离散傅里叶变换后得到的是元素谱密度,谱密度图中高频部分代表数据中的噪声,低频部分代表背景。因此,通过对谱密度图进行分析,就可区分出噪声和背景部分,进一步对噪声进行处理。

2研究区地质背景

研究区域位于青海省大柴旦镇塔塔棱河中下游一带,区内一般海拔3 500 ~ 4 000 m。西至超力本陶勒盖,东至宝里土里盖南,面积约35 km2。行政区划隶属青海省海西州大柴旦镇塔塔棱乡所辖。位于柴北缘造山带( 南部) 与南祁连造山带( 北部) 两个构造单元的交汇处的宗务隆山晚古生代-早中生代裂陷带上( 见图1) 。柴北缘断褶带位于研究去西南侧,出露的地层以古元古界的金水口群结晶片岩系,中奥陶统以碎屑岩建造,及上石炭统变质岩系为主; 次有中-下侏罗统、白垩系、古近系、新近系等陆相碎屑岩。区域所在的宗务隆山裂陷带出露地层主要为下志留统浅海相复理式层系,次为二叠统碎屑岩、灰岩建造。区内有近东西向断裂和北西向次级断裂分布,并有规模不大的破碎蚀变带形成。测区西北部分布有塔塔棱河斑状花岗岩体。

3数据收集与预处理

研究对象是青海省大柴旦镇柴达木山南坡一带1∶ 10 000的土壤地球化学数据,共布施了34条土壤化探剖面,完成土壤地球化学测量6 km2,网度为100 m × 20 m,实际采样点位置如图2所示。

青海大柴旦镇柴达木山南坡金矿的地球化学原始数据主要分析了六种元素。通过对数据的统计和计算,得出了如下各元素地球化学特征参数 ( 表1) 。

注: 含量单位 Au 单位为 ng/g,其余为 μg/g。

4傅里叶变换降噪处理

在以往的地球化学数据处理中,降噪滤波常常以经验判断为准,导致该滤去的没有滤掉,而不该滤去的反而滤掉了,此次采用傅里叶变换后做出功率频谱图,通过对频谱图的分析就能准确找到滤波的截止幅强度。由于在未经过降噪处理的情况下,原始数据显得条理不清,尤其在Au与Cu的表现上更加复杂,几乎整个地区元素都富集。在这种异常难以选取的情况下,会使得成矿预测和找矿工作没有目的性。Au、Cu功率频谱图由图3和图4给出。

数据个数总共有3 307个,从成矿元素Au、Cu两种元素的频谱图中可以看出在Au元素中特高值为5 000以上,普遍的数据幅强度都在2 000以下, 因此将2 000作为滤波时的截止幅值,而Cu元素的特高值也同样超过了5 000,且数据大小差异不是十分明显,为避免遗失数据的真实性,以1 800作为滤波的截止强度。最后得到经过降噪后的数据,再做出波形图,与没有降噪的数据波形图对比分析如下 ( 见图5、图6) 。

通过对比降噪处理前后的波形图,不难看出原本杂乱的波形变得较为平稳并且相对光滑,原本幅值高的地方在降噪处理后仍然高,低的地方依旧很低,而且降噪处理后更加凸显高的部分,但是从观察可以看出,经过降噪后虽然特高值与各个峰值出现的位置没变,但是大小都发生了各自的变化,都相对变小,这说明降噪处理削弱了数据的能量值,但总体来看降噪效果良好。

将原始数据和降噪后的数据分别进行累积频率统计分级,按85% 、90% 、95% 共分为高中低三个等级,统计结果可表示为表2。同时,根据累积频率结果作出等值线图( 图7、图8) ,不难看出元素异常位置大都没有发生变化,但是降噪处理后的异常分布图明显异常更小更集中,尤其是在高浓度的地带。

从图7中Au元素的对比可以看出,经过降噪处理后,原本居中的大片的化学异常背景变小,最高的几个异常区域没有发生变化,而在下方和左下方出现了很小的异常区域,这说明傅里叶变换降噪处理后,噪声得到充分压制,而被噪声影响的数据能找到原本被噪声“淹没”的异常。图8中Cu元素数据经过降噪后,异常分布变得零散,在左上角,左下角有最为集中且连在一起的异常地带,尤以下部的异常最为丰富且能成片存在,而在四周有许多的异常的点集,在右上角的原本连一起的异常也变成了富集的点。其余四种元素经过降噪处理后,与Au、Cu一样,在此不再分述[9]。

5结论

经过傅里叶变换降噪处理,对比处理前后的元素异常分布图,不难看出经过降噪处理后元素异常等值图上含量值高的地带变得相对集中,元素的数据特征得到进一步体现,对异常的划分显得更加客观真实,这说明了噪声对原始数据确实有影响[10]。

地球化学数据经傅里叶变换降噪处理后,噪声得到充分压制,在被噪声影响的数据中能找到原本被“淹没”的异常。傅里叶变换充分体现了元素的空间分布特征,为成矿预测等工作扫清了障碍,澄清了地球化学场的结构。

摘要：消除地球化学数据的噪声是充分提取地球化学信息的关键。快速傅里叶变换常应用于信号处理等方面,在地球化学数据降噪中的应用还处于尝试阶段。对青海省大柴旦镇柴达木山南坡一带金矿的地球化学数据进行快速傅里叶变换,得到Au、Cu、Pb、Zn、As、Sb等六种元素的功率频谱图,进而确定地球化学区域背景和局域背景的滤波截止幅强度。使用幅强度对原始数据做滤波处理,采用低通滤波剔除数据中的噪声。对比降噪前后的元素波形图和元素异常分布图,降噪后的波形更平稳光滑,异常分布更集中,充分降低了噪声的干扰,证明快速傅里叶变换在地球化学数据降噪处理中有良好的应用前景。

快速傅里叶变换算法第10篇

互联网和数字化世界的发展,使得数字图像的安全性保护变得越来越重要。基于互联网数字产品的通信和传递十分便利,因此数字图像的安全存储和传输在许多领域变得越来越重要,如付费电视节目,医学图像系统,军事图像数据库以及秘密视频会议。为了使合法用户传递秘密的图像,而不被第三方发现,图像加密技术研究有很大的意义。优秀的图像加密保护必须完全地混乱原始图像,使其看起来完全像是随机的噪音,这样,对于非法的浏览者,图像是完全不可理解和使用。近来,许多优秀的图像安全保护的算法已经被提出。与此同时,针对提出的图像保护算法的密码分析工作也在发展,因此,从密码分析的角度来说,一些存在的算法已经变得不再安全。由于密码学和混沌系统之间的微妙关系,混沌系统作为一种加密手段广泛地用于图像安全保护。

作为一个新的时频分析工具,分数傅立叶变换(FRFT)已经引起了广泛的关注。从Mendlovic和Ozaktas在1993年完成了FRFT的光学实施开始,它已被广泛应用于光学领域[1]。2000年,Unnikrishnan G等第一次把FRFT应用于图像加密。因为FRFT的阶以及相加的性质,可以提供更多自由的图像加密方案和更大的密钥空间。在图像加密领域,FRFT已受到了越来越多的关注[2,3]。在文献[4]中,离散分数傅里叶变换(DFRFT)用于图像加密;但是,当正变换和逆变换的FRFT阶为闭合时,原始图像的部分信息可从加密图像获得。所以只用FRFT加密图像不足够安全。为了提高安全性,FRFT应与其他方法结合起来。混沌系统因其良好的性能适合数据加密,如随机性,非周期和初始条件的敏感性和控制系数[5]。混沌加密技术,该技术结合了混沌系统和其他现有的加密算法,是被视为一个有前途的加密算法[6]。因此,结合FRFT和混沌系统是一种优秀的图像加密算法,这已成为最近研究的一个热点[7,8]。文献[7]提出了一种图像加密算法基于多参数DFRFT和混沌的Logistic映射。在文献[8]中,FRFT和混沌相结合一起加密图像,并最终形成一个彩色图像用于传输。通过综合分析这些方法的优点和缺点,本文提出了一种新的基于FRFT、混沌系统和RGB映射的图像加密算法,取得了良好的加密结果。

基于混沌系统映射、分数傅里叶变换(FRFT)以及RGB映射,本文提出了新的图像安全保护算法,图像保护算法包含三个步骤:第一步是利用混沌系统产生的扰乱矩阵来打乱原始图像,第二步是通过分数傅里叶变换对图像进行加密,第三步是使用RGB映射规则转换成可视的RGB图像。提出的算法的安全性依赖于相位掩码的随机性、FRFT阶和混沌系统对初始值的敏感性。理论分析和实验结果表明,该算法在鲁棒性、安全性等方面均表现出了好的保护效果。

1 分数傅里叶变换和混沌系统的基本理论

1.1 分数傅里叶变换的基本理论

在信号处理和分析中,傅里叶变换(FT)作为时频分析工具,是最有价值和使用最频繁的工具。分数傅里叶变换是经典傅里叶变换的一般化。FRFT依赖于参数α (α=pπ/2),FRFT能被解释为信号在时域平面上做一个角度α旋转或者信号在线性调频上的分解。信号的FRFT定义如下[9]:

$X_{p} (u) = F^{p} [x (t)] = \int_{- \infty}^{+ \infty} x (t) Κ_{p} (u, t) d t (1)$

式(1)中分数傅里叶变换核函数为:

$Κ_{p} (u, t) = {\begin{cases} \sqrt{\frac{1 - j c o t α}{2 π}} \exp [j (\frac{u^{2} + t^{2}}{2} \cot α - \frac{u t}{\sin α})] α \neq n π \\ δ (t + u) α = (2 n + 1) π \\ δ (t - u) α = 2 n π \end{cases} (2)$

这里,α=pπ/2表示分数傅里叶变换的角度。当p=0时,信号经过分数傅里叶变换后,仍是信号本身。当p=1时,分数傅里叶变换就成了普通的傅里叶变换。分数傅里叶变换的逆变换就是-p次FRFT。我们用Fp表示分数傅里叶变换,这里p表示变换次数,则有加法规则:

Fp1Fp2 = Fp1 + p2

对于二维的数字图像,应当使用二维的离散分数傅里叶变换。二维的连续分数傅里叶变换其定义如下:

$\begin{array}{l} X_{p 1, p 2} (u, v) = F^{p 1, p 2} [x (s, t)] \\ = \int_{- \infty}^{+ \infty} \int_{- \infty}^{+ \infty} x (s, t) Κ_{^{p 1, p 2}} (s, t, u, v) d s d t (3) \end{array}$

这里, x(s, t)表示原始的二维图像信号。

二维分数傅里叶变换核函数为:

$\begin{array}{l} Κ_{p 1, p 2} (s, t, u, v) = \sqrt{\frac{1 - j c o t α}{2 π}} \sqrt{\frac{1 - j c o t β}{2 π}} \cdot \exp [j (\frac{s^{2} + u^{2}}{2} \cot α - \\ \frac{s u}{\sin α} \end{array}$

expj $\frac{t^{2} + v^{2}}{2} \cot β - \frac{t v}{\sin β})] (4)$

利用二维FRFT核函数可分性的特点,二维FRFT能被分解成两个一维的变换,或者分别在列方向上和行方向上的一维分数傅里叶变换。分数傅里叶变换可以表示为张量积的形式:FpX=MpX(Mp)T,Mp是p阶离散分数傅里叶变换(DFRFT)的矩阵,X是数字图像的矩阵表示[10]。分数傅里叶变换的逆变换表示为:X=M-p[FpX](M-p)T。因此,这种计算方式是一种节省时间的方法,其时间复杂度为o(nlog(n))。

1.2 混沌系统

近代科学对混沌系统的定义具有不准确性。混沌系统的数学定义是由Tien-Yien Li和James A.Yorke在1975年提出的,他们第一次给出了混沌系统准确的科学含义。Li-York混沌系统[11]定义如下:

假设存在一个连续自映射f ∶iI⊂R, I是R集合里面的一个闭区间。如果有一个不可数的集合S⊂I满足下面的三个条件:

(1) S中不含有周期点;

(2) 对任意的两个变量X1,X2∈S(X1≠X2)有:

${\begin{cases} \lim_{i \infty} \sup | f^{i} (X_{1}) - f^{i} (X_{2}) | > 0 \\ \lim_{i \infty} \inf | f^{i} (X_{1}) - f^{i} (X_{2}) | = 0 \end{cases} (5)$

这里fi=f(f(f()))。

(3) 对于任意X1∈S和周期点P∈I有:

$\lim_{i \infty} \sup | f^{i} (X_{1}) - f^{i} (Ρ) | > 0 (6)$

则称f在集合S中是混沌的,或者说f是一个混沌映射。

混沌信号可以由各种迭代方程产生。目前对于迭代方程的研究,主要集中于Logistic混沌序列方程。Logistic映射定义如下:

Xn+1 = μXn(1-Xn) (7)

这里μ (0 μ 4)表示混沌系统的控制系数, X0(0 X0 1)是系统的初始状态条件。根据Lyapunov指数可以知道,当3.5699 μ 4,该映射处于混沌状态,[3.5699,4]称为混沌区域。也就是说,混沌区域内产生的相应迭代序列是非周期的、非收敛的,且对初始条件极其的敏感[12]。Logistic映射系统表明对不同的μ值将会有不同的行为表现,且会逐渐变成倍周期分叉。

2 图像保护加密步骤

2.1 图像加密过程

流程图如图1所示,展示了本文提出的加密算法的加密过程。

让 $Ο = {o_{i, j} | i = 1, \dots, a$ ;j=1,,b}表示大小为ab的原始图像,图像加密步骤描述如下:

Step1 利用混沌系统扰乱图像的像素:通过式(7)和密钥{μ0,X00}得到混沌序列X,序列X的长度为L=ab。将X按从小到大排序成新的序列Xs。把二维的原始图像矩阵转换成一维的长度L的序列S。通过下面的式(8)利用Xs将S扰乱成S′。

S′(Xs(i))=S(i) i=1,2,,ab (8)

在经过Logistic映射后,原始图像像素的时域位置被扰乱。

Step2 在X方向进行α阶DFRFT:把Step1获得的矩阵看成一个行向量O′=[o1,o2,,on],并对这行向量实施α阶离散分数傅里叶变换,于是一个被加密的复数矩阵被获得。

Step3 利用混沌系统扰乱图像的像素:扰乱规则类似步骤1。Logistic映射后,原始图像像素的FRFT域位置被扰乱。

Step4 在Y方向进行β阶DFRFT:把Step3获得的矩阵看成一个列向量O′′=[o1,o2,,on]T,并对这列向量实施β阶离散分数傅里叶变换,于是又一个被加密的复数矩阵被获得。

Step5 RGB映射[8]:上述步骤获得的加密图像矩阵是一个复数矩阵,这是难以在实际应用中传输和存储。设经变换并置乱后的矩阵 N 的复数数据为N (s, t) =Re(s, t) +iIm(s,t)。RGB分量的映射规则为:把被加密图像矩阵的实部和虚部的映射成RGB图像的形式。形成一幅用于传输的彩色图像。映射步骤如下:

(1) 复数矩阵分成两部分:实部和虚部分别存储在彩色图像的分量R和G中。

(2) 剩下的分量B作为控制标志。分量B=(b8b7b6b5b4b3b2b1)binary,每个比特位的值bi代表一个控制字。Fr和Fi表示实部和虚部符号控制位,而Mr和Mi分别表示实部和虚部的缩放倍数,如表1所示。

(3) 如果实部为正,Fr设置为0,否则设置Fr为1;同理,对 Fi 位也作类似的处理。接着判断实部绝对值是否小于 255,若是,则把其绝对值直接存入 R 中相应位置, Mr位置为0;否则, Re= Re mod 255+ Re /255,则 $| R e | m o d 255$ 的值存入 R中相应位置,同时 $| R e | / 255$ 转为二进制存入Mr。对虚部也作类似的处理。

(4) 重复第3步直到整个加密图像数据处理完。矩阵的相同位置上R、G、B 分量组合表示了原灰度图像加密后的数据 ,即被用作传输和保存彩色的图像。

2.2 图像解密过程

解密是加密的逆过程,解密的步骤可以简要介绍如下:

Step1 对彩色图像每个像素位的R、G和B做逆映射操作,恢复复数矩阵。

Step2 对复数矩阵在Y方向实施β阶DFRFT。

Step3 利用密钥(μ1,X01)实施逆混沌乱扰。

Step4 对复数矩阵在X方向实施α阶DFRFT。

Step5 利用密钥(μ0,X00)实施逆混沌乱扰。于是,获得解密后的图像。

3 模拟实验和安全性分析

本文采用标准测试灰度图像Camerman(256256)作为仿真图像,在Matlab 7. 0平台上进行下列仿真实验。

3.1 图像的FRFT实验

对原始图像(Camerman.tif)实施不同阶的离散分数傅里叶变换,模拟结果显示,如图2所示。

从上面的仿真结果,可以清楚地看到,α值较大,图像在水平方向上有更大的变化时;当β值是较大,图像在垂直方向有很大改变。因此,本文分别在X和Y方向进行了FRFT,从而也简化了FRFT的计算。

3.2 加密和解密结果

为了验证加密和解密的效果,我们利用提出的算法与密钥K=(3.19,0.84,0.75,3.91,0.45,0.31)加密标准测试图像Camerman。各种加密图像的效果如图3所示。拥有正确密钥的解密图像如图4所示。

从图3和图4中,我们可以看到,加密图像和原图像之间有一个很大的区别,而解密图像和原图像之间的差异是很轻微的。

3.3 统计特性分析

对密码系统来说,密文的统计分析是非常重要的。一个理想的密码应该具有很强的鲁棒性才能反对任何统计攻击。为了证明本文提出的图像加密算法的安全性,进行统计分析。

如图5所示,原始图像和加密图像之间的直方图是极其的不同。因此,通过对加密图像做统计分析几乎不可能解密图像。图5(b)表明,加密图像的直方图是平坦的,这意味着本文的加密算法有一个良好的统计特性。这是因为混沌系统的随机性,对图像像素位置的扰乱达到了良好的扰乱效果。更重要的是,每个像素值可以同时在时间域和FRFT域(实施FRFT后)得到极好的扰乱,分数阶Fourier域后,不仅仅是在时域。

3.4 密钥空间和敏感度分析

本文所提出算法的密钥包括α、β(FRFT在X和Y方向上的分数阶)和X00、X01、μ0、μ1(混沌系统初始条件和控制系数),取值范围为-2α,β2,3.5699μ4,0X00、X011。根据计算机的双浮点精度,8位和15个有效数字进行了分析。假设,试图进行一次解密密文花费时间是在秒级单位,尝试完密钥空间内所有的密钥则需要大约(91010101441014)2/(3.15107)=4.111071年。从理论上讲,它是一个非常大的密钥空间,足以确保信息安全。

尽管一个大的密钥空间,良好的图像加密算法也应对钥匙敏感。如图6所示,正确的密钥解密的图像显示在图6(a),使用错误的密钥,即使有点偏移,将导致不正确的解密图像,如图6(b)-(f)所示。因此,我们可以得出结论,加密算法对钥匙是非常敏感的。

作为一个解密图像的对象方法,MSE(平均方误差),用于模拟。MSE是定义如式(9),它直接反映两个图像之间质量差异。

$Μ S E = \frac{1}{Ν Μ} \sum_{i = 1}^{Ν} \sum_{j = 1}^{Μ} [Ο (i, j) - E (i, j)]^{2} (9)$

这里,O、E分别原始图像和加密图像。正确密钥K=(3.19,0.84,0.75,3.91,0.45,0.31)下,MSE曲线如图(7)所示。当解密的密钥和正确的密钥之间的每个参数的差异是0.01或以上时,MSE的差异将相当大。这就是说,当密钥的一个参数误差不小于0.01时,解密将失败,它可以有效地保护信息。

3.5 鲁棒性分析

对闭塞攻击的鲁棒性进行了分析。正如图8所示,即使宿主图像的某些部分被剪切,整个隐藏的图像信息,仍可以重建,但也有一些额外的噪音。加密图像剪切部分越多,解密后会有更多的噪音。但毕竟,一些图像处理技术可以改善图像品质。

4 结语

基于分数傅里叶变换和混沌系统,本文提出一个新的图像加密算法。图像分别在X和Y方向都分别进行了分数傅里叶变换,扰乱都是在时域和分数阶Fourier域。且在时域和FRFT域都进行了扰乱操作,因此像素位置和像素值取得很好的扰乱效果。实验结果表明,该算法具有相当的稳定性和安全性,很容易实现,加密速度快,密钥空间大,且对密钥敏感。

摘要：针对现有图像加密算法的缺点,提出一种新的图像加密算法。该算法采用混沌系统映射、分数傅里叶变换(FRFT)以及RGB映射。整个图像加密包括三个方面:第一是利用混沌系统产生的扰乱矩阵来打乱原始图像,第二是通过分数傅里叶变换对图像进行加密,第三是使用RGB映射规则转换成可视的RGB图像。由于相位掩码的随机性、FRFT阶和混沌系统对初始值的敏感性,算法的安全性得到了极大的保证。理论分析和实验结果表明,该算法安全性好,具有良好的研究价值。

关键词：图像保护,分数傅里叶变换,混沌系统,RGB映射

参考文献

[1]David Mendlovic,Haldun M Ozaktas.Fractional Fourier Transforms and Their Optical Implementation I[J].J.Opt.Soc.AM.A,1993,10:1875-1881.

[2]Namias V.The fractional order Fourier transform and its application to quantum mechanics[J].J.Inst.Math.Appl,1980,25:241-265.

[3]Naveen Kumar Nishchal,Joby Joseph,Kehar Singh.Securing Informa-tion using Fractional Fourier Transform in Digital Holograph[J].Optics Communication,2004,235:253-259.

[4]He Junfa,Li Jun,Wang Hongxia.Digital image encryption transform implemented by asymmetric discrete fractional Fourier transform[J]Optical Technique,2005,31:410-412.

[5]L.Kocarev.Chaos-based cryptography:a brief overview[J].IEEE Cir-cuits and Systems Magazine,2001,1:6-21.

[6]Michael R Lyu.Handbook of software reliability engineering[M].New York:IEEE Computer Society Press,1996.

[7]Jun Lang,Ran Tao,Yue Wang.Image encryption based on the muti-ple-parameter discrete fractional Fourier transform and chaos function[J].Optic Communications,2010,283:2092-2096.

[8]杨倬,冯久超,方勇.一种基于混沌和分数阶傅里叶变换的图像加密算法[J].计算机科学,2008,35:239-240.

[9]Tao Ran,Deng Bing,Wang Yue.Fractional Fourier Transform and Its Applications[M].Beijing:Tsinghua University Press,2009.

[10]Pei S C,Tseng C C,Yeh M H,et al.Discrete Fractional Hartley and Fourier Transforms[J].IEEE Trans.Circuits Syst.II,June1998,45:665-675.

[11]Tien Yien Li,James A Yorke.Period Three Implies Chaos[J].The A-merican Mathematical Monthly,December1975,82:985-992.

快速傅里叶变换算法

快速傅里叶变换算法第1篇

快速傅里叶变换算法第2篇

快速傅里叶变换算法第3篇

快速傅里叶变换算法第4篇

分数阶傅里叶变换的数值实现第5篇

基于快速傅里叶变换的形状检索第6篇

利用FFTW实现快速傅里叶变换第7篇

快速傅里叶变换算法第8篇

快速傅里叶变换算法第9篇

快速傅里叶变换算法第10篇

快速傅里叶变换算法

猜你喜欢

栏目

最近更新

热门标签

快速傅里叶变换算法

快速傅里叶变换算法 第1篇

快速傅里叶变换算法 第2篇

快速傅里叶变换算法 第3篇

快速傅里叶变换算法 第4篇

分数阶傅里叶变换的数值实现 第5篇

基于快速傅里叶变换的形状检索 第6篇

利用FFTW实现快速傅里叶变换 第7篇

快速傅里叶变换算法 第8篇

快速傅里叶变换算法 第9篇

快速傅里叶变换算法 第10篇

快速傅里叶变换算法

猜你喜欢

栏目

最近更新

热门标签

快速傅里叶变换算法第1篇

快速傅里叶变换算法第2篇

快速傅里叶变换算法第3篇

快速傅里叶变换算法第4篇

分数阶傅里叶变换的数值实现第5篇

基于快速傅里叶变换的形状检索第6篇

利用FFTW实现快速傅里叶变换第7篇

快速傅里叶变换算法第8篇

快速傅里叶变换算法第9篇

快速傅里叶变换算法第10篇