2019-09-29 09:44:19 weixin_43392489 阅读数 135
  • 携手STM32CubeMX玩转STM32

    本课程教大家如何利用STM32CubeMX玩转STM32(STM32CubeMX支持的所有型号MCU都适用于本课程)。课程内容分为基础内容和扩展内容,例如:讲解串口时会扩展开讲Xmodem通信协议,讲解ADC/DAC时会扩展讲傅里叶计算,讲解完FLASH操作会扩展将bootloader的编写,讲解完M3的bootloader编写会扩展讲解M0的bootloader...... 内容绝对实在,对于学习以及工作都会有很大的帮助。最终的目的就是让大家学会快速开发STM32并收获与STM32有关的实用技术知识。

    795 人正在学习 去看看 李凯龙

 图像的傅里叶频谱特性分析

图像傅里叶频谱关于(/,/)的对称性

图像傅里叶频谱特性及其频谱图

傅里叶变换在图像处理中的应用

 

 

2017-03-25 16:47:59 yimingsilence 阅读数 1942
  • 携手STM32CubeMX玩转STM32

    本课程教大家如何利用STM32CubeMX玩转STM32(STM32CubeMX支持的所有型号MCU都适用于本课程)。课程内容分为基础内容和扩展内容,例如:讲解串口时会扩展开讲Xmodem通信协议,讲解ADC/DAC时会扩展讲傅里叶计算,讲解完FLASH操作会扩展将bootloader的编写,讲解完M3的bootloader编写会扩展讲解M0的bootloader...... 内容绝对实在,对于学习以及工作都会有很大的帮助。最终的目的就是让大家学会快速开发STM32并收获与STM32有关的实用技术知识。

    795 人正在学习 去看看 李凯龙

1.图像的傅里叶频谱的意义

之前的博文其实已经归纳过这方面的内容了。我们常用的图像平滑处理,其实就是一个低通滤波,一定程度上去除高频信号,可以使得图像变得柔和(也就是平滑)。但是,在去除周期性噪声时候,空间域内的滤波(卷积)就不是那么好操作了。所以,这里时候,无论是理解起来方便,还是其他原因,都需要在频域内进行滤波。 
详细的叙述还是在下面的博文里面啦!!!! 

2. 傅里叶频谱的计算

这部分的内容,主要就是使用OpenCV自带的函数 
void cvDFT( const CvArr* src, CvArr* dst, int flags, int nonzero_rows=0 ) 
去求取图像的傅里叶变换。这里,输出结果CvArr* dst由两个通道组成,分别代表了实部与虚部。我们再根据如下算式,就可以得到傅里叶频谱了。 
|F(u,v)|=R2(u,v)+F2(u,v)2 
我自己也参考了很多人的代码,然后实现的代码如下。

IplImage* fft2(IplImage* image_input)
{
    int dftWidth  = getOptimalDFTSize(image_input->width);
    int dftHeight = getOptimalDFTSize(image_input->height);

    //cout<< " Width" <<  image_input->width << "    " <<  dftWidth  << "\n";
    //cout<< "Height" << image_input->height << "    " <<  dftHeight << "\n";

    IplImage* image_padded = cvCreateImage(cvSize(dftWidth,dftHeight),
                                           IPL_DEPTH_8U,
                                           1);
    cvCopyMakeBorder( image_input, image_padded, cvPoint(0,0), IPL_BORDER_CONSTANT,cvScalarAll(0)); 

    IplImage *image_Re =0 , *image_Im = 0, *image_Fourier = 0; 

    image_Re = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
    image_Im = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
    image_Fourier = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,2);

    //image_Re <--- image_padded 
    cvConvertScale(image_padded,image_Re);   
    //image_Im <--- 0
    cvZero(image_Im);                 
    //image_Fourier[0] <--- image_Re
    //image_Fourier[1] <--- image_Im
    cvMerge(image_Re,image_Im,0,0,image_Fourier); 

    cvDFT(image_Fourier,image_Fourier,CV_DXT_FORWARD);

    //image_Fourier[0] ---> image_Re
    //image_Fourier[1] ---> image_Im
    cvSplit(image_Fourier,image_Re,image_Im,0,0);

    //Mag = sqrt(Re^2 + Im^2)
    cvPow(image_Re,image_Re,2.0);
    cvPow(image_Im,image_Im,2.0);
    cvAdd(image_Re,image_Im,image_Re);
    cvPow(image_Re,image_Re,0.5);

    // log (1 + Mag)
    cvAddS(image_Re,cvScalar(1),image_Re ); 
    cvLog (image_Re,image_Re); 

    //  |-----|-----|           |-----|-----|   
    //  |  1  |  3  |           |  4  |  2  |
    //  |-----|-----|   --->    |-----|-----|
    //  |  2  |  4  |           |  3  |  1  |
    //  |-----|-----|           |-----|-----|

    IplImage *Fourier = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
    cvZero(image_Fourier);

    int cx = image_Re->width/2;
    int cy = image_Re->height/2;

    cvSetImageROI(image_Re,cvRect( 0, 0,cx,cy));  // 1 
    cvSetImageROI( Fourier,cvRect(cx,cy,cx,cy));  // 4 
    cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);

    cvSetImageROI(image_Re,cvRect(cx,cy,cx,cy));  // 4 
    cvSetImageROI( Fourier,cvRect( 0, 0,cx,cy));  // 1 
    cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);

    cvSetImageROI(image_Re,cvRect(cx, 0,cx,cy));  // 3 
    cvSetImageROI( Fourier,cvRect( 0,cy,cx,cy));  // 2 
    cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);

    cvSetImageROI(image_Re,cvRect( 0,cy,cx,cy));  // 2 
    cvSetImageROI( Fourier,cvRect(cx, 0,cx,cy));  // 3 
    cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);

    cvResetImageROI(image_Re);
    cvResetImageROI( Fourier);

    cvNormalize(Fourier,Fourier,1,0,CV_C,NULL);

    return(Fourier);
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67
  • 68
  • 69
  • 70
  • 71
  • 72
  • 73
  • 74
  • 75
  • 76
  • 77
  • 78
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67
  • 68
  • 69
  • 70
  • 71
  • 72
  • 73
  • 74
  • 75
  • 76
  • 77
  • 78

从这里开始,还是简单的分析一下代码吧。

