6第六章频率域图像增强

1、第六章第六章 频率域图像增强频率域图像增强 为什么要在频率域研究图像增强为什么要在频率域研究图像增强 l 可以利用频率成分和图像外表之间的对应关系。一些在空间域可以利用频率成分和图像外表之间的对应关系。一些在空间域 表述困难的增强任务，在频率域中变得非常普通。表述困难的增强任务，在频率域中变得非常普通。 l滤波在频率域更为直观，它可以解释空间域滤波的某些性质滤波在频率域更为直观，它可以解释空间域滤波的某些性质 l给出一个问题，寻找某个滤波器解决该问题，频率域处理对于试给出一个问题，寻找某个滤波器解决该问题，频率域处理对于试 验、迅速而全面地控制滤波器参数是一个理想工具验、迅速而全面地控制滤波器

2、参数是一个理想工具 l一旦找到一个特殊应用的滤波器，通常在空间域采用硬件实现它一旦找到一个特殊应用的滤波器，通常在空间域采用硬件实现它 频域图像增强是指通过对图像进行傅立叶变频域图像增强是指通过对图像进行傅立叶变 换，将图像从空间域变换到频域，并对图像换，将图像从空间域变换到频域，并对图像 的频率成分进行相应处理，从而实现图像增的频率成分进行相应处理，从而实现图像增 强的功能；强的功能； 傅立叶变换是频域图像增强的基础工具；傅立叶变换是频域图像增强的基础工具； 通过傅立叶变换重建可以不丢失任何信息；通过傅立叶变换重建可以不丢失任何信息； 6.16.1频率域介绍频率域介绍 6.26.2背景背景

