图像傅里叶分析

转载自http://blog.csdn.net/songzitea/article/details/16344441

引言

由于图像处理关心的是采样后的数据，我们需要可以处理这些数据的一种傅里叶变换。傅里叶分析可以分析不同滤波器的频率特征。本节，我们将阐述如何通过傅里叶分析认识图像的频率信息。

基本理论

对于一个给定的滤波器，如果使用正弦信号s(x)和脉冲响应为h(x)的滤波器进行卷积后，我们可以得到另一个频率相

同但幅度A 和相位φ 不同的正弦波 ,ο(x) =h(x)*s(x) = Asin(ωx+φ), 如图3.23所示。实际中，傅里叶变换使用的更复杂的符号，即正弦波的复数值形式

在这种情况下，它可以简写为

傅里叶变换是对每个频率的幅度各相位响应的简单罗列，即我们可以理解为：一个恰当地比喻是将傅里叶变换看作一个玻离棱镜。棱镜是可以将光分解为不同颜色的物理仪器，每个成分的颜色由波长/频率来决定。傅里叶变换可以看作是数学上的棱镜，将函数基于频率分解为不同的成分。当我们考虑光时,讨论它的光谱或频率谱。同样, 傅立叶变换使我们能通过频率成分来分析一个函数(此处内容来自于[4]书中)。

即频率为ω的复数正弦波通过滤波器h(x)的响应。

上述公式并没有给出实际计算傅里叶变换的公式。不过它给出了一个方法，即不断用正弦波和滤波器卷积，观察幅度和相位差。在连续域和离散域中都存在傅里叶变换公式分别如下所示：

傅里叶变换的离散形式(3.53)也称为“离散傅里叶变换(DFT)”。注意，虽然k取任何值，式(3.53)都是可以计算的，但是只有当k 属于[-N/2,N/2]时，它才有意义。

现在我们已经定义了傅里叶变换，下面将讨论它的性质和应用。表3.1 列出了一些有用的性质。

我们详细描述这些性质：

• 叠加(superposition)： 几个信号和的傅里叶变换等于每个信号傅里叶变换的和。因此，傅里叶变换是线性算子。
• 平移(shift)： 将一个信号平移后的傅里叶变换等于原始信号的傅里叶变换乘以线性移相。
• 反向(reversal)：反向信号的傅里叶变换等于原信号的傅里叶变换的共轭。
• 卷积(convolution)： 两个信号的卷积的傅里叶变换等于它们的傅里叶变换的乘积。
• 相关(correlation): 两个信号的相关的傅里叶变换等于第一个信号的傅里叶变换与第二个信号的傅里叶变换的共轭的乘积。(注：关于卷积与相关两个专业语术的更详细地请见[4]一书中第46-50页,作者详细地阐述了).
• 乘(multiplication):两个信号的乘积的傅里叶变换等于它们的傅里叶变换的卷积。
• 微分(differentiation): 一个信号的微分的傅里叶变换等于些信号的傅里叶变换与频率的乘积，即：微分使高频部分线性放大。
• 定义域缩放(domain scaling): 将一个信号拉伸后的傅里叶变换等于压扁原信号的傅里叶变换。
• 实值图像(real Images) :一个实数信号的傅里叶变换关于原点变换对称。也可以在变换过程中通过对信号的实部与虚部交替扫描使图像的FFT速度提高一半。
• Parseval定理(Parseval's Thm.): 一个信号能量(信号值的平方和)等于它的傅里叶变换的能量。
  

关注一些常用的滤波器和信号的傅里叶变对，参见表3.2：

在[3]文中提到是：Fourier theory states that any signal, in our case visual images, can be expressed as a sum of a series of sinusoids. In the case of imagery, these are sinusoidal variations in brightness across the image。所以，我们可以将正弦模式在单傅里叶中由三个分量编码：频率(f)、幅值(A)、相位(γ) 这三个value可以描述正弦图像中的所有信息。频率(frequency)：我们可以理解为在空间域上可由亮度调节；幅值(amplitude/magnitude)：正弦函数的幅值用于描述对比度或说是图像中最明和最暗的峰值之间的差(在数学域上，正弦函数是最大峰值A与最低峰值-A之间的差，其中一个负幅表示一个对比逆转，即明暗交换)；相位(phase)：表示相对于原始信号/波形，这个信号/波形的偏移量。

一个傅里叶变换是一系正弦曲线的编码。它们的频率从0开始（相位为0)到Nyquist频率。傅里叶变换同时将图像中所有频率进行编码：一个只包含一个频率f的信号在频谱上横坐标N为f的点处绘制一个单峰值，峰值高度等于对应的幅度amplitude,或正弦曲线信号的高度，如下图所示(注：此以下两张图片均来自于[3])：