int dftWidth  = getOptimalDFTSize(image_input->width);
int dftHeight = getOptimalDFTSize(image_input->height);
IplImage* image_padded = cvCreateImage(cvSize(dftWidth,dftHeight),
                                       IPL_DEPTH_8U,
                                       1);
cvCopyMakeBorder( image_input, image_padded, cvPoint(0,0), IPL_BORDER_CONSTANT,cvScalarAll(0)); 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6

这里参考了文献[2]中的说法,在尺寸数为2,3,5的倍数的场合,计算的速度是最快的。所以使用函数getOptimalDFTSize()来寻找最匹配的尺寸,然后再同伙cvCopyMakeBorder()进行多余部分的填充,这里选的配置是将图放在从点(0,0)开始的位置,其余不足的地方,用0进行填充。

IplImage *image_Re =0 , *image_Im = 0, *image_Fourier = 0; 

image_Re = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
image_Im = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
image_Fourier = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,2);

//image_Re <--- image_padded 
cvConvertScale(image_padded,image_Re);   
//image_Im <--- 0
cvZero(image_Im);                 
//image_Fourier[0] <--- image_Re
//image_Fourier[1] <--- image_Im
cvMerge(image_Re,image_Im,0,0,image_Fourier); 

cvDFT(image_Fourier,image_Fourier,CV_DXT_FORWARD);

//image_Fourier[0] ---> image_Re
//image_Fourier[1] ---> image_Im
cvSplit(image_Fourier,image_Re,image_Im,0,0);
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19

其实这里的很好理解的,将填充到最适尺寸的图像赋值给image_Re,将image_Im赋值为0。让后将这两层图复制到image_Fourier的两个通道里,然后使用函数cvDFT()进行傅里叶变换。得到结果还是存在于image_Fourier的两个通道里,分别代表实部与虚部,然后通过cvSplit()将其抽出到image_Re与image_Im里。

//Mag = sqrt(Re^2 + Im^2)
cvPow(image_Re,image_Re,2.0);
cvPow(image_Im,image_Im,2.0);
cvAdd(image_Re,image_Im,image_Re);
cvPow(image_Re,image_Re,0.5);

// log (1 + Mag)
cvAddS(image_Re,cvScalar(1),image_Re ); 
cvLog (image_Re,image_Re); 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9

以上代码,实现了以下计算。 
|F(u,v)|=R2(u,v)+F2(u,v)2 
还有就是进行了一个对数变换,这个也没的说,看傅里叶频谱的标配操作。

//  |-----|-----|           |-----|-----|   
//  |  1  |  3  |           |  4  |  2  |
//  |-----|-----|   --->    |-----|-----|
//  |  2  |  4  |           |  3  |  1  |
//  |-----|-----|           |-----|-----|

IplImage *Fourier = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
cvZero(image_Fourier);

int cx = image_Re->width/2;
int cy = image_Re->height/2;

cvSetImageROI(image_Re,cvRect( 0, 0,cx,cy));  // 1 
cvSetImageROI( Fourier,cvRect(cx,cy,cx,cy));  // 4 
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);

cvSetImageROI(image_Re,cvRect(cx,cy,cx,cy));  // 4 
cvSetImageROI( Fourier,cvRect( 0, 0,cx,cy));  // 1 
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);

cvSetImageROI(image_Re,cvRect(cx, 0,cx,cy));  // 3 
cvSetImageROI( Fourier,cvRect( 0,cy,cx,cy));  // 2 
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);

cvSetImageROI(image_Re,cvRect( 0,cy,cx,cy));  // 2 
cvSetImageROI( Fourier,cvRect(cx, 0,cx,cy));  // 3 
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);

cvResetImageROI(image_Re);
cvResetImageROI( Fourier);

cvNormalize(Fourier,Fourier,1,0,CV_C,NULL);

return(Fourier);
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34

其实重头戏在这里,这里需要一个交换操作。至于为何所求得的傅里叶频谱为什么需要交换的原因是,这个代码求得的结果其实是范围[0,2π]内的傅里叶变换。而我们所需要的是[π,π]内的结果。详细的原因,还是在我以前的博客里有说明。 
[数字图像处理]频域滤波(1)–基础与低通滤波器

这里,我使用了ROI操作与cvAddWeighted()函数进行了实现。其运行的结果如下所示。 
这里写图片描述
恩,基本可以看出来,直流分量也被我移动到了中心,以上代码实现了傅里叶频谱的计算与显示。

3. 不用交换操作的代码

使用MATLAB去求取尺寸为M×N图像f的傅里叶频谱时候,通常会用fft2(f,2*M,2*N)。使用此函数求得的福利叶变换,其实还是[0,2π]范围内的傅里叶变换。要想使得傅里叶频谱的DC分量位于中心的话,其实还是需要一些别的操作的。冈萨雷斯的《数字图像处理》一书中,对于这个问题,是利用了傅里叶变换的评议特性,即 
f(x,y)=f(x,y)×(1)x+y 
然后再对函数f(x,y)进行福利叶变换,所得到结果,就是所需要的是[π,π]内傅里叶变换的结果。具体的原理与推导,还是参看冈萨雷斯的《数字图像处理》英文原本的p258页左右,中文译本的p148左右。当然,嫌麻烦的话,还可以看我的博文,博文中也简明推导了平移特性。 
[数字图像处理]频域滤波(1)–基础与低通滤波器 
为此,实现的代码变为了如下形式。

