北大医学数字图像处理2.1傅里叶变换(Fourier-Transform)

第2章图像变换我们将在第二、三章讨论频率域处理法，也就是变换域处理法。从图像工程角度来看，变换域处理是指对数字图像进行预处理，图像处理中的图像增强、图像恢复、图像编码压缩、图像分析描述等都是以图像变换为基础的。基本概念是：把图像信号从空间域变换到频率域，从而从另一个角度分析图像信号的特性。图像频域处理具有高运算速度的特点，同时，可以利用已有的二维数字滤波技术进行多种图像处理，得到了广泛的应用。频域处理的关键是要首先将图像从空间域变换到频域，再将在频率上进行了不同处理的结果反变换到空间域。由于图像系统可以看作一个线性系统，通常我们进行的线性正交变换运算。线性变换：其基本线性运算式是严格可逆的，且满足一定的正交条件，即酉变换（unitarytransform）。考察多维空间坐标的基轴方向互相正交[1]？设有一组基向量集合11121212221212,,,nnnnnaaaaaaaaa⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦nnaaa彼此正交，该组基向量组成正交矩阵1111212122212nnnnnnaaaaaaAaaa⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦满足TTAAAAI==。其中为A的转置矩阵，TA1112112121221222TTaaaaAaaaa⎡⎤⎡==⎤⎢⎥⎢⎥⎣⎦⎣⎦。首先复习了酉矩阵、酉变换和正交变换，图像处理的预处理就是酉变换。●一维的正交变换矩阵A对任一数据向量f进行运算gAf=，由1fAg−=AfAgτ⎯⎯⎯⎯⎯→=若为正交矩阵恢复f。●酉矩阵（UnitaryMatrix，么正矩阵，一元矩阵）若A的每一个元素可能为复数，则正交的条件可用酉矩阵满足者为1⎯⎯⎯→⎯=∗−TAA●酉变换用酉矩阵A对任意向量f的运算和恢复，即若A为酉矩阵，2则满足下式101**0()(,)()01()(,)()01NnNTkgAfgkaknfnkNfAgfnakngknN−=−===≤≤==≤≤∑∑−−如前所述，若A不是复数时，就为正交变换，TAA=−1，故各种正交变换都是酉变换。3本章内容傅里叶变换余弦变换沃尔什变换哈达玛变换K-L变换与图像压缩42.1傅里叶变换(FourierTransform)Fourier’book,LaTheorieAnalytiquedelaChaleur(TheAnalyticTheoryofHeat),is“agreatmathematicalpoem”.-----JamesClarkMaxwellFourier’sideabuiltintothecommonsenseofoursociety.------T.W.Korner不规则数据和复杂的曲线中包含了信号的全部信息，这些信息隐藏的地方是我们的智力所达不到的。而傅立叶分析却把这些信息翻译成了简单明了的形式。-----YvesMeyerIn1807Fouriershowedthatanyfunction(evenirregularone)canbeexpressedasthesumofaseriesofsinesand/orcosines…数学：傅里叶变换是最早(1807)研究和应用最广泛的酉变换。他在进行热的研究和解微分方程的努力中，在牛顿之后的近150年，Fourier为解决所有类型的数字式方程式和线性偏微分方程提供了一个实用的方法。其思想在数学领域里主宰了近100多年，并发展了大量的分支，甚至包括数论和概率论学科。傅里叶变换是一种刻画函数空间，求解微分方程，进行数值运算的有效工具。它可以把许多常见的微分、积分和卷积运算化简成代数5运算。工程技术：Fourier思想被广泛用于线性程序编程，结晶学，电话及无线电的设备及医院的X光机上,IP长途电话和遥远星系的射线分析。量子力学:将一个()rψ按平面波()exp.Cipr作傅立叶展开()().33/21(),2priprpedpψφπ==∫k这里()pφ为波函数按平面波展开后的波幅，()2pφ表示()rψ含有平面波(exp.Cipr)的份额。其逆表达式为()().33/21()2pripredrφψπ−=∫在Fourier变换的“物理空间”，对应微观粒子的几何位置；在Fourier变换的“傅立叶空间”一边，对应的粒子的动量或粒子的波。量子力学中的“测不准原理”：在微观领域里，我们无法同时精确测量微观粒子的位置和动量，或者说，无法同时测定微观粒子的能量和时间。这在的傅立叶分析中就变成自然结果，因为我们不可能同时精确地同时知道傅立叶变换的两个方面。随着量子力学的发现，Fourier分析成了自然本身的语言。所以，正如一位美国科学家所说的，长期以来，研究应用FT的科学家不厌其烦地将正余弦相加以建立信号。他们“读”Fourier系数以得到他们想要的信息，就像一些音乐家通过读乐谱而能静静地听到音乐一样。他们会花很长时间愉悦地在傅立叶空间里畅游而不进入物理空间。62.1.1One-DFourierFunction设为实变量x的连续、可积函数，则其傅氏变换为)(xf[]∫+∞∞−−=dxuxjxfuFπ2exp)()(=()()()()juRujIuFueΦ+=，其中1j=−；u为频域变量,x为空域变量;为实数部分,为虚数部分;)(uR)(uI()Fu为幅度函数,)(uΦ为相角。)(xf的傅里叶谱：[]2122)()()(uIuRuF+=，相位：)()()(1uRuItgu−=Φ能量谱：)()()()(222uIuRuFuE+==傅氏变换中，中包含了无限多项的正、余弦之和()(uF[]uxjuxuxjπππ2sin2cos2exp−=−),这里每一个u值对应正弦与余弦的频域。定义傅立叶逆变换，[]∫+∞∞−=dxuxjuFxfπ2exp)()(傅里叶变换对：)()(uFxf⇔72.1.2Two-D连续傅里叶变换（ContinuousFourierTransforminTwoDimension）对2－D的连续可积函数(,)fxy，其傅氏变换对(FourierTransformPair)[][]),(),()(2exp),(),()(2exp),(),(vuFyxfdudvvyuxjvuFyxfdxdyvyuxjyxfvuF⇔⇒+=+−=∫∫∫∫∞+∞−+∞∞−ππ傅里叶谱：[]2122),(),(),(vuIvuRvuF+=，相位：),(),(),(1vuRvuItgvu−=Φ傅里叶能量谱：),(),(),(22vuIvuRvuE+=82.1.3One-D离散傅里叶变换（DFT,DiscreteFourierTransform）●傅里叶变换对(TransformPairfortheFouriertransform)的离散化采样：用N个互相间隔xΔ单位采样，使其成为系列)(xf{}0000(),(),(2),,([1])fxfxxfxxfxNx+Δ+Δ+−Δ规定)()(0xxxfxfΔ+=，这时x为离散值1,,2,1,0−=Nx，上述序列表示为{})1(,),2(),1(),0(−Nffff这时傅里叶变换对为[][]∑∑−=−==−=1010/2exp)()(/2exp)(1)(NuNxNuxjuFxfNuxjxfNuFππ式中1,,2,1,0)()(−=Δ=NuuuFuF，相应于uNuuΔ−ΔΔ)1.(,2,,0且)(1xNuΔ=ΔTransformPair:)()(uFxf⇔离散傅里叶变换满足正交条件[][]1121201if1exp2/exp2/0elsewhereNxuujuxNjuxNNππ−==⎧−=⎨⎩∑实用上，常令次幂。mN2=9●矩阵表示法例、的原信号序列4=N{})3(),2(),1(),0()(ffffxf=的FT[][]∑∑−=−=−=−=140104/2exp)(41/2exp)(1)(xNxuxjxfNuxjxfNuFππ展开：[]⎥⎦⎤⎢⎣⎡+++==⎥⎦⎤⎢⎣⎡+++==⎥⎦⎤⎢⎣⎡+++==+++==−−−−−−−−−292623026242202322200000)3()2()1()0(41)3(,3)3()2()1()0(41)2(,2)3()2()1()0(41)1(,1)3()2()1()0(41)0(,0πππππππππjjjjjjjjjefefefefFuefefefefFuefefefefFuefefefefFu把上面的e指数项可写成矩阵形式，引入了DFT的矩阵表（N＝4）：000023022224602223690222(0)(1)(2))(3jjjjjjjjjeeeefeeeeffeeeefeeeeπππππππππ−−−−−−−−−⎡⎤⎢⎥00003022020(014)(1)(2)(3)14jjjjjjFFFFeeeeeeeeeeeeeeππππππ−−−−−−⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎣⎦⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦=322(0)(1)(2)(3)jjjffffeeπππ−−−⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎣⎦记为,10FWf=其W为变换核矩阵。由此，傅里叶变换就表示为向量、矩阵的简式。为了方便矩阵计算，重写DFT表式为[]10101()()exp2/1()NxNuxNxFufxjuxNfxWNπ−=−==−=∑∑N也就是说，上述矩阵中的每一个矩阵元表示成2juxuxNNWeπ−=。对前述N＝4情况，000102034444101112134444202122234444303132334444111111111111⎡⎤⎡⎤⎢⎥⎢⎥−−⎢⎥⎢⎥==⎢⎥⎢⎥−−⎢⎥⎢⎥−−⎢⎥⎣⎦⎣⎦。112.1.4快速傅里叶变换FFT(FastFourierTransform)问题的提出，当时有快速算法。nN2=对F＝Wf或1100121()()exp()NNuxNxxuxFufxjfxNNNπ−−==⎡⎤=−=⎢⎥⎣⎦∑∑W，包括了N次乘法和（N－1）次复数加法。又考虑有N个值，故全部计算需次复数乘法和)(uF)(uF2N)1(−×NN次加法。1965年CooleyandTukey提出的算法结果是乘法次数：NNN22log2→加法次数：NNNN2log)1(→−N很大时，效果显著，如NDFT计算量FFT计算量DFTFFT计算量计算量20484，194，30422，528186．2人们据此设计FFT硬件单元，速度比软件快了数百倍。软件算法有不断的新算法。下面介绍基于的算法。nN2=FFT算法包括两大类：时域分组：将2juxuxNNWeπ−=中的x不断分解为奇偶表达式频域分组：将2juxuxNNWeπ−=中的u不断分解为奇偶表达式●时域分组（按时间抽取的FFT算法）按的x＝偶数部分和奇数部分)(xf)2(xf(21),0,1,2,,12Nfxx+=−12分成两组，的DFT可用两个)(xf2N采样点的DFT计算：1011222(001F()()122(2)(21)2NuxNxNNuxuxNNxxufxWNfxWfxWNN−=−−+===⎡⎤⎢⎥=++⎢⎥⎣⎦∑∑∑21)2222/2/21122/2/200()()122(2)(21)2eojuxjuxuxuxNNNNNNuxuxuNNxxFuFuWeeWfxWfxWWNNππ−−−−==⎛⎞===⎜⎟⎝⎠⎡⎤⎢⎥=++⎢⎥⎢⎥⎢⎥⎣⎦∑∑∵N1()()()0,1,2,,/212ueNoFuFuWFuuN⎡⎤⇒=+=⎣⎦−（1）这里，和分别对应N个抽取值序列中的偶数点序列和奇数点序列，而)(uFe)(uFo)2()2(21)2(2⎥⎦⎤⎢⎣⎡+++=++NuFWNuFNuFoNuNeu以N/2为周期（从和的原式可见）)(uFe)(uFo(u))2((u))2(ooeeFNuFFNuF=+=+Becauseof222expexp()2NNuuuuNNNNNNππ+⎡⎤==−=−=−⎢⎥⎣⎦uN，then131()()()22ueNoNFuFuWFu⎡⇒+=−⎣⎤⎦（2）只要分别计算出0到