上图中DC term直流信号对应于频率为0的点，表示整幅图像的平均亮度，如果直流信号DC=0就表示整幅图像平均亮度的像素点个数=0，可推出 灰度图中，正弦曲线在正负值之间交替变化，但是由于灰度图中没有负值，所以所有的真实图像都有一个正的DC term。

Actually, for mathematical reasons beyond the scope of this tutorial, the Fourier transform also plots a mirror-image of the spatial frequency plot reflected across the origin, with spatial frequency increasing in both directions from the origin. For mathematical reasons beyond the scope of this explanation, these two plots are always mirror-image reflections of each other, with identical peaks at f and -f as shown below.即：因在数学分析原因，经常把傅里叶用mirror-image表示，在原点的两端频率都是增加的方向，同时具有幅值， 如下图所示：

上述都是一维信号，一个二维傅里叶变换是一维傅里叶变换在每一个行扫描线和列扫描线上的傅里叶变换的叠加。傅里叶谱图上的每一个像素点都代表一个频率值，幅值由像素点亮度变码而得。最中心的亮点是指直流分量，傅里叶谱图中越亮的点，对应于灰度图中对比越强烈（对比度越大）的点。由于每一列扫描线上没有变化，所以相应的傅里叶频谱图上行向量为0, 每一行扫描线上有常量，所以有频率幅值。

         

图1(a)  原始图                 图1(b)傅里叶频谱图             图1(c)傅里叶逆变换

          

图2(a)  原始图                 图2(b)傅里叶频谱图             图2(c)傅里叶逆变换

比较图2(b)与图1(b)傅里叶频谱图，不难看出，图2(a)频率比图1(a)频率上面的小，相应的亮点比图1(b)也集中。同时与作者[3]给出图像傅里叶频谱图一致。（说明本人对作者思想基本理解到位，呵呵呵）。

图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。如：大面积的沙漠在图像中是一片灰度变化缓慢的区域，对应的频率值很低；而对于地表属性变换剧烈的边缘区域在图像中是一片灰度变化剧烈的区域，对应的频率值较高。傅立叶变换在实际中有非常明显的物理意义，设f是一个能量有限的模拟信号，则其傅立叶变换就表示f的谱。从纯粹的数学意义上看，傅立叶变换是将一个函数转换为一系列周期函数来处理的。从物理效果看，傅立叶变换是将图像从空间域转换到频率域，其逆变换是将图像从频率域转换到空间域。换句话说，傅立叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，傅立叶逆变换是将图像的频率分布函数变换为灰度分布函数傅立叶变换以前，图像（未压缩的位图）是由对在连续空间（现实空间）上的采样得到一系列点的集合，我们习惯用一个二维矩阵表示空间上各点，则图像可由z=f(x,y)来表示。由于空间是三维的，图像是二维的，因此空间中物体在另一个维度上的关系就由梯度来表示，这样我们可以通过观察图像得知物体在三维空间中的对应关系。

为什么要提梯度？因为实际上对图像进行二维傅立叶变换得到频谱图，就是图像梯度的分布图，当然频谱图上的各点与图像上各点并不存在一一对应的关系，即使在不移频的情况下也是没有。傅立叶频谱图上我们看到的明暗不一的亮点，实际上图像上某一点与邻域点差异的强弱，即梯度的大小，也即该点的频率的大小（可以这么理解，图像中的低频部分指低梯度的点，高频部分相反）。

一般来讲，梯度大则该点的亮度强，否则该点亮度弱。这样通过观察傅立叶变换后的频谱图，也叫功率图，我们首先就可以看出，图像的能量分布，如果频谱图中暗的点数更多，那么实际图像是比较柔和的（因为各点与邻域差异都不大，梯度相对较小），反之，如果频谱图中亮的点数多，那么实际图像一定是尖锐的，边界分明且边界两边像素差异较大的。对频谱移频到原点以后，可以看出图像的频率分布是以原点为圆心，对称分布的。将频谱移频到圆心除了可以清晰地看出图像频率分布以外，还有一个好处，它可以分离出有周期性规律的干扰信号，比如正弦干扰，一副带有正弦干扰，移频到原点的频谱图上可以看出除了中心以外还存在以某一点为中心，对称分布的亮点集合，这个集合就是干扰噪音产生的，这时可以很直观的通过在该位置放置带阻滤波器消除干扰。

           

图3(a)原始图(白色图像)         图3(b)相对应的频谱图                   图3(c)原始图        图3(d)相对应的频谱图

 

参考代码

OpenCV1.x版