IplImage* fft2_New(IplImage* image_input)
{
    int dftWidth  = getOptimalDFTSize(image_input->width);
    int dftHeight = getOptimalDFTSize(image_input->height);

    cout<< " Width" <<  image_input->width << "    " <<  dftWidth  << "\n";
    cout<< "Height" << image_input->height << "    " <<  dftHeight << "\n";

    IplImage* image_padded = cvCreateImage(cvSize(dftWidth,dftHeight),
                                           IPL_DEPTH_8U,
                                           1);
    cvCopyMakeBorder( image_input, image_padded, cvPoint(0,0), IPL_BORDER_CONSTANT,cvScalarAll(0)); 

    IplImage *image_Re =0 , *image_Im = 0, *image_Fourier = 0; 

    image_Re = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
    image_Im = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
    image_Fourier = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,2);

    //image_Re = image_padded .* (-1)^(x+y);  
    double pixel;
    for(int y=0;y<image_padded->height;y++)
    {
        for(int x=0;x<image_padded->width;x++)
        {
            pixel = cvGetReal2D(image_padded,x,y);
            pixel = ((x+y)%2 == 0)?(pixel):((-1)*pixel);
            cvSetReal2D(image_Re,x,y,pixel);
        }
    }

    //image_Im <--- 0
    cvZero(image_Im);                 
    //image_Fourier[0] <--- image_Re
    //image_Fourier[1] <--- image_Im
    cvMerge(image_Re,image_Im,0,0,image_Fourier); 

    cvDFT(image_Fourier,image_Fourier,CV_DXT_FORWARD);

    //image_Fourier[0] ---> image_Re
    //image_Fourier[1] ---> image_Im
    cvSplit(image_Fourier,image_Re,image_Im,0,0);

    //Mag = sqrt(Re^2 + Im^2)
    cvPow(image_Re,image_Re,2.0);
    cvPow(image_Im,image_Im,2.0);
    cvAdd(image_Re,image_Im,image_Re);
    cvPow(image_Re,image_Re,0.5);

    // log (1 + Mag)
    cvAddS(image_Re,cvScalar(1),image_Re ); 
    cvLog (image_Re,image_Re); 

    cvNormalize(image_Re,image_Re,1,0,CV_C,NULL);

    return(image_Re);
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57

在这里,由于考虑到计算的原因,我将f(x,y)=f(x,y)×(1)x+y这个计算,用下面代码去进行了实现。

for(int y=0;y<image_padded->height;y++)
{
    for(int x=0;x<image_padded->width;x++)
    {
        pixel = cvGetReal2D(image_padded,x,y);
        pixel = ((x+y)%2 == 0)?(pixel):((-1)*pixel);
        cvSetReal2D(image_Re,x,y,pixel);
    }
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9

其实也就相当于,x+y是偶数的话f(x,y)不变,x+y是奇数的话f(x,y)变为负。用这样一个简单的判断去实现。其运行的结果,如下所示。

原图的傅里叶频谱 
原图的结果
使用11×11高斯滤波器后的傅里叶频谱 
这里写图片描述

从实验结果看来,可以看出以下两点

  1. 对于原图而言两个函数的结果基本一致,两个函数都得到了正确的结果。
  2. 使用平滑处理后,频谱的高频成分明显变小,对于空间域的图像而言,图像变得模糊了。


2018-06-07 14:25:09 daduzimama 阅读数 1188
  • 携手STM32CubeMX玩转STM32

    本课程教大家如何利用STM32CubeMX玩转STM32(STM32CubeMX支持的所有型号MCU都适用于本课程)。课程内容分为基础内容和扩展内容,例如:讲解串口时会扩展开讲Xmodem通信协议,讲解ADC/DAC时会扩展讲傅里叶计算,讲解完FLASH操作会扩展将bootloader的编写,讲解完M3的bootloader编写会扩展讲解M0的bootloader...... 内容绝对实在,对于学习以及工作都会有很大的帮助。最终的目的就是让大家学会快速开发STM32并收获与STM32有关的实用技术知识。

    795 人正在学习 去看看 李凯龙

                     图像的傅里叶变换在频谱居中上的误区

 

   如果你用MATLAB去计算图像的傅里叶变换那么你一定会用到FFTSHIFT这一函数为了保证计算后的频谱能够居中化。如下图。

                               

Matlab代码:

clear all;
close all;
I = im2double(imread('cameraman.tif'));
imshowpair(I,log(abs(fftshift(fft2(I)))+1),'montage')

 

     但是你如果去看matlab的介绍或者在计算FFT之后调用FFTSHIFT函数,那么经过该函数处理过的图像会变成如下形式。如下图。一定要注意:fftshift这一函数的调用是放在傅里叶变换之后的,这和频谱居中的真实计算原理(即,利用傅里叶变换的属性)根本就不一样,可以说是一种投机取巧的方法,但也不失为一个好的方法。

 

Matlab代码:

Ishift = fftshift(I);
imshow(Ishift,[])

并且MATLAB还给出了如下图所示的示意图,即交换四个象限。注意:所以MATLAB的这种做法是一种投机取巧的做法。

而实际情况不是这样的,根据傅里叶变换的属性或者说是傅里叶变换对。

     要想对频域的信号进行移位,时域的信号需要乘以一个对应的指数函数。这条属性叫做频域移位特性。图像信号在进行变换之前需要对每一个像素乘以(-1)^x+y次方。这是根据欧拉公式计算出来的(注意:e^jπ = -1)。

Matlab代码:

Ie = zeros(size(I));
for i = 1:size(I,1)
    for j = 1:size(I,2)
        Ie(i,j) = (-1)^(i+j).*I(i,j);
    end
end
imshowpair(I,Ie,'montage')

计算出来的图像如下,左边是原始图像,右图是乘以(-1)为底的幂函数后的图像。

     简单点说就是一正一负,每隔一个像素乘以一个负号,得到的图像就好像黑白间隔的棋盘一样。如下图放大后的图像细节所示。而不是那种四个象限互换后的结果。

      最后我用空间域图像乘以指数函数即(-1)^(x+y)的方式对图像的频谱进行居中化的操作,能够得到和MATLAB自带的函数fftshift一模一样的结果。如下图。

      Wherefore,不论是一维信号的频谱还是二维信号的频谱,一定要记得真正的频谱居中不只是简单的fftshift或是四个象限的对调而是在进行变换之前,用一个指数函数乘以原始信号。而这个指数函数根据著名的欧拉公式e^jπ+1=0,最终变成了以(-1)为底的幂函数。

 

(全文完)

 

鸣谢:

【1】Digital image processing, Third Edition.  冈萨雷斯

【2】signals and systems,Second Edition. 奥本海姆

【3】Matlab 2017a

 

谢谢收看!

再见!

 

《圣经》箴言 14章12节   有一条路,人以为正,至终成为死亡之路。

                                                                                                 *配图与本文无关*

 

2015-05-21 10:30:12 thnh169 阅读数 10559
  • 携手STM32CubeMX玩转STM32

    本课程教大家如何利用STM32CubeMX玩转STM32(STM32CubeMX支持的所有型号MCU都适用于本课程)。课程内容分为基础内容和扩展内容,例如:讲解串口时会扩展开讲Xmodem通信协议,讲解ADC/DAC时会扩展讲傅里叶计算,讲解完FLASH操作会扩展将bootloader的编写,讲解完M3的bootloader编写会扩展讲解M0的bootloader...... 内容绝对实在,对于学习以及工作都会有很大的帮助。最终的目的就是让大家学会快速开发STM32并收获与STM32有关的实用技术知识。

    795 人正在学习 去看看 李凯龙

1.图像的傅里叶频谱的意义

之前的博文其实已经归纳过这方面的内容了。我们常用的图像平滑处理,其实就是一个低通滤波,一定程度上去除高频信号,可以使得图像变得柔和(也就是平滑)。但是,在去除周期性噪声时候,空间域内的滤波(卷积)就不是那么好操作了。所以,这里时候,无论是理解起来方便,还是其他原因,都需要在频域内进行滤波。
详细的叙述还是在下面的博文里面啦!!!!
[数字图像处理]频域滤波(1)–基础与低通滤波器
[数字图像处理]频域滤波(2)–高通滤波器,带阻滤波器与陷波滤波器

2. 傅里叶频谱的计算

这部分的内容,主要就是使用openCV自带的函数
void cvDFT( const CvArr* src, CvArr* dst, int flags, int nonzero_rows=0 )
去求取图像的傅里叶变换。这里,输出结果CvArr* dst由两个通道组成,分别代表了实部与虚部。我们再根据如下算式,就可以得到傅里叶频谱了。
|F(u,v)|=R2(u,v)+F2(u,v)2
我自己也参考了很多人的代码,然后实现的代码如下。

IplImage* fft2(IplImage* image_input)
{
    int dftWidth  = getOptimalDFTSize(image_input->width);
    int dftHeight = getOptimalDFTSize(image_input->height);

    //cout<< " Width" <<  image_input->width << "    " <<  dftWidth  << "\n";
    //cout<< "Height" << image_input->height << "    " <<  dftHeight << "\n";

    IplImage* image_padded = cvCreateImage(cvSize(dftWidth,dftHeight),
                                           IPL_DEPTH_8U,
                                           1);
    cvCopyMakeBorder( image_input, image_padded, cvPoint(0,0), IPL_BORDER_CONSTANT,cvScalarAll(0)); 

    IplImage *image_Re =0 , *image_Im = 0, *image_Fourier = 0; 

    image_Re = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
    image_Im = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
    image_Fourier = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,2);

    //image_Re <--- image_padded 
    cvConvertScale(image_padded,image_Re);   
    //image_Im <--- 0
    cvZero(image_Im);                 
    //image_Fourier[0] <--- image_Re
    //image_Fourier[1] <--- image_Im
    cvMerge(image_Re,image_Im,0,0,image_Fourier); 

    cvDFT(image_Fourier,image_Fourier,CV_DXT_FORWARD);

    //image_Fourier[0] ---> image_Re
    //image_Fourier[1] ---> image_Im
    cvSplit(image_Fourier,image_Re,image_Im,0,0);

    //Mag = sqrt(Re^2 + Im^2)
    cvPow(image_Re,image_Re,2.0);
    cvPow(image_Im,image_Im,2.0);
    cvAdd(image_Re,image_Im,image_Re);
    cvPow(image_Re,image_Re,0.5);

    // log (1 + Mag)
    cvAddS(image_Re,cvScalar(1),image_Re ); 
    cvLog (image_Re,image_Re); 

    //  |-----|-----|           |-----|-----|   
    //  |  1  |  3  |           |  4  |  2  |
    //  |-----|-----|   --->    |-----|-----|
    //  |  2  |  4  |           |  3  |  1  |
    //  |-----|-----|           |-----|-----|

    IplImage *Fourier = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
    cvZero(image_Fourier);

    int cx = image_Re->width/2;
    int cy = image_Re->height/2;

    cvSetImageROI(image_Re,cvRect( 0, 0,cx,cy));  // 1 
    cvSetImageROI( Fourier,cvRect(cx,cy,cx,cy));  // 4 
    cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);

    cvSetImageROI(image_Re,cvRect(cx,cy,cx,cy));  // 4 
    cvSetImageROI( Fourier,cvRect( 0, 0,cx,cy));  // 1 
    cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);

    cvSetImageROI(image_Re,cvRect(cx, 0,cx,cy));  // 3 
    cvSetImageROI( Fourier,cvRect( 0,cy,cx,cy));  // 2 
    cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);

    cvSetImageROI(image_Re,cvRect( 0,cy,cx,cy));  // 2 
    cvSetImageROI( Fourier,cvRect(cx, 0,cx,cy));  // 3 
    cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);

    cvResetImageROI(image_Re);
    cvResetImageROI( Fourier);

    cvNormalize(Fourier,Fourier,1,0,CV_C,NULL);

    return(Fourier);
}