3、傅立叶分析傅立叶分析 周期函数可以表示为不同频率周期函数可以表示为不同频率 的正弦和的正弦和/ /或余弦和的形式（傅或余弦和的形式（傅 立叶级数）立叶级数） 非周期函数可以用正弦非周期函数可以用正弦/ /或余弦或余弦 乘以加权函数的积分来表示乘以加权函数的积分来表示 （傅立叶变换）（傅立叶变换） 简单的周期运动简单的周期运动 :)sin(tAy(谐波函数) ( A为振幅, 复杂的周期运动复杂的周期运动 :)sin( 1 0n n n tnAAy tnAtnA nnnn sincoscossin 令 , 2 0 0 A a ,sin nnn Aa,cos nnn Abxt 得函数项级数得函数项级

4、数)sincos( 2 1 0 xnbxna a nn n 为角频率,为初相 ) (谐波迭加) 称上述形式的级数为称上述形式的级数为三角级数三角级数. .1傅立叶级数. 6 三角函数的傅里叶级数三角函数的傅里叶级数: 1 1 2 T )sincos()( 11 1 01 tnbtnaatf nn n 直流 分量 基波分量 n =1 谐波分量 n1 1 n 7 10 0 ).( 1 1 0 Tt t dttf T a 10 0 .cos).( 2 1 1 Tt t n dttntf T a dttntf T b Tt t n .sin).( 210 0 1 1 直流 系数 余弦分

5、量 系数 正弦分量 系数 y x 例例1. 设 f (x) 是周期为 2 的周期函 数 , 它在 上的表达式为 ), 0,1 0,1 )( x x xf 解解: 先求傅里叶系数 dcos)( 1 xnxxfan 0 0 dcos1 1 dcos) 1( 1 xnxxnx ),2,1,0(0n 将 f (x) 展成傅里叶级数. O 1 1 dsin)( 1 xnxxfbn 0 0 11 ( 1)sind1 sin d nxxnx x 0 1cos nx n 0 1cos nx n 2 1 cos n n 2 1 ( 1) n n 4 , n ,0 ,5,3,1n当 ,6,4,2n当 4 ( )s

6、in f xxx3sin 3 1 xk k ) 12sin( 12 1 (,0,2 ,)xx y x 1 1 O ),2,0,(xx 7 7sin x 9 9sin x 1) 根据收敛定理可知, 时,级数收敛于 0 2 11 2) 傅氏级数的部分和逼近 3 3sin sin 4 )( x xxf 5 5sin x 说明说明: ), 2, 1, 0(kkx当 f (x) 的情况见右图. O y x 11 指数形式的傅里叶级数指数形式的傅里叶级数 10 0 1 )( 1 1 Tt t tjn n dtetf T F tjn n enFtf 1 )()( 1 复指数傅立叶级数有一个统一的傅立叶系数，

7、复指数傅立叶级数有一个统一的傅立叶系数， 在信号处理和系统分析中更容易使用，余弦傅在信号处理和系统分析中更容易使用，余弦傅 立叶级数在周期信号和幅度谱，相位谱中更加立叶级数在周期信号和幅度谱，相位谱中更加 直观。直观。 .3傅里叶变换和频率域介绍傅里叶变换和频率域介绍 傅里叶变换：非周期函数可以用正弦和傅里叶变换：非周期函数可以用正弦和/ /或余弦乘以加权函或余弦乘以加权函 数的积分来表示。数的积分来表示。 傅里叶反变换：函数特征可以通过反变换来重建，不丢失任傅里叶反变换：函数特征可以通过反变换来重建，不丢失任 何信息。何信息。 一维傅里叶变换及其反变换一维傅里叶变换及其反变换

8、 二维二维DFTDFT及其反变换及其反变换 频率域滤波频率域滤波 空间域滤波和频率域滤波之间的对应关系空间域滤波和频率域滤波之间的对应关系 一维傅里叶变换及其反变换一维傅里叶变换及其反变换 dxexfuF xuj 2 设 x：空间变量（实变量） f(x)：实变量x的连续函数 u：频率变量（实变量） F(u)：频率函数（有实部和虚部） 傅里叶正变换为： 若已知F(u), 则利用傅里叶反变换，可求得f(x) 1 2 jdueuFxf xuj Mxuj M x M exfu / 2 1 0 1 )()(F Mxuj M u euFxf / 2 1 0 )()( 1, 2 , 1 , 0Mu 1, 2

9、 , 1 , 0Mx 离散形式： 例：两个简单一维函数的傅里叶谱例：两个简单一维函数的傅里叶谱 、当曲线下的面积在x域加倍时， 频率谱的高度也加倍。 、当函数的长度加倍时，相同间 隔下频谱中零点的数量也加倍。 dxdyeyxfvuF vyuxj)(2 ),(),( 定义: 若f(x,y)是连续图像函数 正变换: 反变换: dudvevuFyxf vyuxj)(2 ),(),( 变换对: ),(),(vuFyxf 二维二维DFTDFT及其反变换（及其反变换（2DFT2DFT） 定义: 若f(x,y)是离散图像函数（尺寸M*N） 正变换: 反变换: 1, 1 , 0 1, 1 , 0 ),( 1

10、),( 1 0 1 0 )/(2 Nv Mu eyxf MN vuF M x N y NvyMuxj 1, 1 , 0 1, 1 , 0 ),(),( 1 0 1 0 )/(2 Ny Mx evuFyxf M u N v NvyMuxj 一般一般F(u,v)是复函数是复函数,即即: ),( ),(),(),(),( vuj evuFvujIvuRvuF ),(),(),( 22 vuIvuRvuF 幅度谱: ),( ),( ),( 1 vuR vuI tgvu 相位谱: 功率谱： ),(),(),( 22 vuIvuRvuP 6.2.3傅里叶变换和频率域的介绍傅里叶变换和频率域的介绍 移中性

11、n DFT取的区间是0,N-1，在这 个区间内频谱是由两个背靠背的 半周期组成的 n 要显示一个完整的周期， 必须将变换的原点移至 u=N/2点。 图像移中后进行傅里叶变换，则变换后主要能量（低频分图像移中后进行傅里叶变换，则变换后主要能量（低频分 量）集中在频率平面的中心量）集中在频率平面的中心（M/2,N/2）； DFT的的原点，即原点，即F(0,0)被设置在被设置在u=M/2和和v=N/2上；上； 如果是一幅图像，在原点的傅里叶变换如果是一幅图像，在原点的傅里叶变换F(0,0)等于图像的等于图像的 平均灰度级，也称作频率谱的直流成分。平均灰度级，也称作频率谱的直流成分。 1 0 1 0

12、),( 1 )0 ,0( M x N y yxf MN F FT 移中的变换： 能量分布于四角 能量集中于中心 移中FT 原图像f（x，y） 移中FT Matlab显示频谱 I=imread(lena2.png); I=im2double(I); fcoef=fft2(I); spectrum=fftshift(fcoef); temp=log(1+abs(spectrum); subplot(1,2,1); imshow(temp,); subplot(1,2,2); imshow(I); .4直接计算DFT的问题及改进途 径 设x(n)为N点有限长序列，其DFT为 一般来说

13、，x(n)和 都是复数，X(k)也是复数， 因此每计算一个X(k)值，需要N次复数乘法以及 （N-1)次复数加法。而X(k)一共有N个点，所以 完成整个DFT运算总共需要 次复数乘法及 N(N-1)次复数加法。 ），（11.31,.1 ,0,)()( 1 0 NkWnxkX nk N N n nk N w 2 N 4N时时 共需16次乘法，12次加法。 1）直接计算DFT的问题及改进途径 由于一次复数乘法需用四次实数加法；一次复数 加法则需二次实数加法。因此每运算一个X(k)需 要4N次实数乘法及2N+2(N-1)=2(2N-1)次实数 加 法。所以整个DFT运算总共需要 次实数乘法和 次加

14、法。 例如：N=1024时，DFT需要复乘1,048,576次。 所 以，直接计算DFT对实时性很强的信号处理来说， 对计算速度要求是太高了。 2 4N ) 12(2) 12(2NNNN 1）直接计算DFT的问题及改进途径 仔细观察DFT的运算就可以看出，利用系数 的 以下固有特性，就可以减小DFT的运算量： 1、 的周期性 符号 表示取nk除以N所得之余数。 此特性的另一种表达式为： 2、 的对称性 因为 ，于是得到 nk N W nk N W ），（31.3 )( N nknk WW nk N W ），（61 . 3 ) 2 ( nk N nk WW ），（），（51 . 341 . 3

15、)()( nk N nNk N nk N kNn N WWWW 1 2 N W N nk )( 1）直接计算DFT的问题及改进途径 这样，(1)利用这些特性，使DFT运算中有些项可 以合并；(2)利用 的对称性和周期性，可以将 长序列的DFT分解为短序列的DFT。前面已经说 过，DFT的运算量与 成正比，所以N越小越有 利。 快速傅立叶变换就是在这种思路下发展起来的。 下面将详细介绍。 nk N W 2 N 2）按时间抽取的FFT算法 先设序列长度为 ,M为整数。 将 的序列x(n)(n=0,1,.,N-1),先按 n的奇偶分成两组 r=0,1, -1 则将DFT化为 M N 2 M N 2

16、)()12( )()2( 2 1 rxrx rxrx 2N ) 12 . 3( ,)(12()(2()( )()( 1 0 1 2 0 1 2 0 22 N n N r N r rk N k N rk N nk N WrxWWrxWnx nxDFTkX 由于 则上式(3.2-1)可表示成 )22 . 3( , 2 ) 2 ( 2 ) 2 (2 2 N WeeW N j N j N ），（32 . 3)()( ) 12()2()( 2 1 2 0 2 1 2 0 kHWkG WrxWWrxkX k N rk N N r k N rk N N r . 1 2 , 1 , 0 N k 2）按时间抽取

17、的FFT算法 式中 及 分别是 及 的 点DFT 由上式可以看出，一个N点DFT已分解为两个 点DFT。 )(kG)(kH)( 1 rx)( 2 rx 2 N ），（42 .3)2()( 1 2 0 2 N r rk N WrxkG ），（52 .3)12()( 1 2 0 2 N r rk N WrxkH 2 N 2）按时间抽取的FFT算法 但是， 以及 都是 点的序列即r,k满足r,k=0,1, -1。而 却有N点，而用上面的式子计算得到的只是 的前一半项数的结果，要用 来 表达全部的 值，还必须应用系数的周期 性，即 这样可得到 )(),( 21 rxrx)(),(kHkG 2 N 2

18、N )(kX )(kX)(),(kHkG )(kX ) 2 ( 22 N kr N rk N WW )62 . 3( , )()()() 2 ( 1 2 0 1 2 0 2 1 ) 2 ( 2 1 N r N r rk N k N r N kGWrxWrxk N G 2）按时间抽取的FFT算法 同理可得： 上两式说明了后半部分k值( ) 所对应的 分别等于前 半部分k值( )所对应 的 。 再考虑到 的对称性 这样，就可以得到 分为前后两部分得 表达： 前半部分 ），（72 . 3)() 2 (kHk N H 1 2 Nk N )(),( 21 kXkX 1 2 0 N k)(),( 21 k

19、XkX k N Wk N k N N N k N N WWWW 2 ) 2 ( )(kX )82 . 3( , 1 2 ,.1 , 0),()()( : ) 1 2 ,.,1 , 0)( N kkHWkGkX N kkX k N 2）按时间抽取的FFT算法 后半部分 这样，只要求出0到( )区间得所有 值，即可求出0到（N1）区间内的所有 值 ，这就大大节省了运算。 )92 . 3( , 1 2 ,.1 , 0),()( ) 2 () 2 () 2 ( : ) 1,., 2 )( N kkHWkG N kHW N kG N kX N N kkX k N k N 1 2 N )(),(kHkG

20、)(kX 2）按时间抽取的FFT算法 上面的式子可以用下面的蝶形图符号表示： 可以看出，每个蝶形运算，需要一次复数乘法 及两次复数加(减)法。据此，一个N点DFT 分解为两个 点DFT只需要 次复数乘法， 次复数加法，两个 点 DFT共需 次复 )(kG )(kH )()(kHWkG k N )()(kHWkG k N nk N W 1 2 N 4 ) 2 ( 2 2 NN )1 2 ( 2 NN 2 N 2 ) 2 (2 2 2 NN 2）按时间抽取的FFT算法 数乘法和 次复数加法。此外把两个 点 DFT合成N点DFT时，有 个蝶形运算，还需要 次复数乘法及 次复数加法。因而通过这一步 分

21、解后，总共需要 次复数 乘法和 次复数加法，而原始的直接DFT方 法需 要 N*N次复数乘法，N(N-1)次复数加法。 。 )1 2 ( N N 2 N 2 N 2 N N N N M N 2 log 22 NNMN 2 log 2）按时间抽取的FFT算法 例题分析（设N=4) 6.3 6.3 傅里叶变换和频率域滤波的介绍傅里叶变换和频率域滤波的介绍 频率域滤波频率域滤波 一、频率域的基本性质一、频率域的基本性质 图像图像变化平缓的部分变化平缓的部分靠近频率平面的圆心，这个区域为靠近频率平面的圆心，这个区域为 低频区域（如房间中的墙和地板）低频区域（如房间中的墙和地板）； 图像中的图像中的边、