1. IplImage* ImgTest::ImgDFT(IplImage*pImg){  
2.     int dft_M, dft_N;  
3.     CvMat* dft_A, tmp;  
4.     double m, M;  
5.   
6.     if( !pImg ||pImg->nChannels ==3){  
7.         printf("Error: Load Gray File in ImgDFT(...).\n");  
8.         exit(EXIT_FAILURE);  
9.     }  
10.   
11.     IplImage * realInput = cvCreateImage( cvGetSize(pImg), IPL_DEPTH_64F, 1);  
12.     IplImage *imaginaryInput = cvCreateImage( cvGetSize(pImg), IPL_DEPTH_64F, 1);  
13.     IplImage *complexInput = cvCreateImage( cvGetSize(pImg), IPL_DEPTH_64F, 2);  
14.   
15.     cvScale(pImg, realInput, 1.0, 0.0);  
16.     cvZero(imaginaryInput);  
17.     cvMerge(realInput, imaginaryInput, NULL, NULL, complexInput);  
18.   
19.     dft_M = cvGetOptimalDFTSize( pImg->height - 1 );  
20.     dft_N = cvGetOptimalDFTSize( pImg->width - 1 );  
21.   
22.     dft_A = cvCreateMat( dft_M, dft_N, CV_64FC2 );  
23.     IplImage * image_Re = cvCreateImage( cvSize(dft_N, dft_M), IPL_DEPTH_64F, 1);  
24.     IplImage * image_Im = cvCreateImage( cvSize(dft_N, dft_M), IPL_DEPTH_64F, 1);  
25.   
26.     // copy A to dft_A and pad dft_A with zeros  
27.     cvGetSubRect( dft_A, &tmp, cvRect(0,0, pImg->width, pImg->height));  
28.     cvCopy( complexInput, &tmp, NULL );  
29.     cvGetSubRect(dft_A,&tmp,cvRect(pImg->width,0,dft_A->cols-pImg->width,pImg->height));  
30.   
31.     // no need to pad bottom part of dft_A with zeros because of  
32.     // use nonzero_rows parameter in cvDFT() call below  
33.     cvDFT( dft_A, dft_A, CV_DXT_FORWARD, complexInput->height );  
34.     //Split Fourier in real and imaginary parts  
35.     cvSplit( dft_A, image_Re, image_Im, 0, 0 );   
36.   
37.     // Mag = sqrt(Re^2 + Im^2)  
38.     cvPow( image_Re, image_Re, 2.0);  
39.     cvPow( image_Im, image_Im, 2.0);  
40.     cvAdd( image_Re, image_Im, image_Re, NULL);  
41.     cvPow( image_Re, image_Re, 0.5 );  
42.   
43.     cvAddS( image_Re, cvScalarAll(1.0), image_Re, NULL ); // 1 + Mag  
44.     cvLog( image_Re, image_Re ); // log(1 + Mag)  
45.   
46.     ShiftDFT( image_Re, image_Re );//此函数的功能是以图像中心为原点，调整傅立叶变换图像的四个象  
47.                                    //限区，即第一和第三象限交换,第二和第四象限交换   
48.     cvMinMaxLoc(image_Re, &m, &M, NULL, NULL, NULL);  
49.     cvScale(image_Re, image_Re, 1.0/(M-m), 1.0*(-m)/(M-m));  
50.     cvShowImage("Magnitude", image_Re);  
51.       
52.     cvWaitKey(0);  
53.     return 0;  
54. }  

注：ShiftDFT( ) 函数可以根据注释思路完成此函数的功能。

python版