从这里开始,还是简单的分析一下代码吧。

int dftWidth  = getOptimalDFTSize(image_input->width);
int dftHeight = getOptimalDFTSize(image_input->height);
IplImage* image_padded = cvCreateImage(cvSize(dftWidth,dftHeight),
                                       IPL_DEPTH_8U,
                                       1);
cvCopyMakeBorder( image_input, image_padded, cvPoint(0,0), IPL_BORDER_CONSTANT,cvScalarAll(0)); 

这里参考了文献[2]中的说法,在尺寸数为2,3,5的倍数的场合,计算的速度是最快的。所以使用函数getOptimalDFTSize()来寻找最匹配的尺寸,然后再同伙cvCopyMakeBorder()进行多余部分的填充,这里选的配置是将图放在从点(0,0)开始的位置,其余不足的地方,用0进行填充。

IplImage *image_Re =0 , *image_Im = 0, *image_Fourier = 0; 

image_Re = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
image_Im = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
image_Fourier = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,2);

//image_Re <--- image_padded 
cvConvertScale(image_padded,image_Re);   
//image_Im <--- 0
cvZero(image_Im);                 
//image_Fourier[0] <--- image_Re
//image_Fourier[1] <--- image_Im
cvMerge(image_Re,image_Im,0,0,image_Fourier); 

cvDFT(image_Fourier,image_Fourier,CV_DXT_FORWARD);

//image_Fourier[0] ---> image_Re
//image_Fourier[1] ---> image_Im
cvSplit(image_Fourier,image_Re,image_Im,0,0);

其实这里的很好理解的,将填充到最适尺寸的图像赋值给image_Re,将image_Im赋值为0。让后将这两层图复制到image_Fourier的两个通道里,然后使用函数cvDFT()进行傅里叶变换。得到结果还是存在于image_Fourier的两个通道里,分别代表实部与虚部,然后通过cvSplit()将其抽出到image_Re与image_Im里。