22、噪音、变化陡峻的部分边、噪音、变化陡峻的部分，以放射方向离开，以放射方向离开 频率平面的圆心，这个区域为频率平面的圆心，这个区域为高频区域高频区域。 u v 6.3.1 傅里叶变换和频率域的介绍傅里叶变换和频率域的介绍 边缘、噪音、变化陡峭部分 变化平缓部分 6.3.2傅里叶变换和频率域的介绍傅里叶变换和频率域的介绍 频域中滤波基础 频率域滤波步骤： 用(-1)x+y乘以输入图像来进行中心变换； 由(1)计算图像的DFT，即F(u,v)； 用滤波器函数H(u,v)乘以F(u,v); 计算(3)中结果的反DFT; 得到(4)中结果的实部； 1. 用(-1)x+y乘以(5)中的结果。 ),(),(

23、),(vuFvuHvuG ),( 1 vuG 滤波后的图像 空间域滤波和频率域滤波之间的对应关系空间域滤波和频率域滤波之间的对应关系 的离散卷积为：和),(),(yxhyxf 将图像的模板在图像中逐像素移动，并对每个像素进将图像的模板在图像中逐像素移动，并对每个像素进 行指定数量计算的过程就是卷积过程。行指定数量计算的过程就是卷积过程。 1 0 1 0 ),( ),( 1 ),(),( M m N n nymxhnmf MN yxhyxf 卷积：1、翻转；2、移动；3、乘积； 4、求和 空间域滤波和频率域滤波之间的对应关系空间域滤波和频率域滤波之间的对应关系 ),(),(),(),(vuHvu

24、Fyxhyxf ),(),(),(),(vuHvuFyxhyxf 卷积定理： 6.4频率域滤波器频率域滤波器 空间平滑滤波器 消除或减弱图像中灰度值具有较大较快变化部分的影响，这些部分对应频域中的高频分 量，所以可用频域低通滤波来实现 空间锐化滤波器 消除或减弱图像中灰度值缓慢变化的部分，这些部分对应频域中的低频分量，所以可用 频域高通滤波来实现 频率域滤波模型: ),(),(),(vuFvuHvuG F(u,v)为含有噪声原图像的傅里叶变换 H(u,v)为低通滤波器的传递函数 G(u,v)为经低通滤波后输出图像的傅里叶变换 6.4.1 图像平滑图像平滑低通滤波器低通滤波器 u uF ( )x

25、f ( ) x 00 低通滤波器低通滤波器空域通过用带空域通过用带正系数正系数的模版实现的模版实现低通低通滤波；滤波； 频域越频域越宽宽，滤除的频率成份越滤除的频率成份越少少，空域越空域越窄窄，模板越模板越小小，平平 滑作用越弱滑作用越弱。 频域越频域越窄窄，滤除的频率成份越，滤除的频率成份越多多，空域越，空域越宽宽，模板越，模板越大大，图图 像越模糊像越模糊； 1、理想低通滤波 器 2、高斯低通滤波 器 6.4 .2低通滤波器低通滤波器理想低通滤波器理想低通滤波器(ILPF) (ILPF) 到频率矩形原点的距离 截止频率 ),(),( ),(; 0 ),(; 1 ),( 0 0 0 vuvu

26、D D DvuD DvuD vuH 2/1 22 2/2/),(NvMuvuD 理想低通滤波器作用理想低通滤波器作用 D0半径内的频率分量无损通过，圆外的频率分量会被滤除 若滤除的高频分量中含有大量的边缘信息，会发生图像边缘模糊现象。 标准截止频率：是通过计算半径为r的圆包围的图像功率 的百分比决定的。 1 0 1 0 ),( M u N v T vuPP 半径为r的圆包含a%的功率 u T v PvuPa/ ),(100 在谱中叠加的圆周分别有5，15，30，80，230像素的半径。 这些圆周包围的图像功率的百分比分别为92.0%，94.6%，96.4%，98%，99.5%。 理理 想想 低低 通通 滤滤 波波 器器 5 15 30 80 230 ),(),(),(vuFvuHvuG 频域： 空间域： ),(*),(),(yxfyxhy