1、图像的傅里叶变换图像的傅里叶变换 Fourier Transformation For Image 时域分析只能反映信号的幅值随时间的变化时域分析只能反映信号的幅值随时间的变化情况,除单频率分量的简谐波外,很难明确揭示情况,除单频率分量的简谐波外,很难明确揭示信号的频率组成和各频率分量大小。信号的频率组成和各频率分量大小。图例:受噪声干扰的多频率成分信号图例:受噪声干扰的多频率成分信号 时间时间幅值幅值频率频率时域分析时域分析频域分析频域分析信号频谱信号频谱X(f)X(f)代表了信号在不同频率分量成分的大小,能够提代表了信号在不同频率分量成分的大小,能够提供比时域信号波形更直观,丰富的信息。供
2、比时域信号波形更直观,丰富的信息。一维一维FT及其反变换及其反变换n连续函数连续函数f(x)的傅立叶变换的傅立叶变换F(u):n傅立叶变换傅立叶变换F(u)的反变换的反变换:dxexfuFuxj2)()(dueuFxfuxj2)()(一一维维DFTDFT及其及其反变换反变换n离散函数离散函数f(x)(其中其中x,u=0,1,2,N-1)的傅立叶变换的傅立叶变换:10/2)()(NxNuxjexfuF10/2)(1)(NxNuxjeuFNxfF(u)的反变换的反变换的反变换的反变换:计算计算F(u):1)在指数项中代入在指数项中代入u=0,然后将所有,然后将所有x 值值相加,得到相加,得到F(0
3、);2)u=1,复对所有,复对所有x 的的相加,得到相加,得到F(1);3)对所有对所有M 个个u 重复重复此过程,得到全部完此过程,得到全部完整的整的FT。n离散傅里叶变换及其反变换总存在。离散傅里叶变换及其反变换总存在。n用欧拉公式得用欧拉公式得sincosjej10/2sin/2)cos()(NxNuxjNuxxfuF每个每个F(u)由由f(x)与对应频率的正弦和余弦乘积和组成与对应频率的正弦和余弦乘积和组成;u 值决定了变换的频率成份,因此,值决定了变换的频率成份,因此,F(u)覆盖的域覆盖的域(u值值)称为称为频率域频率域,其中每一项都被称为,其中每一项都被称为FT 的的频率频率分量
4、分量。与。与f(x)的的“时间域时间域”和和“时间成份时间成份”相对应。相对应。傅里叶变换的作用傅里叶变换的作用傅里叶变换将信号分成不同频率成份傅里叶变换将信号分成不同频率成份。类似光学中的。类似光学中的分色棱镜把白光按波长分色棱镜把白光按波长(频率频率)分成不同颜色,称数学分成不同颜色,称数学棱镜棱镜。傅里叶变换的成份:傅里叶变换的成份:直流分量和交流分量直流分量和交流分量信号变化的快慢与频率域的频率有关信号变化的快慢与频率域的频率有关。噪声、边缘、噪声、边缘、跳跃部分代表图像跳跃部分代表图像的的高频高频分量;分量;背景区域和慢变部分背景区域和慢变部分代表图像代表图像的的低频低频分量分量二维
5、二维DFT傅里叶变换傅里叶变换n一个图像尺寸为一个图像尺寸为MN的函数的函数f(x,y)的离散傅立叶变换的离散傅立叶变换F(u,v):nF(u,v)的反变换的反变换:112(/)00(,)(,)MNjux Mvy NxyF u vf x y e112(/)001(,)(,)MNjux Mvy Nuvf x yF u v eMN二维二维DFT傅里叶变换傅里叶变换n(u,v)=(0,0)位置的傅里叶变换值为位置的傅里叶变换值为1010),(),(1)0,0(MxNyyxfyxfMNF即即f(x,y)的均值,原点的均值,原点(0,0)的傅里叶变换是图像的的傅里叶变换是图像的平均灰度平均灰度。F(0,
6、0)称为频率谱的称为频率谱的直流分量直流分量(系数系数),其它其它F(u,v)值称为值称为交流分量交流分量(交流系数交流系数)。n 二维连续傅里叶变换二维连续傅里叶变换1)定义定义dxexfuFuxj2)()(dydxeyxfvuFvyuxj)(2),(),(2)逆傅里叶变换逆傅里叶变换dueuFxfuxj2)()(dvduevuFyxfvyuxj)(2),(),(3)傅里叶变换特征参数傅里叶变换特征参数),(),(),(vujIvuRvuF频谱频谱/幅度谱幅度谱/模模),(),(),(22vuIvuRvuF能量谱能量谱/功率谱功率谱),(),(),(),(222vuIvuRvuFvuP相位谱
7、相位谱),(),(arctan),(vuRvuIvu傅里叶变换中出现的变量傅里叶变换中出现的变量u和和v通常称为频率变量,通常称为频率变量,空间空间频率频率可以理解为等相位线在可以理解为等相位线在x,y坐标投影的截距的倒数。坐标投影的截距的倒数。xy0XY相应的空间频率分别为相应的空间频率分别为cos1,cos1YvXu 对图像信号而言,空间频率是指单位长度内亮度作对图像信号而言,空间频率是指单位长度内亮度作周期性变化的次数。周期性变化的次数。思考:噪声、线、细节、思考:噪声、线、细节、背景或平滑区域对应的空背景或平滑区域对应的空间频率特性?间频率特性?傅里叶变换的意义傅里叶变换的意义傅里叶变
8、换好比一个玻璃棱镜傅里叶变换好比一个玻璃棱镜棱镜是可以将光分成不同颜色的物理仪棱镜是可以将光分成不同颜色的物理仪器,每个成分的颜色由波长决定。器,每个成分的颜色由波长决定。傅里叶变换可看做是傅里叶变换可看做是“数学中的棱镜数学中的棱镜”,将函数基于频率分成不同的成分。将函数基于频率分成不同的成分。一些图像的傅里叶变换一些图像的傅里叶变换 是是g(x,y)的频谱,的频谱,物函数物函数g(x,y)可以看作不同方可以看作不同方向传播的单色平面波分量的线性叠加向传播的单色平面波分量的线性叠加。为权重因子。空间频率为权重因子。空间频率 表示了单色表示了单色平面波的传播方向。平面波的传播方向。(,)G(,
9、)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离散的情况下,傅里叶变换和逆傅里叶变换始终存在。离散的情况下,傅里叶变换和逆傅里叶变换始终存在。例例 设一函数如图(设一函数如图(a)
10、所示,如果将此函数在自变量)所示,如果将此函数在自变量25.1,00.1,75.0,5.03210 xxxx并重新定义为图(并重新定义为图(b)离散函数,求其傅里叶变换。)离散函数,求其傅里叶变换。取样取样(a)(b))3()2()1()0(111111111111)(41ffffjjjjuFxy1-1j-j图像的频谱幅度随频率增大而迅速衰减图像的频谱幅度随频率增大而迅速衰减 许多图像的傅里叶频谱的幅度随着频率的增大而迅速减小,这使许多图像的傅里叶频谱的幅度随着频率的增大而迅速减小,这使得在显示与观察一副图像的频谱时遇到困难。但以图像的形式显示它得在显示与观察一副图像的频谱时遇到困难。但以图像
11、的形式显示它们时,其高频项变得越来越不清楚。们时,其高频项变得越来越不清楚。解决办法:解决办法:对数化对数化2526主极大的值用主极大的值用Fmax表示,表示,第一个旁瓣的峰值用第一个旁瓣的峰值用Fmin表示表示)1(log)1(logminmaxKFKFRaa例题:对一幅图像实施二维例题:对一幅图像实施二维DFT,显示并观察其频谱。,显示并观察其频谱。解:源程序及运行结果如下:解:源程序及运行结果如下:%对单缝进行快速傅里叶变换,以三种方式显示频谱,对单缝进行快速傅里叶变换,以三种方式显示频谱,%即:直接显示(坐标原点在左上角);把坐标原点平即:直接显示(坐标原点在左上角);把坐标原点平%移
12、至中心后显示;以对数方式显示。移至中心后显示;以对数方式显示。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);%把频谱坐标原点由左上角移至屏幕中央把频谱坐标原点由左上角移至屏幕中央subplot(223)Fd=abs(Fc);ims
13、how(Fd,)ratio=max(Fd(:)/min(Fd(:)%ratio=2.3306e+007,动态范围太大,显示器无法正常显动态范围太大,显示器无法正常显示示title(幅度谱(频谱坐标原点在屏幕中央)幅度谱(频谱坐标原点在屏幕中央))S2=log(1+abs(Fc);subplot(224)imshow(S2,)title(以对数方式显示频谱以对数方式显示频谱)运行上面程序后,结果如下:运行上面程序后,结果如下:单 狭 缝 图 像幅 度 谱(频 谱 坐 标 原 点 在 坐 上 角)幅 度 谱(频 谱 坐 标 原 点 在 屏 幕 中 央)以 对 数 方 式 显 示 频 谱n 二维离散
14、傅里叶变换的性质二维离散傅里叶变换的性质 线性性线性性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%该程序验证了二维该程序验证了二维DFT的线性性质的线性性质f=imread(D:chenpcdatathrychpt4Fig4.04(a).jpg)
15、;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=fftshift(fft2(f+g);title(DFT(f)+DFT(g)subplot(224)imshow(log(abs(FG)
16、,)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 u v eeMNF u y eM120120,vyNjNyuxMjMxF x vf x y eyDFTF u vF x v exDFT
17、方向的方向的1201201,1,vyNjNvuxMjMuF u yF u v eyIDFTNf x yF u y exIDFTM方向的方向的其中:其中:例题:编程验证二维离散傅里叶变换可分离为两个一维离例题:编程验证二维离散傅里叶变换可分离为两个一维离散傅里叶变换。散傅里叶变换。解:解:%myseparable.m%该程序验证了二维该程序验证了二维DFT的可分离性质的可分离性质%该程序产生了冈萨雷斯该程序产生了冈萨雷斯数字图像处理数字图像处理(第二版)(第二版)%P125 图图4.4f=imread(D:chenpcdatathrychpt4Fig4.04(a).jpg);subplot(21
18、1)imshow(f,)title(原图原图)F=fftshift(fft2(f);subplot(223)imshow(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实 现 二 维 离 散 傅 里 叶 变 换用
19、 fft实 现 二 维 离 散 傅 里 叶 变 换 平移性平移性0000200200,u xv yjMNuxvyjMNf 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)移至显示屏的中心点只要将移
20、至显示屏的中心点只要将f(x,y)乘以乘以(-1)x+y因子再进行傅因子再进行傅里叶变换即可实现。里叶变换即可实现。例题:利用例题:利用(-1)x+y对单缝图像对单缝图像f(x,y)进行调制,实现把频谱进行调制,实现把频谱坐标原点移至屏幕正中央的目标。坐标原点移至屏幕正中央的目标。2,200NvMu当当yxyxjNyvMxujee)1()()/(200解:完成本题的源程序为:解:完成本题的源程序为:%在傅里叶变换之前,把函数乘以在傅里叶变换之前,把函数乘以(-1)x+y,相当于把频谱,相当于把频谱%坐标原点移至屏幕窗口正中央。坐标原点移至屏幕窗口正中央。f(512,512)=0;f=mat2g
21、ray(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=fft2(g);subplot(224),imshow(log(1+abs(G),),title(g(x
22、,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)曲线进行调制,达到平移频域坐标原点曲线进行调制,达到平移频域坐标原点至屏幕正中央的目的。至屏幕正中央
23、的目的。%以一维情况为例,说明空域调制对应着频域坐标原点移位。以一维情况为例,说明空域调制对应着频域坐标原点移位。f(1:512)=0;f(251:260)=1;%产生宽度为产生宽度为10的窗口函数的窗口函数subplot(221),plot(f),title(宽度为宽度为10 的窗口函数的窗口函数)F=fft(f,512);%进行快速傅里叶变换,延拓周期周期为进行快速傅里叶变换,延拓周期周期为512subplot(222)plot(abs(F)%绘幅度频谱(频谱坐标原点在左边界处)绘幅度频谱(频谱坐标原点在左边界处)title(幅度谱(频谱坐标原点在左边界处)幅度谱(频谱坐标原点在左边界处)
24、)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(幅度谱(频谱坐标原点在中央)幅度谱(频谱坐标原点在中央))figuref(1:512)=0;f(251:270)=1;%产生
25、宽度为产生宽度为20的窗口函数的窗口函数subplot(221),plot(f),title(宽度为宽度为20 的窗口函数的窗口函数)F=fft(f,512);%进行快速傅里叶变换,延拓周期周期为进行快速傅里叶变换,延拓周期周期为512subplot(222)plot(abs(F)%绘幅度频谱(频谱坐标原点在左边界处)绘幅度频谱(频谱坐标原点在左边界处)title(幅度谱(频谱坐标原点在左边界处)幅度谱(频谱坐标原点在左边界处))x=251:270;f(251:270)=(-1).x;%把曲线把曲线f(x)乘以乘以(-1)x,可以把频谱坐标原点移至,可以把频谱坐标原点移至屏幕正中央屏幕正中央s
26、ubplot(223),plot(f),title(宽度为宽度为20 的调制窗口函数的调制窗口函数)F=fft(f,512);%进行快速傅里叶变换进行快速傅里叶变换subplot(224);plot(abs(F)%直接显示幅度频谱(频谱坐标原点在正中央)直接显示幅度频谱(频谱坐标原点在正中央)title(幅度谱(频谱坐标原点在中央)幅度谱(频谱坐标原点在中央))020040060000.51宽 度 为 10 的 窗 口 函 数02004006000510幅 度 谱(频 谱 坐 标 原 点 在 左 边 界 处)0200400600-1-0.500.51宽 度 为 10 的 调 制 窗 口 函 数
27、02004006000510幅 度 谱(频 谱 坐 标 原 点 在 中 央)020040060000.51宽 度 为 20 的 窗 口 函 数020040060005101520幅 度 谱(频 谱 坐 标 原 点 在 左 边 界 处)0200400600-1-0.500.51宽 度 为 20 的 调 制 窗 口 函 数020040060005101520幅 度 谱(频 谱 坐 标 原 点 在 中 央)(2)空域移位:)空域移位:),(00yyxxfDFT0000000000001120000112000011222,uxvyMNjMNxyu x xxv y yyMNjMNxyuxvyuxvyM
28、x 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 vFuv证明:证明:(1)周期性:)周期性:1120011200,1,uxvyMNjMNxyuxvyMNjMNuvF u vf x yef x yF u veMN21jme,F umM vnNF u vf xmM ynNf x y
29、(2)共轭对称性:共轭对称性:11200,ux vyMNjMNxyF u vf x y e*11200*,u xv yMNjMNxyf x y eFuv,|,|,|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
30、注注:为看清问题的实质、简化旋转不变性的证明,以:为看清问题的实质、简化旋转不变性的证明,以上用二维连续傅里叶变换进行证明。实际上,由连续上用二维连续傅里叶变换进行证明。实际上,由连续积分公式进行离散化处理,即可得到离散公式,证明积分公式进行离散化处理,即可得到离散公式,证明可参照连续情况进行。可参照连续情况进行。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(原图的频谱原图的频谱)
31、f=imrotate(f,45,bilinear,crop);subplot(223)imshow(f,)title(旋转旋转450图图)Fc=fftshift(fft2(f);subplot(224);imshow(log(1+abs(Fc),)title(旋转图的频谱旋转图的频谱)原 图原 图 的 频 谱旋 转 450图旋 转 图 的 频 谱 离散卷积定理离散卷积定理daaxgafxgxf)()()(*)(例例1 其它0101)(xxf其它0102/1)(xxg求以下两个函数的卷积求以下两个函数的卷积1)连续卷积)连续卷积dudvvyuxgvufyxgyxf),(),(),(*),(其它0
32、212/1102/)(*)(xxxxxgxf2 2)离散卷积定理)离散卷积定理1010),(),(),(*),(MmNnnymxgnmfyxgyxf 19)()(*)(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 Tfxygx
33、y1100110011200,MNmnMNmnm un vMNjMNmnD F TfmngxmynfmnD F TgxmynfmneGu vFu vGu 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 离散的卷积原理基本上是和连续卷积相同,其差别离散的卷积原理基本上是和连续卷积相同,其差别仅仅是在与抽样间隔对应的离散增增量处发生位移,仅仅是在与抽样间隔对应的离散增增量
34、处发生位移,用求和代替微分,用求和代替微分,由于离散傅里叶变换和它的逆傅里由于离散傅里叶变换和它的逆傅里叶变换都是周期函数,叶变换都是周期函数,那么离散卷积定理应该和这个那么离散卷积定理应该和这个周期联系起来,周期联系起来,就是让在计算卷积时让这两个离散函就是让在计算卷积时让这两个离散函数具有同样的周期,否则将产生错误。数具有同样的周期,否则将产生错误。注意注意:利用:利用FFT计算卷积时,为防止频谱混叠误差,计算卷积时,为防止频谱混叠误差,需对离散的二维函数补零,即需对离散的二维函数补零,即周期延拓周期延拓,对两个函数,对两个函数同时添加零,使它们具有相同的周期。同时添加零,使它们具有相同的
35、周期。0200400)(mf20080030200400)(mh28000200400)(mh 2800周期延拓周期延拓周期延拓周期延拓),(yxf 的大小为的大小为BA),(yxg 的大小为的大小为DC1,1010,10),(),(NyBMxAByAxyxfyxfe1,1010,10),(),(NyDMxCDyCxyxgyxge1,1DCNCAM),(*),(),(yxgyxfyxzeee空间域滤波和频域滤波的关系空间域滤波和频域滤波的关系1010)0,0(),(),(MxNysyxyxsMNeyxMNMxNyNvyMuxj1),(11010)/(2MNyxhnymxhnmyxhyxMmNn
36、),(),(),(),(*),(1010),(),(),(*),(vuHvuFyxhyxf),(),(*),(vuHyxhyx),(),(vuHyxh空间域和频域空间域和频域的滤波器构成的滤波器构成傅里叶变换对傅里叶变换对 相关定理相关定理ddyxgfyxgyxf),(),(),(),(),(*),(),(),(vuGvuFyxgyxf),(),(),(*),(vuGvuFyxgyxf证明:证明:,D F Tfxygxy1100110011200*11200*,MNmnMNmnm unvMNjMNmnm unvMNjMNmnD F Tfm ngxmynfm nD F Tgxmynfm neGu vfm neGu vFu vGu v