//Mag = sqrt(Re^2 + Im^2)
cvPow(image_Re,image_Re,2.0);
cvPow(image_Im,image_Im,2.0);
cvAdd(image_Re,image_Im,image_Re);
cvPow(image_Re,image_Re,0.5);

// log (1 + Mag)
cvAddS(image_Re,cvScalar(1),image_Re ); 
cvLog (image_Re,image_Re); 

以上代码,实现了以下计算。
|F(u,v)|=R2(u,v)+F2(u,v)2
还有就是进行了一个对数变换,这个也没的说,看傅里叶频谱的标配操作。

//  |-----|-----|           |-----|-----|   
//  |  1  |  3  |           |  4  |  2  |
//  |-----|-----|   --->    |-----|-----|
//  |  2  |  4  |           |  3  |  1  |
//  |-----|-----|           |-----|-----|

IplImage *Fourier = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
cvZero(image_Fourier);

int cx = image_Re->width/2;
int cy = image_Re->height/2;

cvSetImageROI(image_Re,cvRect( 0, 0,cx,cy));  // 1 
cvSetImageROI( Fourier,cvRect(cx,cy,cx,cy));  // 4 
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);

cvSetImageROI(image_Re,cvRect(cx,cy,cx,cy));  // 4 
cvSetImageROI( Fourier,cvRect( 0, 0,cx,cy));  // 1 
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);

cvSetImageROI(image_Re,cvRect(cx, 0,cx,cy));  // 3 
cvSetImageROI( Fourier,cvRect( 0,cy,cx,cy));  // 2 
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);

cvSetImageROI(image_Re,cvRect( 0,cy,cx,cy));  // 2 
cvSetImageROI( Fourier,cvRect(cx, 0,cx,cy));  // 3 
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);

cvResetImageROI(image_Re);
cvResetImageROI( Fourier);

cvNormalize(Fourier,Fourier,1,0,CV_C,NULL);

return(Fourier);

其实重头戏在这里,这里需要一个交换操作。至于为何所求得的傅里叶频谱为什么需要交换的原因是,这个代码求得的结果其实是范围[0,2π]内的傅里叶变换。而我们所需要的是[π,π]内的结果。详细的原因,还是在我以前的博客里有说明。
[数字图像处理]频域滤波(1)–基础与低通滤波器

这里,我使用了ROI操作与cvAddWeighted()函数进行了实现。其运行的结果如下所示。
这里写图片描述
恩,基本可以看出来,直流分量也被我移动到了中心,以上代码实现了傅里叶频谱的计算与显示。

3. 不用交换操作的代码

使用MATLAB去求取尺寸为M×N图像f的傅里叶频谱时候,通常会用fft2(f,2*M,2*N)。使用此函数求得的福利叶变换,其实还是[0,2π]范围内的傅里叶变换。要想使得傅里叶频谱的DC分量位于中心的话,其实还是需要一些别的操作的。冈萨雷斯的《数字图像处理》一书中,对于这个问题,是利用了傅里叶变换的评议特性,即
f(x,y)=f(x,y)×(1)x+y
然后再对函数f(x,y)进行福利叶变换,所得到结果,就是所需要的是[π,π]内傅里叶变换的结果。具体的原理与推导,还是参看冈萨雷斯的《数字图像处理》英文原本的p258页左右,中文译本的p148左右。当然,嫌麻烦的话,还可以看我的博文,博文中也简明推导了平移特性。
[数字图像处理]频域滤波(1)–基础与低通滤波器
为此,实现的代码变为了如下形式。

IplImage* fft2_New(IplImage* image_input)
{
    int dftWidth  = getOptimalDFTSize(image_input->width);
    int dftHeight = getOptimalDFTSize(image_input->height);

    cout<< " Width" <<  image_input->width << "    " <<  dftWidth  << "\n";
    cout<< "Height" << image_input->height << "    " <<  dftHeight << "\n";

    IplImage* image_padded = cvCreateImage(cvSize(dftWidth,dftHeight),
                                           IPL_DEPTH_8U,
                                           1);
    cvCopyMakeBorder( image_input, image_padded, cvPoint(0,0), IPL_BORDER_CONSTANT,cvScalarAll(0)); 

    IplImage *image_Re =0 , *image_Im = 0, *image_Fourier = 0; 

    image_Re = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
    image_Im = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
    image_Fourier = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,2);

    //image_Re = image_padded .* (-1)^(x+y);  
    double pixel;
    for(int y=0;y<image_padded->height;y++)
    {
        for(int x=0;x<image_padded->width;x++)
        {
            pixel = cvGetReal2D(image_padded,x,y);
            pixel = ((x+y)%2 == 0)?(pixel):((-1)*pixel);
            cvSetReal2D(image_Re,x,y,pixel);
        }
    }

    //image_Im <--- 0
    cvZero(image_Im);                 
    //image_Fourier[0] <--- image_Re
    //image_Fourier[1] <--- image_Im
    cvMerge(image_Re,image_Im,0,0,image_Fourier); 

    cvDFT(image_Fourier,image_Fourier,CV_DXT_FORWARD);

    //image_Fourier[0] ---> image_Re
    //image_Fourier[1] ---> image_Im
    cvSplit(image_Fourier,image_Re,image_Im,0,0);

    //Mag = sqrt(Re^2 + Im^2)
    cvPow(image_Re,image_Re,2.0);
    cvPow(image_Im,image_Im,2.0);
    cvAdd(image_Re,image_Im,image_Re);
    cvPow(image_Re,image_Re,0.5);

    // log (1 + Mag)
    cvAddS(image_Re,cvScalar(1),image_Re ); 
    cvLog (image_Re,image_Re); 

    cvNormalize(image_Re,image_Re,1,0,CV_C,NULL);

    return(image_Re);
}

在这里,由于考虑到计算的原因,我将f(x,y)=f(x,y)×(1)x+y这个计算,用下面代码去进行了实现。

for(int y=0;y<image_padded->height;y++)
{
    for(int x=0;x<image_padded->width;x++)
    {
        pixel = cvGetReal2D(image_padded,x,y);
        pixel = ((x+y)%2 == 0)?(pixel):((-1)*pixel);
        cvSetReal2D(image_Re,x,y,pixel);
    }
}

其实也就相当于,x+y是偶数的话f(x,y)不变,x+y是奇数的话f(x,y)变为负。用这样一个简单的判断去实现。其运行的结果,如下所示。

