第七讲图像的傅里叶变换

1、第七讲第七讲 傅里叶变换傅里叶变换 Lecture 7 Fourier Transformation重点内容重点内容n 傅里叶变换的定义傅里叶变换的定义n 空间频率空间频率(Spatial Frequency)的理解的理解n 离散傅里叶变换离散傅里叶变换(DFTDiscrete Fourier Transformation)的计算的计算思考问题：思考问题： 为什么要进行傅里叶变换为什么要进行傅里叶变换 什么是图像什么是图像 傅里叶（傅里叶（FourierFourier），法国数学及物理学家，傅里叶级数（三角级），法国数学及物理学家，傅里叶级数（三角级数）创始人。数）创始人。 18011801年

2、任伊泽尔省地方长官，年任伊泽尔省地方长官，18171817年当选科学院院士，年当选科学院院士，18221822年任年任该院终身秘书，后又任法兰西学院终身秘书和理工科大学校务委员会该院终身秘书，后又任法兰西学院终身秘书和理工科大学校务委员会主席。主席。 主要贡献：在研究热的传播时创立了一套数学理论，主要贡献：在研究热的传播时创立了一套数学理论，18071807年向巴年向巴黎科学院呈交了黎科学院呈交了热的传播论文热的传播论文，推导著名的热传导方程，并在求，推导著名的热传导方程，并在求解该方程时发现函数可由三角函数构成的级数形式表示，从而提出任解该方程时发现函数可由三角函数构成的级数形式表示，从而提