1. import sys  
2.   
3. # Rearrange the quadrants of Fourier image so that the origin is at  
4. # the image center  
5. # src & dst arrays of equal size & type  
6.     def ShiftDFT(src_arr, dst_arr ):  
7.    
8.         size = cvGetSize(src_arr)  
9.         dst_size = cvGetSize(dst_arr)  
10.        
11.         if(dst_size.width != size.width or  
12.                 dst_size.height != size.height) :  
13.             cvError(CV_StsUnmatchedSizes, "ShiftDFT", "Source and Destination arrays   
14.                         must have equal sizes", __FILE__, __LINE__ )     
15.        
16.         if(src_arr is dst_arr):  
17.             tmp = cvCreateMat(size.height/2, size.width/2, cvGetElemType(src_arr))  
18.            
19.         cx = size.width/2  
20.         cy = size.height/2 # image center  
21.         q1 = cvGetSubRect( src_arr, cvRect(0,0,cx, cy) )  
22.         q2 = cvGetSubRect( src_arr, cvRect(cx,0,cx,cy) )  
23.         q3 = cvGetSubRect( src_arr, cvRect(cx,cy,cx,cy) )  
24.         q4 = cvGetSubRect( src_arr, cvRect(0,cy,cx,cy) )  
25.         d1 = cvGetSubRect( src_arr, cvRect(0,0,cx,cy) )  
26.         d2 = cvGetSubRect( src_arr, cvRect(cx,0,cx,cy) )  
27.         d3 = cvGetSubRect( src_arr, cvRect(cx,cy,cx,cy) )  
28.         d4 = cvGetSubRect( src_arr, cvRect(0,cy,cx,cy) )  
29.        
30.         if(src_arr is not dst_arr):  
31.             if( not CV_ARE_TYPES_EQ( q1, d1 )):  
32.                 cvError(CV_StsUnmatchedFormats, "ShiftDFT", "Source and Destination   
33.                             arrays must have the same format", __FILE__, __LINE__ )     
34.                
35.             cvCopy(q3, d1)  
36.             cvCopy(q4, d2)  
37.             cvCopy(q1, d3)  
38.             cvCopy(q2, d4)  
39.            
40.         else:  
41.             cvCopy(q3, tmp)  
42.             cvCopy(q1, q3)  
43.             cvCopy(tmp, q1)  
44.             cvCopy(q4, tmp)  
45.             cvCopy(q2, q4)  
46.             cvCopy(tmp, q2)  
47.        
48.     if __name__ == "__main__":  
49.            
50.         im = cvLoadImage( sys.argv[1], CV_LOAD_IMAGE_GRAYSCALE)    
51.         realInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 1)  
52.         imaginaryInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 1)  
53.         complexInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 2)  
54.        
55.         cvScale(im, realInput, 1.0, 0.0)  
56.         cvZero(imaginaryInput)  
57.         cvMerge(realInput, imaginaryInput, None, None, complexInput)  
58.        
59.         dft_M = cvGetOptimalDFTSize( im.height - 1 )  
60.         dft_N = cvGetOptimalDFTSize( im.width - 1 )  
61.         dft_A = cvCreateMat( dft_M, dft_N, CV_64FC2 )  
62.         image_Re = cvCreateImage( cvSize(dft_N, dft_M), IPL_DEPTH_64F, 1)  
63.         image_Im = cvCreateImage( cvSize(dft_N, dft_M), IPL_DEPTH_64F, 1)  
64.    
65.         # copy A to dft_A and pad dft_A with zeros  
66.         tmp = cvGetSubRect( dft_A, cvRect(0,0, im.width, im.height))  
67.         cvCopy( complexInput, tmp, None )  
68.         if(dft_A.width > im.width):  
69.             tmp = cvGetSubRect(dft_A, cvRect(im.width,0,dft_N - im.width, im.height))  
70.             cvZero( tmp )  
71.        
72.         # no need to pad bottom part of dft_A with zeros because of  
73.         # use nonzero_rows parameter in cvDFT() call below  
74.        
75.         cvDFT( dft_A, dft_A, CV_DXT_FORWARD, complexInput.height )  
76.        
77.         cvNamedWindow("win", 0)  
78.         cvNamedWindow("magnitude", 0)  
79.         cvShowImage("win", im)  
80.        
81.         # Split Fourier in real and imaginary parts  
82.         cvSplit( dft_A, image_Re, image_Im, None, None )  
83.        
84.         # Compute the magnitude of the spectrum Mag = sqrt(Re^2 + Im^2)  
85.         cvPow( image_Re, image_Re, 2.0)  
86.         cvPow( image_Im, image_Im, 2.0)  
87.         cvAdd( image_Re, image_Im, image_Re, None)  
88.         cvPow( image_Re, image_Re, 0.5 )  
89.        
90.         # Compute log(1 + Mag)  
91.         cvAddS( image_Re, cvScalarAll(1.0), image_Re, None ) # 1 + Mag  
92.         cvLog( image_Re, image_Re ) # log(1 + Mag)  
93.            
94.         # Rearrange the quadrants of Fourier image so that the origin is at  
95.         # the image center  
96.         cvShiftDFT( image_Re, image_Re )  
97.        
98.         min, max, pt1, pt2 = cvMinMaxLoc(image_Re)  
99.         cvScale(image_Re, image_Re, 1.0/(max-min), 1.0*(-min)/(max-min))  
100.         cvShowImage("magnitude", image_Re)  
101.        
102.         cvWaitKey(0)  

参考资料

[1] Fourier Transform from Wikipedia, the free encyclopedia.

[2] EE261: The Fourier Transform and its Applications by Professor Brad G. Osgood(Stanford Univeristy).

[3] Steven Lehar,"An Intuitive Explanation of Fourier Theory".

[4] Rafael C. Gonzalez, Rechard E. Woods Steven L. Eddins, "Digital Image Processing Using MatLab (2nd Edition)", Gatesmark Publishing 2009,

[5] Tom Irvine, "Shock and Vibration Signal Analysis", September 12, 2005.

[6] Richard Szeliski,"Computer Vision:Algorithms and Applications", Microsoft Research, 2010.

[7] Steven W. Smith,"The Scientist and Engineer's Guide to Digital Signal Processing",http://www.dspguide.com/pdfbook.htm.

关于Image Engineering & Computer Vision的更多讨论与交流，敬请关注本博和新浪微博songzi_tea.