原图的傅里叶频谱
原图的结果
使用11×11高斯滤波器后的傅里叶频谱
这里写图片描述

从实验结果看来,可以看出以下两点

  1. 对于原图而言两个函数的结果基本一致,两个函数都得到了正确的结果。
  2. 使用平滑处理后,频谱的高频成分明显变小,对于空间域的图像而言,图像变得模糊了。

原文发于博客:http://blog.csdn.net/thnh169/

参考文献

[1]opencv 中 傅里叶变换 FFT :http://blog.csdn.net/abcjennifer/article/details/7359952
[2]OpenCV实现基于傅里叶变换的旋转文本校正 : http://johnhany.net/2013/11/dft-based-text-rotation-correction/#imageclose-380
[3]学习OpenCV范例(八)——离散傅立叶变换 : http://blog.csdn.net/chenjiazhou12/article/details/21240647

=============更新日志===================

2015 - 5 - 19 初版
2015 - 5 - 21 修改了某些叙述

2019-07-06 10:42:52 qq_33208851 阅读数 1847
  • 携手STM32CubeMX玩转STM32

    本课程教大家如何利用STM32CubeMX玩转STM32(STM32CubeMX支持的所有型号MCU都适用于本课程)。课程内容分为基础内容和扩展内容,例如:讲解串口时会扩展开讲Xmodem通信协议,讲解ADC/DAC时会扩展讲傅里叶计算,讲解完FLASH操作会扩展将bootloader的编写,讲解完M3的bootloader编写会扩展讲解M0的bootloader...... 内容绝对实在,对于学习以及工作都会有很大的帮助。最终的目的就是让大家学会快速开发STM32并收获与STM32有关的实用技术知识。

    795 人正在学习 去看看 李凯龙

图像处理系列笔记: https://blog.csdn.net/qq_33208851/article/details/95335809


傅里叶变换是一种函数在空间域和频率域的变换,从空间域到频率域的变换是傅里叶变换,而从频率域到空间域是傅里叶的反变换
时域与频域

  • 频域(frequency domain)
    是指在对函数或信号进行分析时,分析其和频率有关部份,而不是和时间有关的部份,和时域一词相对。
  • 时域
    是描述数学函数或物理信号对时间的关系。例如一个信号的时域波形可以表达信号随着时间的变化。若考虑离散时间,时域中的函数或信号,在各个离散时间点的数值均为已知。若考虑连续时间,则函数或信号在任意时间的数值均为已知。在研究时域的信号时,常会用示波器将信号转换为其时域的波形。
  • 两者相互间的变换
    时域(信号对时间的函数)和频域(信号对频率的函数)的变换在数学上是通过积分变换实现。对周期信号可以直接使用傅立叶变换,对非周期信号则要进行周期扩展,使用拉普拉斯变换。

信号在频率域的表现
在频域中,频率越大说明原始信号 变化速度越快;频率越小说明原始信号越平缓。当频率为0时,表示直流信号,没有变化。因此,频率的 大小反应了信号的变化快慢。高频分量解释信号的突变部分,而低频分量决定信号的整体形象。
在 图像处理中,频域反应了图像在空域灰度变化剧烈程度,也就是图像灰度的变化速度,也就是图像的梯度大小。对图像而言,图像的边缘部分是突变部分,变化较 快,因此反应在频域上是高频分量;图像的噪声大部分情况下是高频部分;图像平缓变化部分则为低频分量。也就是说,傅立叶变换提供另外一个角度来观察图像, 可以将图像从灰度分布转化到频率分布上来观察图像的特征。书面一点说就是,傅里叶变换提供了一条从空域到频率自由转换的途径。对图像处理而言,以下概念非 常的重要:

  • 图像高频分量:图像突变部分;在某些情况下指图像边缘信息,某些情况 下指噪声,更多是两者的混合;
  • 低频分量:图像变化平缓的部分,也就是图像轮廓信息
  • 高通滤波器:让图像使低频分量抑制,高频分量通过
  • 低通滤波器:与高通相反,让图像使高频分量抑制,低频分量通过
  • 带通滤波器:使图像在某一部分 的频率信息通过,其他过低或过高都抑制
  • 还有个带阻滤波器,是带通的反。

1. 傅里叶变换及其反变换

1.0 什么是傅里叶变换

  1. 什么是傅里叶变换?
    也称作傅立叶变换,表示能将满足一定条件的某个函数表示成三角函数(正弦和/或余弦函数)或者它们的积分的线性组合。在不同的研究领域,傅立叶变换具有多种不同的变体形式,如连续傅立叶变换和离散傅立叶变换。傅里叶变换是一种分析信号的方法,它可分析信号的成分,也可用这些成分合成信号。许多波形可作为信号的成分,比如正弦波、方波、锯齿波等,傅里叶变换用正弦波作为信号的成分。
    傅里叶变换的实质是将一个信号分离为无穷多多正弦/复指数信号的加成,也就是说,把信号变成正弦信号相加的形式——既然是无穷多个信号相加,那对于非周期信号来说,每个信号的加权应该都是零——但有密度上的差别,你可以对比概率论中的概率密度来思考一下——落到每一个点的概率都是无限小,但这些无限小是有差别的所以,傅里叶变换之后,横坐标即为分离出的正弦信号的频率,纵坐标对应的是加权密度
  2. 傅里叶变换有什么用呢?
    举例说明:傅里叶变换可以将一个时域信号转换成在不同频率下对应的振幅及相位,其频谱就是时域信号在频域下的表现,而反傅里叶变换可以将频谱再转换回时域的信号。最简单最直接的应用就是时频域转换,比如在移动通信的LTE系统中,要把接收的信号从时域变成频域,就需要使用FFT(快速傅里叶变换)。又例如对一个采集到的声音做傅立叶变化就能分出好几个频率的信号。比如南非世界杯时,南非人吹的呜呜主拉的声音太吵了,那么对现场的音频做傅立叶变化(当然是对声音的数据做),会得到一个展开式,然后找出呜呜主拉的特征频率,去掉展开式中的那个频率的sin函数,再还原数据,就得到了没有呜呜主拉的嗡嗡声的现场声音。而对图片的数据做傅立叶,然后增大高频信号的系数就可以提高图像的对比度。同样,相机自动对焦就是通过找图像的高频分量最大的时候,就是对好了。