3、出任意函数可以展成三角函数的无穷级数。意函数可以展成三角函数的无穷级数。 数学与图像处理数学与图像处理 空域与频域的桥梁空域与频域的桥梁 第二种语言第二种语言n 二维连续傅里叶变换二维连续傅里叶变换1） 定义定义dxexfuFuxj2)()(dydxeyxfvuFvyuxj)(2),(),(2） 逆傅里叶变换逆傅里叶变换dueuFxfuxj2)()(dvduevuFyxfvyuxj)(2),(),(3） 傅里叶变换特征参数傅里叶变换特征参数),(),(),(vujIvuRvuF频谱频谱/模模),(),(),(22vuIvuRvuF能量谱能量谱/ /功率谱功率谱),(),(),(),(222vu

4、IvuRvuFvuP相位角相位角),(),(arctan),(vuRvuIvu傅里叶变换中出现的变量傅里叶变换中出现的变量u和和v通常称为频率变量，通常称为频率变量，空间空间频率频率可以理解为等相位线在可以理解为等相位线在x,y坐标投影的截距的倒数。坐标投影的截距的倒数。xy0XY相应的空间频率分别为相应的空间频率分别为cos1,cos1YvXu 对图像信号而言，空间频率是指单位长度内亮度作对图像信号而言，空间频率是指单位长度内亮度作周期性变化的次数。周期性变化的次数。思考：噪声、线、细节、思考：噪声、线、细节、背景或平滑区域对应的空背景或平滑区域对应的空间频率特性？间频率特性？傅里叶变换的意

5、义傅里叶变换的意义傅里叶变换好比一个玻璃棱镜傅里叶变换好比一个玻璃棱镜棱镜是可以将光分成不同颜色的物理仪棱镜是可以将光分成不同颜色的物理仪器，每个成分的颜色由波长决定。器，每个成分的颜色由波长决定。傅里叶变换可看做是傅里叶变换可看做是“数学中的棱镜数学中的棱镜”，将函数基于频率分成不同的成分。将函数基于频率分成不同的成分。一些图像的傅里叶变换一些图像的傅里叶变换选自MATLAB 是是g(x,y)的频谱，的频谱，物函数物函数g(x,y)可以看作不同方可以看作不同方向传播的单色平面波分量的线性叠加向传播的单色平面波分量的线性叠加。 为权重因子。空间频率为权重因子。空间频率 表示了单色表示了单色平面

6、波的传播方向。平面波的传播方向。( , )G ( , )Gd d coscos,对于对于xy平面上一点的复振幅分布平面上一点的复振幅分布g(x,y)可由逆傅里叶可由逆傅里叶变换表示成：变换表示成：( , )( , )exp 2 ()g x yGjxy d d n 二维离散傅里叶变换二维离散傅里叶变换1） 定义定义1010)/(2),(1),(MxNyNvyMuxjeyxfMNvuF2） 逆傅里叶变换逆傅里叶变换1010)/(2),(),(MuNvNvyMuxjevuFyxf1, 1 , 01, 1 , 0NvMu1, 1 , 01, 1 , 0NyMx离散的情况下，傅里叶变换和逆傅里叶变换始终

7、存在。离散的情况下，傅里叶变换和逆傅里叶变换始终存在。例例 设一函数如图（设一函数如图（a）所示，如果将此函数在自变量）所示，如果将此函数在自变量25. 1,00. 1,75. 0, 5 . 03210 xxxx并重新定义为图（并重新定义为图（b）离散函数，求其傅里叶变换。）离散函数，求其傅里叶变换。取样取样（a）（b）)3()2() 1 ()0(111111111111)(41ffffjjjjuFxy1-1j-j图像的频谱幅度随频率增大而迅速衰减图像的频谱幅度随频率增大而迅速衰减 许多图像的傅里叶频谱的幅度随着频率的增大而迅速减小，这使许多图像的傅里叶频谱的幅度随着频率的增大而迅速减小，这使

8、得在显示与观察一副图像的频谱时遇到困难。但以图像的形式显示它得在显示与观察一副图像的频谱时遇到困难。但以图像的形式显示它们时，其高频项变得越来越不清楚。们时，其高频项变得越来越不清楚。解决办法：解决办法： 对数化对数化主极大的值用主极大的值用Fmax表示，表示，第一个旁瓣的峰值用第一个旁瓣的峰值用Fmin表示表示)1 (log)1 (logminmaxKFKFRaa例题：对一幅图像实施二维例题：对一幅图像实施二维DFT，显示并观察其频谱。，显示并观察其频谱。解：源程序及运行结果如下：解：源程序及运行结果如下：%对单缝进行快速傅里叶变换，以三种方式显示频谱，对单缝进行快速傅里叶变换，以三种方式显

9、示频谱，%即：直接显示（坐标原点在左上角）；把坐标原点平即：直接显示（坐标原点在左上角）；把坐标原点平%移至中心后显示；以对数方式显示。移至中心后显示；以对数方式显示。f=zeros(512,512);f(246:266,230:276)=1;subplot(221),imshow(f,),title(单狭缝图像单狭缝图像)F=fft2(f); %对图像进行快速傅里叶变换对图像进行快速傅里叶变换S=abs(F); subplot(222)imshow(S,) %显示幅度谱显示幅度谱title(幅度谱（频谱坐标原点在坐上角）幅度谱（频谱坐标原点在坐上角）)Fc=fftshift(F); %把频谱

10、坐标原点由左上角移至屏幕中央把频谱坐标原点由左上角移至屏幕中央subplot(223)Fd=abs(Fc);imshow(Fd,)ratio=max(Fd(:)/min(Fd(:) %ratio = 2.3306e+007,动态范围太大，显示器无法正常显动态范围太大，显示器无法正常显示示title(幅度谱（频谱坐标原点在屏幕中央）幅度谱（频谱坐标原点在屏幕中央）)S2=log(1+abs(Fc);subplot(224)imshow(S2,)title(以对数方式显示频谱以对数方式显示频谱)运行上面程序后，结果如下：运行上面程序后，结果如下：单 狭 缝 图 像幅 度 谱 （ 频 谱 坐 标 原

11、 点 在 坐 上 角 ）幅 度 谱 （ 频 谱 坐 标 原 点 在 屏 幕 中 央 ）以 对 数 方 式 显 示 频 谱n 二维离散傅里叶变换的性质二维离散傅里叶变换的性质 线性性线性性111 122112222,fx yF u vc fx yc fx yc F u vc F u vfx yF u v证明：证明：1122,DFTc fx yc fx y112112200111122112200001122,uxvyMNjMNxyuxvyuxvyMNMNjjMNMNxyxyc fx yc fx yecfx y ecfx y ec Fu vc Fu v %imagelinear.m%该程序验证了二

12、维该程序验证了二维DFT的线性性质的线性性质f=imread(D:chenpcdatathrychpt4Fig4.04(a).jpg);g=imread(D:chenpcdatathrychpt4Fig4.30(a).jpg);m,n=size(g);f(m,n)=0;f=im2double(f);g=im2double(g);subplot(221)imshow(f,)title(f)subplot(222)imshow(g,)title(g)F=fftshift(fft2(f);G=fftshift(fft2(g);subplot(223)imshow(log(abs(F+G),)FG=f

13、ftshift(fft2(f+g);title(DFT(f)+DFT(g)subplot(224)imshow(log(abs(FG),)title(DFT(f+g)fgDFT(f)+DFT(g)DFT(f+g) 可分离性可分离性二维二维DFT可视为由沿可视为由沿x,y方向的两个一维方向的两个一维DFT所构成。所构成。11200,ux vyMNjMNxyF u vf x y e112200120,vyuxMNjjNMxyuxMjMxf x y eeF x v e 112001,ux vyMNjMNuvf x yF u v eMN11220012011,1,vyuxMNjjNMuvuxMjMuF

14、 u v eeMNF u y eM120120,vyNjNyuxMjMxF x vf x y eyDFTF u vF x v exDFT方向的方向的1201201,1,vyNjNvuxMjMuF u yF u v eyIDFTNf x yF u y exIDFTM方向的方向的其中：其中：The discrete Fourier transform1010)/(2),(1),(MxNyNvyMuxjeyxfMNvuFcan be expressed in the separable form10/210/210/2),(),(111),(MxMuxjNyNvyjMxMuxjevxFeyxfeMN

15、MvuFSeparability例题：编程验证二维离散傅里叶变换可分离为两个一维离例题：编程验证二维离散傅里叶变换可分离为两个一维离散傅里叶变换。散傅里叶变换。解：解：%myseparable.m%该程序验证了二维该程序验证了二维DFT的可分离性质的可分离性质%该程序产生了冈萨雷斯该程序产生了冈萨雷斯数字图像处理数字图像处理（第二版）（第二版）%P125 图图4.4f=imread(D:chenpcdatathrychpt4Fig4.04(a).jpg);subplot(211)imshow(f,)title(原图原图)F=fftshift(fft2(f);subplot(223)imshow

16、(log(1+abs(F),)title(用用fft2实现二维离散傅里叶变换实现二维离散傅里叶变换)m,n=size(f);F=fft(f); %沿沿x方向求离散傅里叶变换方向求离散傅里叶变换G=fft(F); %沿沿y方向求离散傅里叶变换方向求离散傅里叶变换F=fftshift(G);subplot(224)imshow(log(1+abs(F),)title(用用fft实现二维离散傅里叶变换实现二维离散傅里叶变换)原 图用 fft2实 现 二 维 离 散 傅 里 叶 变 换用 fft实 现 二 维 离 散 傅 里 叶 变 换 平移性平移性0000200200,u xv yjMNuxvyjM

17、Nf x yeF uu vvf x yF u vf xxyyF u ve证明：证明：（1）频域移位）频域移位002,uxvyjMND F Tfxye00001122001120000,uxvyu xv yMNjjMNMNxyuuxvvyMNjMNxyfxyeefxyeFuuvv 结论：结论：00200,uxvyjMNfxyeFuuvv ,1,22xyMNfxyFuv即如果需要将频域的坐标原点从显示屏起始点（即如果需要将频域的坐标原点从显示屏起始点（0，0）移至显示屏的中心点只要将移至显示屏的中心点只要将f(x,y)乘以乘以(-1)x+y因子再进行傅因子再进行傅里叶变换即可实现。里叶变换即可实现

18、。例题：利用例题：利用(-1)x+y对单缝图像对单缝图像f(x,y)进行调制，实现把频谱进行调制，实现把频谱坐标原点移至屏幕正中央的目标。坐标原点移至屏幕正中央的目标。2,200NvMu当当yxyxjNyvMxujee) 1()()/(200Application of translation propertyWhen u0 = M/2 and v0 = N/2, it follows thatyxyxjNyvMxujee) 1()()/(200In this case)2/, 2/() 1)(,(NvMuFyxfyxsimilarlyvuvuFNyMxf) 1)(,()2/, 2/(It i

19、s common practice to multiply the input image function by (-1)x+y prior to computing the Fourier transform./2)/2,(1)(,(The operation means to shift the origin of F(u,v) from (0,0) to coordinates (M/2,N/2).Shift factor (-1)x+y解：完成本题的源程序为：解：完成本题的源程序为：%在傅里叶变换之前，把函数乘以在傅里叶变换之前，把函数乘以(-1) x+y，相当于把频谱，相当于把频谱

20、%坐标原点移至屏幕窗口正中央。坐标原点移至屏幕窗口正中央。f(512,512)=0;f=mat2gray(f);Y,X=meshgrid(1:512,1:512);f(246:266,230:276)=1;g=f.*(-1).(X+Y);subplot(221),imshow(f,),title(原图像原图像f(x,y)subplot(222),imshow(g,),title(空域调制图像空域调制图像g(x,y)=f(x,y)*(-1)x+y)F=fft2(f);subplot(223),imshow(log(1+abs(F),),title(f(x,y)的傅里叶的傅里叶频谱频谱)G=fft

21、2(g);subplot(224),imshow(log(1+abs(G),),title(g(x,y)的傅里叶的傅里叶频谱频谱)原 图 像 f(x,y)空 域 调 制 图 像 g(x,y)=f(x,y)*(-1)x+yf(x,y)的 傅 里 叶 频 谱g(x,y)的 傅 里 叶 频 谱（a） 在0 N-1周期中有两个背靠背半周期（b） 同一区间内有一个完整的周期 这就意味着，坐标原点移到了频谱图像的中间位置，这一点十分重要，这就意味着，坐标原点移到了频谱图像的中间位置，这一点十分重要，尤其是对以后的图像显示和滤波处理。尤其是对以后的图像显示和滤波处理。例题：利用例题：利用(-1)x对对f(x

22、)曲线进行调制，达到平移频域坐标原点曲线进行调制，达到平移频域坐标原点至屏幕正中央的目的。至屏幕正中央的目的。%以一维情况为例，说明空域调制对应着频域坐标原点移位。以一维情况为例，说明空域调制对应着频域坐标原点移位。f(1:512)=0;f(251:260)=1; %产生宽度为产生宽度为10的窗口函数的窗口函数subplot(221),plot(f),title(宽度为宽度为10 的窗口函数的窗口函数)F=fft(f,512); %进行快速傅里叶变换，延拓周期周期为进行快速傅里叶变换，延拓周期周期为512subplot(222)plot(abs(F) %绘幅度频谱（频谱坐标原点在左边界处）绘幅

23、度频谱（频谱坐标原点在左边界处）title(幅度谱（频谱坐标原点在左边界处）幅度谱（频谱坐标原点在左边界处）)x=251:260;f(251:260)=(-1).x; %把曲线把曲线f(x)乘以乘以(-1)x，可以把频谱，可以把频谱 %坐标原点移至屏幕正中央坐标原点移至屏幕正中央subplot(223),plot(f),title(宽度为宽度为10 的调制窗口函数的调制窗口函数)F=fft(f,512); %进行快速傅里叶变换进行快速傅里叶变换subplot(224);plot(abs(F) %直接显示幅度频谱（频谱坐标原点在正中央）直接显示幅度频谱（频谱坐标原点在正中央）title(幅度谱（

24、频谱坐标原点在中央）幅度谱（频谱坐标原点在中央）)figuref(1:512)=0;f(251:270)=1; %产生宽度为产生宽度为20的窗口函数的窗口函数subplot(221),plot(f),title(宽度为宽度为20 的窗口函数的窗口函数)F=fft(f,512); %进行快速傅里叶变换，延拓周期周期为进行快速傅里叶变换，延拓周期周期为512subplot(222)plot(abs(F) %绘幅度频谱（频谱坐标原点在左边界处）绘幅度频谱（频谱坐标原点在左边界处）title(幅度谱（频谱坐标原点在左边界处）幅度谱（频谱坐标原点在左边界处）)x=251:270;f(251:270)=(

25、-1).x; %把曲线把曲线f(x)乘以乘以(-1)x，可以把频谱坐标原点移至，可以把频谱坐标原点移至屏幕正中央屏幕正中央subplot(223),plot(f),title(宽度为宽度为20 的调制窗口函数的调制窗口函数)F=fft(f,512); %进行快速傅里叶变换进行快速傅里叶变换subplot(224);plot(abs(F) %直接显示幅度频谱（频谱坐标原点在正中央）直接显示幅度频谱（频谱坐标原点在正中央）title(幅度谱（频谱坐标原点在中央）幅度谱（频谱坐标原点在中央）)020040060000.51宽 度 为 10 的 窗 口 函 数02004006000510幅 度 谱 （

26、 频 谱 坐 标 原 点 在 左 边 界 处 ）0200400600-1-0.500.51宽 度 为 10 的 调 制 窗 口 函 数02004006000510幅 度 谱 （ 频 谱 坐 标 原 点 在 中 央 ）020040060000.51宽 度 为 20 的 窗 口 函 数020040060005101520幅 度 谱 （ 频 谱 坐 标 原 点 在 左 边 界 处 ）0200400600-1-0.500.51宽 度 为 20 的 调 制 窗 口 函 数020040060005101520幅 度 谱 （ 频 谱 坐 标 原 点 在 中 央 ）（2）空域移位：）空域移位：),(00yyx

27、xfDFT0000000000001120000112000011222,uxvyMNjMNxyu x xxv y yyMNjMNxyuxvyuxvyMx NyjjMNMNxxyyuxvyjMNyf xxyyef xxyyeef x yeef x y 00112002,uxvyMNjMNxuxvyjMNeeF u v 周期性和共轭对称性周期性和共轭对称性 周期性：周期性：,( , ),F u vF umM vnNf x yf xmM ynN,0, 1, 2,m n 共轭对称性共轭对称性：*,|,| |,|F u vFuvF u vFuvPeriodicity and conjugate sym

28、metryThe discrete Fourier transform has the following periodicity properties: F(u, v) = F(u+M, v) = F(u, v+N) = F(u+M, v+N)The inverse transform also is periodic: f(x, y) = f(x+M, y) = f(x, y+N) = f(x+M, y+N)The idea of conjugate symmetry was introduced in previous section, and is repeated here for

29、convenience:F(u, v) = F*(-u, -v)The spectrum also is symmetric about the origin:),(),(Conjugate symmetry证明：证明： （1）周期性：）周期性：1120011200,1,uxvyMNjMNxyuxvyMNjMNuvF u vf x yef x yF u veMN21jme,F umM vnNF u vf xmM ynNf x y(2) 共轭对称性：共轭对称性：11200,ux vyMNjMNxyF u vf x y e*11200*,u xv yMNjMNxyf x y eFuv ,| |,|

30、,|F u vFuvF u v ，即关于原点对称 旋转不变性旋转不变性 ,fx yF u v,f rF cossinxryrcossinuv00,f rF 证明：证明：222cos0 0,coscossin ,sin,jux vyjrF u vfx yedxdyuxryrvFf rer dr d 若则：022cos00 0,jrFf rerdrd 0022cos0022cos00 0,jrjrf rerdrdf rerdrd ,2f rf r注注：为看清问题的实质、简化旋转不变性的证明，以：为看清问题的实质、简化旋转不变性的证明，以上用二维连续傅里叶变换进行证明。实际上，由连续上用二维连续傅里

31、叶变换进行证明。实际上，由连续积分公式进行离散化处理，即可得到离散公式，证明积分公式进行离散化处理，即可得到离散公式，证明可参照连续情况进行。可参照连续情况进行。f=zeros(512,512);f(246:266,230:276)=1;subplot(221);imshow(f,)title(原图原图)F=fftshift(fft2(f);subplot(222);imshow(log(1+abs(F),)title(原图的频谱原图的频谱)f=imrotate(f,45,bilinear,crop);subplot(223)imshow(f,)title(旋转旋转450图图)Fc=fftsh

32、ift(fft2(f);subplot(224);imshow(log(1+abs(Fc),)title(旋转图的频谱旋转图的频谱)原 图原 图 的 频 谱旋 转 450图旋 转 图 的 频 谱 离散卷积定理离散卷积定理daaxgafxgxf)()()(*)(例例1 其它0101)(xxf其它0102/1)(xxg求以下两个函数的卷积求以下两个函数的卷积1）连续卷积）连续卷积dudvvyuxgvufyxgyxf),(),(),(*),(其它0212/1102/)(*)(xxxxxgxf2 2）离散卷积定理）离散卷积定理1010),(),(),(*),(MmNnnymxgnmfyxgyxf 19

33、)()(*)(Mmeeeemxgmfxgxfaasbbttysxftswyxg),(),(),(离散卷积定义：离散卷积定义：空间滤波输出：空间滤波输出：结论：空间域进行滤波的过程就是结论：空间域进行滤波的过程就是“卷积卷积”的过程。的过程。,*,1,*,f x yg x yF u vG u vf x yF u vg x yG u vf x yg x yF u vG u vMN证明：（证明：（1）空域卷积和）空域卷积和,*,D F Tfxygxy1100110011200,MNmnMNmnm un vMNjMNmnD F TfmngxmynfmnD F TgxmynfmneGu vFu vGu

34、v （2）频域卷积和：）频域卷积和：,DFTfx ygx y112 0 0112 0 011 0 01,1,1,1,*,u xv yMNjMNuvu xv yMNjMNuvMNuvDFTFuvegx yMNFuvDFTegx yMNFuvG uuvvMNFu vG u vMN 离散的卷积原理基本上是和连续卷积相同，其差别离散的卷积原理基本上是和连续卷积相同，其差别仅仅是在与抽样间隔对应的离散增增量处发生位移，仅仅是在与抽样间隔对应的离散增增量处发生位移，用求和代替微分，用求和代替微分，由于离散傅里叶变换和它的逆傅里由于离散傅里叶变换和它的逆傅里叶变换都是周期函数，叶变换都是周期函数，那么离散卷

35、积定理应该和这个那么离散卷积定理应该和这个周期联系起来，周期联系起来，就是让在计算卷积时让这两个离散函就是让在计算卷积时让这两个离散函数具有同样的周期，否则将产生错误。数具有同样的周期，否则将产生错误。注意注意：利用：利用FFT计算卷积时，为防止频谱混叠误差，计算卷积时，为防止频谱混叠误差，需对离散的二维函数补零，即需对离散的二维函数补零，即周期延拓周期延拓，对两个函数，对两个函数同时添加零，使它们具有相同的周期。同时添加零，使它们具有相同的周期。0200400)(mf20080030200400)(mh28000200400)( mh 2800周期延拓周期延拓周期延拓周期延拓),(yxf 的

36、大小为的大小为BA),(yxg 的大小为的大小为DC1, 1010 , 10),(),(NyBMxAByAxyxfyxfe1, 1010 , 10),(),(NyDMxCDyCxyxgyxge1, 1DBNCAM),(*),(),(yxgyxfyxzeeeThe need for padding (补零补零)Some important facts that need special attention:1. For DFT, the periodicity is a mathematical by-product of the way in which the discrete Fourier transform pair is defined. Periodicity is part of the process