数字图像频域变换.ppt

数字图象处理,北京大学计算机研究所 陈晓鸥,第三节 频域变换,傅立叶变换导言理论基础、连续与离散的傅立叶变换二维傅立叶变换特性可分离性、周期与共轭对称、平移性、旋转特性、线性与相似性、均值性、拉普拉斯、卷积与相关快速傅立叶变换FFT算法、逆向FFT算法、算法实现,第三节 频域变换,2.2.1 傅立叶变换导言理论基础连续与离散的傅立叶变换,第三节 频域变换:理论基础,理论基础线性系统卷积与相关,第三节 频域变换:理论基础,线性系统系统的定义：接受一个输入，并产生相应输出的任何实体。系统的输入是一个或两个变量的函数，输出是相同变量的另一个函数。,系统,x(t)输入,y(t)输出,第三节 频域变换:理论基础,线性系统线性系统的定义：对于某特定系统，有：x1(t)y1(t)x2(t)y2(t)该系统是线性的当且仅当：x1(t)+x2(t)y1(t)+y2(t)从而有：a*x1(t)a*y1(t),第三节 频域变换:理论基础,线性系统线性系统移不变性的定义：对于某线性系统，有：x(t)y(t)当输入信号沿时间轴平移T，有：x(t-T)y(t-T)则称该线性系统具有移不变性,第三节 频域变换:理论基础,卷积卷积的定义离散一维卷积二维卷积的定义离散二维卷积相关的定义,第三节 频域变换:理论基础,卷积的定义对于一个线性系统的输入f(t)和输出h(t)，如果有一个一般表达式，来说明他们的关系，对线性系统的分析，将大有帮助卷积积分就是这样的一般表达式 h(t)=g(t-)f()d 记为：h=g*f-g(t)称为冲激响应函数,第三节 频域变换:理论基础,离散一维卷积 h(i)=f(i)*g(i)=f(j)g(i-j)j二维卷积的定义 h(x,y)=f*g=f(u,v)g(x u,y v)dudv-,第三节 频域变换:理论基础,离散二维卷积h(x,y)=f*g=f(m,n)g(x m,y n)m n相关的定义 h(t)=g(t+)f()d 记为：y=g x-,第三节 频域变换,连续与离散的傅立叶变换一维连续傅立叶变换二维连续傅立叶变换离散傅立叶变换离散傅立叶变换的计算与显示,第三节 频域变换:傅立叶变换,一维连续傅立叶变换：定义设 f(x)为实变量x的连续函数，f(x)的傅立叶变换表示为Ff(x),定义为：Ff(x)=F(u)=f(x)exp(-j2ux)dx 其中 j2=-1-如果给定F(u),f(x)可以由傅立叶逆变换得到：FF(u)=f(x)=F(u)exp(j2ux)du,第三节 频域变换:傅立叶变换,一维连续傅立叶变换：几个概念 假设函数f(x)为实函数。但一个实函数的傅立叶变换可能为复函数：F(u)=R(u)+jI(u)（1）f(x)的傅立叶模记为：|F(u)|F(u)|=R2(u)+I2(u)1/2（2）f(x)的傅立叶模平方记为：P(u)P(u)=|F(u)|2=R2(u)+I2(u),第三节 频域变换:傅立叶变换,一维连续傅立叶变换：几个概念（3）f(x)的傅立叶相位记为：(u)(u)=tan-1(I(u)/R(u)（4）傅立叶变换中的变量u通常称为频率变量 这个名称源于尤拉公式中的指数项 exp-j2ux=cos2ux-jsin2ux 如果把傅立叶变换的积分解释为离散项的和，则易推出F(u)是一组sin和cos函数项的无限和，其中u的每个值决定了其相应cos,sin函数对的频率。,第三节 频域变换:傅立叶变换,二维连续傅立叶变换 如果f(x,y)连续可积，并且F(u,v)可积，则存在以下傅立叶变换对，其中u,v为频率变量：Ff(x,y)=F(u,v)=f(x,y)exp-j2(ux+vy)dxdy-FF(u,v)=f(x,y)=F(u,v)expj2(ux+vy)dudv-,第三节 频域变换:傅立叶变换,二维连续傅立叶变换 二维傅立叶模、相位和模平方分别为：模：|F(u,v)|=R2(u,v)+I2(u,v)1/2 相位：(u,v)=tan-1(I(u,v)/R(u,v)模平方：P(u,v)=|F(u,v)|2=R2(u,v)+I2(u,v),第三节 频域变换:傅立叶变换,离散傅立叶变换 假设连续函数f(x),通过取N个x单位的采样点，被离散化为一个序列：f(x0),f(x0+x),f(x0+2x),f(x0+N1 x)无论将x作为离散的或连续的变量，在子序列中来研究都将是方便的，仅仅依赖于讨论的上下文。为作到此要求定义：f(x)=f(x0+xx),第三节 频域变换:傅立叶变换,离散傅立叶变换 其中假设x现在的离散值是：0,1,2,N-1。f(x0),f(x0+x),f(x0+2x),.,f(x0+N1x)表示相对与连续函数的任意N个统一的空间采样,第三节 频域变换:傅立叶变换,离散傅立叶变换函数f(x0+xx)的离散傅立叶变换对有：N-1F(u)=1/N f(x)exp-j2ux/N x=0u=0,1,2,.N-1 N-1 f(x)=F(u)expj2ux/N u=0 x=0,1,2,.N-1,第三节 频域变换:傅立叶变换,离散傅立叶变换:二维 M-1 N-1F(u,v)=1/MNf(x,y)exp-j2(ux/M+vy/N)x=0 y=0u=0,1,2,M-1;v=0,1,2,.N-1 M-1 N-1 f(x,y)=F(u,v)expj2(ux/M+vy/N)u=0 v=0 x=0,1,2,.N-1;y=0,1,2,.N-1,第三节 频域变换:傅立叶变换,离散傅立叶变换的计算与显示离散傅立叶变换的计算举例离散傅立叶变换的显示,第三节 频域变换:傅立叶变换,离散傅立叶变换的计算举例,x,f(x0)=f(x0+x),0,1,2,3,1,2,3,4,第三节 频域变换:傅立叶变换,F(0)=1/4f(x)exp0=1/4f(0)+f1(1)+f(2)+f(3)=1/4(2+3+4+4)=3.25F(1)=1/4f(x)exp-j2x/4)=1/4(2e0+3e j2/4+4e j22/4+4e j23/4)=1/4(-2+j)F(2)=-1/4(1+j0)F(3)=-1/4(2+j),第三节 频域变换:傅立叶变换,离散傅立叶变换的计算举例因为，函数f(x,y)的傅立叶变换是f(x,y)积分的函数，因此计算每一个傅立叶变换值，原函数f(x,y)的每一个点都需要参予,第三节 频域变换:傅立叶变换,离散傅立叶变换的显示 通过对傅立叶变换模，来显示傅立叶变换图象。由于模的值域大于显示的值域，因此要进行动态值域的压缩D(u,v)=c log(1+|F(u,v)|)其中：c=255/k;k=max(log(1+|F(u,v)|)值域0,k的上限（最大值）,第三节 频域变换:傅立叶变换,离散傅立叶变换的显示,第三节 频域变换:傅立叶变换,离散傅立叶变换的显示,第三节 频域变换:二维傅立叶变换特性,2.2.2 二维傅立叶变换特性可分离性周期与共轭对称平移性旋转特性,线性与相似性 均值性 拉普拉斯 卷积与相关,第三节 频域变换:二维傅立叶变换特性,可分离性二维离散傅立叶变换DFT可分离性的基本思想是：二维DFT可分离为两次一维DFT应用：二维快速傅立叶算法FFT，是通过计算两次一维FFT实现的,第三节 频域变换:二维傅立叶变换特性,可分离性的定义 M-1 N-1 F(u,v)=1/MN f(x,y)e(-j2vy/N)e(-j2ux/M)x=0 y=0u=0,1,2,M-1;v=0,1,2,.N-1 M-1 N-1 f(x,y)=F(u,v)e(j2vy/N)e(j2ux/M)u=0 v=0 x=0,1,2,.N-1;y=0,1,2,.N-1,第三节 频域变换:二维傅立叶变换特性,可分离性成立的推导先对行（y变量）做变换：N-1F(x,v)=1/Nf(x,y)exp(-j2vy/N)y=0然后对列（x变量）进行变换：M-1F(u,v)=1/MF(x,v)exp(-j2ux/M)x=0,第三节 频域变换:二维傅立叶变换特性,先对行做变换：,然后对列进行变换：,f(x,y),（0，0）,（N-1，M-1）,x,y,F(x,v),（0，0）,（N-1，M-1）,x,v,F(x,v),（0，0）,（N-1，M-1）,x,v,F(u,v),（0，0）,（N-1，M-1）,u,v,第三节 频域变换:二维傅立叶变换特性,平移性定理平移性的描述函数自变量的位移的傅立叶变换产生一个复系数 Ff(x-a,y-b)=exp-j2(au+bv)F(u,v),第三节 频域变换:二维傅立叶变换特性,平移性成立的证明用一维函数为例进行证明：设位移为a,f(x-a)的傅立叶变换为：Ff(x-a)=F(u)=f(x-a)exp(-j2ux)dx-将积分乘以 1=exp(-j2au)exp(j2au)且设：v=x-a,dv=dx,第三节 频域变换:二维傅立叶变换特性,平移性成立的证明Ff(x-a)=F(u)=f(x-a)exp(-j2ux)exp(-j2au)exp(j2au)dx-=exp(-j2au)f(x-a)exp(-j2ux)exp(j2ua)dx-,第三节 频域变换:二维傅立叶变换特性,=exp(-j2au)f(x-a)exp-j2u(x-a)dx-=exp(-j2au)f(v)exp-j2uvdv-=exp(-j2au)F(u),第三节 频域变换:二维傅立叶变换特性,周期与共轭对称周期性的描述：离散傅立叶变换DFT和它的逆变换是以N为周期的对于一维傅立叶变换有：F(u)=F(u+N)对于二维傅立叶变换有：F(u,v)=F(u+M,v+N),第三节 频域变换:二维傅立叶变换特性,周期与共轭对称共轭对称性的描述：傅立叶变换结果是以原点为中心的共轭对称函数对于一维傅立叶变换有：F(u)=F*(-u)对于二维傅立叶变换有：F(u,v)=F*(-u,-v),第三节 频域变换:二维傅立叶变换特性,共轭对称性证明以一维傅立叶变换为例证明：F(u)=f(x)exp-j2uxdx=f(x)expj2(-u)xdx=f(x)exp-j2(-u)x*dx（取共轭复数）=F*(-u),第三节 频域变换:二维傅立叶变换特性,旋转特性旋转特性描述：如果f(x,y)旋转了一个角度，那么f(x,y)旋转后的图象的傅立叶变换也旋转了相同的角度。设：a(x,y)=x cos()-y sin()b(x,y)=x sin()+y cos()Ff(a(x,y),b(x,y)F(a(u,v),b(u,v),第三节 频域变换:二维傅立叶变换特性,旋转特性结论：对图象的旋转变换和傅立叶变换的顺序是可交换的FRf(x,y)RFf(x,y),第三节 频域变换:二维傅立叶变换特性,线性与相似性线性的描述：傅立叶变换是线性系统、函数和的傅立叶变换是可分离的设：f(x,y)的傅立叶变换为Ff(x,y)g(x,y)的傅立叶变换为Fg(x,y)有：Ff(x,y)+g(x,y)=Ff(x,y)+Fg(x,y),第三节 频域变换:二维傅立叶变换特性,线性的证明用一维函数为例进行证明：Ff(x)+g(x)=F(u)=(f(x)+g(x)exp-j2uxdx=(f(x)exp-j2ux+g(x)exp-j2ux)dx=f(x)exp-j2uxdx+g(x)exp-j2uxdx=F(u)+G(u),第三节 频域变换:二维傅立叶变换特性,线性与相似性相似性的描述：a f(x,y)a F(u,v)且有：f(ax,by)1/|ab|F(u/a,v/b),第三节 频域变换:二维傅立叶变换特性,相似性的证明用一维函数为例进行证明：f(ax)1/|a|F(u/a)Ff(ax)=F(u)=f(ax)exp-j2uxdx将指数和积分同时乘以 1=a/a并设：v=ax,dv=adxFf(ax)=f(ax)exp-j2ux a/a a/a dx=1/a f(ax)exp-j2u xa/a adx=1/a f(v)exp-j2v(u/a)dv=1/|a|F(u/a),第三节 频域变换:二维傅立叶变换特性,均值性均值性的描述：离散函数的均值等于该函数傅立叶变换在(0,0)点的值 M-1N-1 f(x,y)=1/MNf(x,y)e0 x=0 y=0f(x,y)=F(0,0),第三节 频域变换:二维傅立叶变换特性,拉普拉斯拉普拉斯特性的描述：给出二维拉普拉斯函数的傅立叶变换表达式：拉普拉斯函数：2 f(x,y)=2f/x2+2f/y2其傅立叶变换为：F2 f(x,y)=-42(u2+v2)F(u,v)这个定理将在图象的边界提取中用到,第三节 频域变换:二维傅立叶变换特性,卷积与相关:空域和频域之间的基本联系卷积定理的描述：空域中的卷积等价于频域中的相乘f(x,y)*g(x,y)F(u,v)G(u,v)Ff(x,y)*g(x,y)=F(u,v)G(u,v)同时有：f(x,y)g(x,y)F(u,v)*G(u,v),第三节 频域变换:二维傅立叶变换特性,卷积定理成立的证明用一维函数为例进行证明：Ff(x)*g(x)=f(x)*g(x)exp-j2uxdx=f(t)g(x-t)dt exp-j2uxdx 对于上面这个式子，我们可以看出是一个两个变量t、x的二维积分。交换积分的顺序有：,第三节 频域变换:二维傅立叶变换特性,=f(t)g(x-t)exp-2juxdxdt=f(t)g(x-t)exp-2juxdxdt将 t 视为位移量，由平移定理Gg(x-t)=g(x-t)exp-2juxdx=exp-j2tuG(u)有：=f(t)exp-2jtuG(u)dt=G(u)f(t)exp-2jtu dt=F(u)G(u),第三节 频域变换:二维傅立叶变换特性,卷积与相关:空域和频域之间的基本联系相关定理的描述：空域中f(x,y)与g(x,y)的相关等价于频域中F(u,v)的共轭与G(u,v)相乘f(x,y)g(x,y)F*(u,v)G(u,v)同时有：f*(x,y)g(x,y)F(u,v)G(u,v),第三节 频域变换:快速傅立叶变换,2.2.3 快速傅立叶变换FFT算法逆向FFT算法算法实现,第三节 频域变换:快速傅立叶变换,FFT算法基本思想 FFT算法基于一个叫做递推加倍的方法。为方便起见我们用下式表达离散傅立叶变换公式 N-1F(u)=1/Nf(x)(WN)ux x=0 这里 WN=exp(-j2/N)是一个常数,第三节 频域变换:快速傅立叶变换,FFT算法基本思想 假设N为：N=2n 其中n是一个正整数，因此N可表示为：N=2M 这里M仍然是一个正整数。将N=2M代入上式，得到：,第三节 频域变换:快速傅立叶变换,FFT算法基本思想 2M-1 F(u)=1/(2M)f(x)(W2M)ux x=0 M-1 M-1=1/2 1/Mf(2x)W2Mu(2x)+1/Mf(2x+1)W2Mu(2x+1)x=0 x=0,第三节 频域变换:快速傅立叶变换,FFT算法基本思想由于：WN=exp(-j2/N)W2M2ux=exp-j22ux/2M=exp-j2ux/M=WMux所以：W2M2xu=Wmxu代入上式有：,第三节 频域变换:快速傅立叶变换,M-1 M-1 1/21/Mf(2x)Wmux+1/Mf(2x+1)WMux W2Mu x=0 x=0定义两个符号：M-1 Feven(u)=1/Mf(2x)Wmux 偶数部分 x=0u=0,1,2,M-1 M-1 Fodd(u)=1/Mf(2x+1)Wmux 奇数部分 x=0 u=0,1,2,M-1,第三节 频域变换:快速傅立叶变换,得出FFT 的第一个递推公式：F(u)=1/2(Feven(u)+Fodd(u)W2Mu)该公式说明F(u)可以通过奇部和偶部之和来计算,第三节 频域变换:快速傅立叶变换,另有：WMu+M=exp-2j(u+M)/M=exp-2ju/Mexp-2j=WMuej(-2)=WMu(-1)(-2)=Wmu且：W2Mu+M=exp-2j(u+M)/2M=exp-2ju/2M ej(-1)=W2Mu(-1)(-1)=-W2Mu最后有：WMu+M=Wmu；W2Mu+M=-W2Mu,第三节 频域变换:快速傅立叶变换,因为：WMu+M=Wmu；W2Mu+M=-W2Mu得出u+M 的DFT为：M-1F(u+M)=1/2 1/Mf(2x)WM(u+M)x+x=0 M-1 1/Mf(2x+1)WM(u+M)x W2Mu+M x=0=1/2(Feven(u)-Fodd(u)W2Mu),第三节 频域变换:快速傅立叶变换,得出FFT的第二个递推公式：F(u+M)=1/2(Feven(u)-Fodd(u)W2Mu)该公式说明F(u+M)可以通过F(u)偶部和奇部之差来计算,第三节 频域变换:快速傅立叶变换,得出FFT的两个递推公式：F(u)=1/2(Feven(u)+Fodd(u)W2Mu)F(u+M)=1/2(Feven(u)-Fodd(u)W2Mu),第三节 频域变换:快速傅立叶变换,分析这些表达式得到如下一些有趣的特性：（1）一个N个点的变换，能够通过将原始 表达 式分成两个部分来计算（2）通过计算两个（N/2）个点的变换。得到 Feven(u)和 Fodd(u)（3）奇部与偶部之和得到F(u)的前(N/2)个值。（4）奇部与偶部之差得到F(u)的后(N/2)个值。且不需要额外的变换计算。,第三节 频域变换:快速傅立叶变换,归纳快速傅立叶变换的思想：1）通过计算两个单点的DFT，来计算两个点的DFT，2）通过计算两个双点的DFT，来计算四个点的DFT，以此类推3）对于任何N=2m的DFT的计算，通过计算两个N/2点的DFT，来计算N个点的DFT,第三节 频域变换:快速傅立叶变换,逆向FFT算法算法思想描述：用正向变换计算逆向变换 N-1F(u)=1/N f(x)exp-j2ux/N x=0u=0,1,2,.N-1 N-1 f(x)=F(u)expj2ux/N u=0 x=0,1,2,.N-1,第三节 频域变换:快速傅立叶变换,逆向FFT算法 在离散逆向变换表达式两边同取共轭，并除N N-11/Nf*(x)=1/N F*(u)exp-j2ux/N u=0u=0,1,2,.N-1 用正向变换算法计算，得到1/Nf*(x)，取共轭并乘上N，即得到f(x),第三节 频域变换:快速傅立叶变换,FFT算法实现通过一个实例来体会一下FFT算法：设：有函数f(x)，其 N=23=8 有：f(0),f(1),f(2),f(3),f(4),f(5),f(6),f(7)计算：F(0),F(1),F(2),F(3),F(4),F(5),F(6),F(7),第三节 频域变换:快速傅立叶变换,FFT算法实现首先分成奇偶两组：有：f(0),f(2),f(4),f(6)f(1),f(3),f(5),f(7)为了利用递推特性，再分成两组：有：f(0),f(4)，f(2),f(6)f(1),f(5)，f(3),f(7),第三节 频域变换:快速傅立叶变换,f(0),f(4)f(2),f(6)f(1),f(5)f(3),f(7)F2(0),F2(4)F2(2),F2(6)F2(1),F2(5)F2(3),F2(7)F4(0),F4(4),F4(2),F4(6)F4(1),F4(5),F4(3),F4(7)F8(0),F8(1),F8(2),F8(3),F8(4),F8(5),F8(6),F8(7),第三节 频域变换:快速傅立叶变换,算法实现的几个关键点1）地址的排序：按位倒序规则例如：N=23=8原地址 原顺序 新地址 新顺序 000 f(0)000 f(0)001 f(1)100 f(4)010 f(2)010 f(2)011 f(3)110 f(6)100 f(4)001 f(1)101 f(5)101 f(5)110 f(6)011 f(3)111 f(7)111 f(7),第三节 频域变换:快速傅立叶变换,算法实现的几个关键点2）计算顺序及地址增量 地址+1 地址+2 地址+4 f(0)F2(0)F4(0)f(4)F2(4)F4(4)f(2)F2(2)F4(2)f(6)F2(6)F4(6)f(1)F4(1)F4(1)f(5)F2(5)F4(5)f(3)F2(3)F4(3)f(7)F2(7)F4(7),第三节 频域变换:快速傅立叶变换,算法实现的几个关键点3）复系数的计算：尤拉公式 W2M=exp-j2/2M=exp-j/M=cos(/M)+jsin(/M),第三节 频域变换:快速傅立叶变换,SUBROUTINE FFT(F,LN)COMPLEX F(1024),U,W,T,CMPLXPI=3.1415926N=2*LN/*要计算FFT的函数点数*/NV2=N/2NM1=N-1J=1,第三节 频域变换:快速傅立叶变换,DO 3 I=1,NM1IF(I.GE.J)GOTO 1T=F(J)F(J)=F(I)T=F(I)1K=NV22IF(K.GE.J)GOTO 3 K=K/2GOTO 23 J=J+K/*交换输入函数F(I)的顺序*/,第三节 频域变换:快速傅立叶变换,DO 5 L=1,LNLE=2*LLE1=LE/2/*地址增量计算*/U=(1.0,0.0)/*系数赋初值*/W=CMPLX(COS(PI/LE1),-SIN(PI/LE1)DO 5 J=1,LE1 DO 4 I=J,N,LE IP=I+LE1/*计算地址*/T=F(IP)*U/*奇部乘系数*/F(IP)=F(I)-T/*后半部分计算,第三节 频域变换:快速傅立叶变换,4F(I)=F(I)+T/*后半部分计算*/5U=U*W/*新递推系数计算*/DO 6 I=1,N 6 F(I)=F(I)/FLOAT(N)RETURNEND,第二次作业 FFT算法的实现,1.将宽为2n的正方形图象，用FFT算法从空域变 换到频域；2.将频域图象以中心为原点的四个象限，做水平和垂直镜像，使图象能量中心，对应到几何中心，并用频域图象的模来进行显示。3.将频域图象，通过FFT逆变换到空域，并显示。,请提问,