1.1 为什么要在频率域研究图像增强?

  • 可以利用频率成分和图像外表之间的对应关系。一些在空间域表述困难的增强任务,在频率域中变得非常普通
  • 滤波在频率域更为直观,它可以解释空间域滤波的某些性质
  • 可以在频率域指定滤波器,做反变换,然后在空间域使用结果滤波器作为空间域滤波器的指导
  • 一旦通过频率域试验选择了空间滤波,通常实施都在空间域进行

1.2 傅里叶变换及反转

1.2.1 一维连续傅里叶变换及反变换

单变量连续函数f(x)的傅里叶变换F(u)定义为:
在这里插入图片描述
其中,j = 根号(-1)=±i
给定F(u),通过傅里叶反变换可以得到f(x):
在这里插入图片描述

1.2.2 二维连续傅里叶变换及反变换

二维连续函数f(x,y)的傅里叶变换F(u,v)定义为:
在这里插入图片描述
如果f(x,y)是实函数,它的傅里叶变换是对称的,即

F(u,v) = F(− u,−v)
傅里叶变换的频率谱是对称的
|F(u,v)| =| F(− u,−v)|

给定F(u,v),通过傅里叶反变换可以得到 f(x,y):
在这里插入图片描述

1.2.3 一维离散傅里叶变换及反变换

单变量离散函数f(x)(x=0,1,2,…,M-1)的傅里叶变换F(u)定义为:
在这里插入图片描述
其中,u=0,1,2,…,M-1
从欧拉公式:
在这里插入图片描述
在这里插入图片描述
给定F(u),通过傅里叶反变换可以得到f(x):
在这里插入图片描述
其中,x=0,1,2,…,M-1

1.2.4 二维离散傅里叶变换及反变换

图像尺寸为M×N的函数f(x,y)的DFT为:
在这里插入图片描述
u=0,1,2,…,M-1, v=0,1,2,…,N-1
给出F(u,v),可通过反DFT得到f(x,y):
在这里插入图片描述
x=0,1,2,…,M-1, y=0,1,2,…,N-1
注:u和v是频率变量,x和y是空间域图像变量
F(0,0)表示:
在这里插入图片描述
这说明:假设f(x,y)是一幅图像,在原点的傅里叶变换等于图像的平均灰度级(M*N是总的像素点,f(x,y)是(x,y)点的灰度值,将所有的像素点的灰度值求和然后除以总的个数即为平均灰度值)

1.2.5 傅里叶变换的一维极坐标表示

在这里插入图片描述
幅度或频率谱为:
在这里插入图片描述
R(u)和I(u)分别是F(u)的实部和虚部
相角或相位谱为:
在这里插入图片描述
功率谱为:
在这里插入图片描述
f(x)的离散表示:
在这里插入图片描述
F(u)的离散表示:
在这里插入图片描述

1.2.6 傅里叶变换的二维极坐标表示

二维DFT的极坐标表示:
在这里插入图片描述
幅度或频率谱为:
在这里插入图片描述
R(u,v)和I(u,v)分别是F(u,v)的实部和虚部
相角或相位谱为:
在这里插入图片描述
功率谱为:
在这里插入图片描述
F(u,v)的原点变换:
在这里插入图片描述
用(-1)x+y乘以f(x,y),将F(u,v)原点变换到率坐标下的(M/2,N/2),它是M×N区域的中心
u=0,1,2,…,M-1, v=0,1,2,…,N-1

2. 傅里叶变换的性质

2.1 平移性

以⇔表示函数和其傅里叶变换的对应性
在这里插入图片描述
注:u和v是频率变量,x和y是空间域图像变量
公式(1)表明将f(x,y)与一个指数项相乘就相当于把其变换后的频域中心f(u,v) 移动到新的位置 f(u-uo,v-v0)
公式(2)表明将F(u,v)与一个指数项相乘就相当于把其变换后的空域中心f(x,y) 移动到新的位置 f(x-x0,y-y0)
公式(2)表明对f(x,y)的平移不影响其傅里叶变换的幅值

当u0=M/2且v0=N/2,
在这里插入图片描述
带入(1)和(2),得到
在这里插入图片描述

2.2 分配律

傅里叶变换对加法满足分配律,但对乘法则不满足:
在这里插入图片描述

2.3 尺度变换(缩放)

给定2个标量a和b,可以证明对傅里叶变换下列2个公式成立:
在这里插入图片描述

2.4 旋转性

引入极坐标 x = r cosθ, y = rsinθ,u =ω cosϕ,v =ωsinϕ
将f(x,y)和F(u,v)转换为 f (r,θ ) 和F(ω,ϕ)。将它们带入傅里叶变换对得到:
在这里插入图片描述
f(x,y)旋转角度θ 0,F(u,v)也将转过相同的角度
F(u,v)旋转角度θ 0,f(x,y)也将转过相同的角度

2.5 周期性和共轭对称性

在这里插入图片描述
尽管F(u,v)对无穷多个u和v的值重复出现,但只需根据在任一个周期里的N个值就可以从F(u,v)得到f(x,y)
只需一个周期里的变换就可将F(u,v)在频域里完全确定
同样的结论对f(x,y)在空域也成立

如果f(x,y)是实函数,则它的傅里叶变换具有共轭对称性
在这里插入图片描述
其中,F*(u,v)为F(u,v)的复共轭(当两个复数实部相等,虚部互为相反数时,这两个复数叫做互为共轭复数)

2.6 分离性

在这里插入图片描述
F(x,v)是沿着f(x,y)的一行所进行的傅里叶变换。当x=0,1,…,M-1,沿着f(x,y)的所有行计算傅里叶变换。

二维傅里叶变换的全过程
在这里插入图片描述
 先通过沿输入图像的每一行计算一维变换
 再沿中间结果的每一列计算一维变换
 可以改变上述顺序,即先列后行
 上述相似的过程也可以计算二维傅里叶反变换

2.7 平均值

由二维傅里叶变换的定义:
在这里插入图片描述
所以,
在这里插入图片描述
上式说明:如果f(x,y)是一幅图像,在原点的傅里叶变换即等于图像的平均灰度级

2.8 卷积理论

卷积是空间域过滤和频率域过滤之间的纽带
大小为M×N的两个函数f(x,y)和h(x,y)的离散卷积:
在这里插入图片描述
卷积定理:
在这里插入图片描述

2.9 相关性理论

相关的重要应用在于匹配:确定是否有感兴趣的物体区域
 f(x,y)是原始图像
 h(x,y)作为感兴趣的物体或区域(模板)
 如果匹配,两个函数的相关值会在h找到f中相应点的位置上达到最大
大小为M×N的两个函数f(x,y)和h(x,y)的相关性定义为:
在这里插入图片描述
f* 表示f的复共轭。对于实函数,f*=f
相关定理:
在这里插入图片描述
自相关理论
在这里插入图片描述

3. 快速傅里叶变换(FFT)

采用快速傅里叶变换(FFT)算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少,特别是被变换的抽样点数N越多,FFT算法计算量的节省就越显著。
函数或信号可以透过一对数学的运算子在时域及频域之间转换。和傅里叶变换作用一样。

3.1 为什么需要快速傅里叶变换?

人们想让计算机能处理信号 但由于信号都是连续的、无限的,计算机不能处理,于是就有了傅里叶级数、傅里叶变换,将信号由时域变到频域,把一个信号变为有很多个不同频率不同幅度的正弦信号组成,这样计算机就能处理了,但又由于傅里叶变换中要用到卷积计算,计算量很大,计算机也算不过来,于是就有了快速傅里叶变换,大大降低了运算量,使得让计算机处理信号成为可能。快速傅里叶变换是傅里叶变换的快速算法而已,主要是能减少运算量和存储开销,对于硬件实现特别有利。
在这里插入图片描述

  • 对u的M个值中的每一个都需进行M次复数乘法(将f(x)与 e− j2πux / M 相乘)和M-1次加法,即复数乘法和加法的次数都正比于M2
  • 快速傅里叶变换(FFT)则只需要Mlog2M次运算
  • FFT算法与原始变换算法的计算量之比是log2M/M,如M=1024≈103,则原始变换算法需要106次计算,而FFT需 要104次计算,FFT与原始变换算法之比是1:100
    只考虑一维的情况,根据傅里叶变换的分离性可知,二维傅里叶变换可由连续2次一维傅里叶变换得到

3.2 FFT算法基本思想

FFT算法基于一个叫做逐次加倍的方法。通过推导将原始傅里叶转换成两个递推公式:
在这里插入图片描述

3.3 FFT公式推导

在这里插入图片描述
假设M的形式是
M = 2n, n为正整数。因此,M可以表示为:M = 2K 。将M=2K带入上式:
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
特性:

  • 一个M个点的变换,能够通过将原始表达式分成两个部分来计算
  • 通过计算两个(M/2)个点的变换。得Feven(u)和 Fodd(u)
  • 奇部与偶部之和得到F(u)的前(M/2)个值
  • 奇部与偶部之差得到F(u)的后(M/2)个值。且不需要额外的变换计算

3.4 归纳快速傅立叶变换的思想

(1)通过计算两个单点的DFT,来计算两个点的DFT, (2)通过计算两个双点的DFT,来计算四个点的DFT,…,以此类推
(3)对于任何N=2m的DFT的计算,通过计算两个N/2点的DFT,来计算N个点的DFT

3.5 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)}
解法:
首先分成奇偶两组,有:
{ 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) }
对输入数据的排序可根据一个简单的位对换规则进行:
如用x表示f(x)的1个自变量值,那么它排序后对应的值可通过把x表示成二进制数并对换各位得到。例如N=23,f(6)排序后为f(3),因为6=1102而0112 =3
把输入数据进行了重新排序,则输出结果是正确的次序。反之不把输入数据进行排序,则输出结果需要重新排序才能得到正确的次序
地址的排序:——按位倒序规则
例如:N = 23 = 8
在这里插入图片描述
2)计算顺序及地址增量:2n, n = 0,1,2…
在这里插入图片描述

4. 傅里叶变换的物理意义

  1. 图像的频率是表征图像中灰度变化剧烈程度的指标,是灰度在平面空间上的梯度(灰度变化得快频率就高,灰度变化得慢频率就低)。如:大面积的沙漠在图像中是一片灰度变化缓慢的区域,对应的频率值很低;而对于地表属性变换剧烈的边缘区域在图像中是一片灰度变化剧烈的区域,对应的频率值较高。傅立叶变换在实际中有非常明显的物理意义,设f是一个能量有限的模拟信号,则其傅立叶变换就表示f的谱。从纯粹的数学意义上看,傅立叶变换是将一个函数转换为一系列周期函数来处理的。从物理效果看,傅立叶变换是将图像从空间域转换到频率域,其逆变换是将图像从频率域转换到空间域。换句话说,傅立叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数,傅立叶逆变换是将图像的频率分布函数变换为灰度分布函数
  2. 傅立叶变换以前,图像(未压缩的位图)是由对在连续空间(现实空间)上的采样得到一系列点的集合,我们习惯用一个二维矩阵表示空间上各点,则图像可由z=f(x,y)来表示。由于空间是三维的,图像是二维的,因此空间中物体在另一个维度上的关系就由梯度来表示,这样我们可以通过观察图像得知物体在三维空间中的对应关系。为什么要提梯度?因为实际上对图像进行二维傅立叶变换得到频谱图,就是图像梯度的分布图,当然频谱图上的各点与图像上各点并不存在一一对应的关系,即使在不移频的情况下也是没有。傅立叶频谱图上我们看到的明暗不一的亮点,实际上图像上某一点与邻域点灰度值差异的强弱,即梯度的大小,也即该点的频率的大小(差异/梯度越大,频率越高,能量越低,在频谱图上就越 暗。差异/梯度越小,频率越低,能量越高,在频谱图上就越 亮。换句话说,频率谱上越亮能量越高,频率越低,图像差异越小/平缓)。一般来讲,梯度大则该点的亮度强,否则该点亮度弱频谱图,也叫功率图
    在这里插入图片描述
    在经过频谱中心化(用(-1)x+y乘以输入的图像函数在这里插入图片描述)后的频谱中,中间最亮的点是最低频率,属于直流分量(DC分量)(当频率为0时,表示直流信号,没有变化。在原点(u,v两个频率域变量均为零)的傅里叶变换即等于图像的平均灰度级,F(0,0)称做频 率谱的直流成分)。越往边外走,频率越高。所以,频谱图中的四个角和X,Y轴的尽头都是高频,如下图:
    在这里插入图片描述

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

参考文章:https://blog.csdn.net/EbowTang/article/details/39004979
图像处理系列笔记: https://blog.csdn.net/qq_33208851/article/details/95335809

没有更多推荐了,返回首页