图像处理之图像傅里叶变换

一、傅里叶变换的物理意义

从纯粹的数学意义上看，傅里叶变换是将-一个图像函数转换为一系列周期函数来处理的；从物理效果看，傅里叶变换是将图像从空间域转换到频率域，其逆变换是将图像从频率域转换到空间域。即傅里叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，傅里叶逆变换是将图像的频率分布函数变换为灰度分布函数。实际上对图像进行二维傅里叶变换得到频谱图就是图像梯度的分布图，傅里叶频谱图上看到的明暗不一的亮点，实际上图像上某一点 与邻域点差异的强弱，即梯度的大小，即该点的频率大小。如果频谱图中暗的点数更多，则实际图像是比较柔和的；反之，如果频谱图中亮的点数多，则实际图像是比较尖锐的，边界分明且边界两边像素差异较大。

二、傅里叶变换的定义及性质

1、傅里叶变换的定义

傅里叶变换是在以时间为自变量的“信号”与频率为自变量的“频谱”函数之间的某域研究中较复杂的问题在频域中变得简单起来，从而简化其分析过程；另一方面使信号与系统的物理本质在频域中能更好地被揭示出来。当自变量“时间”或“频率”为连续形式和离散形式的不同组合时，就可以形成各种不同的傅里叶变换对，即“信号”与“频谱”的对应关系。傅里叶变换包含连续傅里叶变换、离散傅里叶变换、快速傅里叶变换和短时傅里叶变换等，在数字图像处理中使用的是二维离散傅里叶变换。

2、二维离散傅里叶变换性质

三、傅里叶变换的matlab实现

在MATLAB软件中，通过函数fft()进行一维离散傅里叶变换，通过函数ifft()进行一维离散傅立叶反变换。

在MATLAB中，通过函数fft2()进行二维离散傅里叶变换，通过函数ifft2()进行二维离散傅里叶反变换。函数fft()和fft2()的关系为fft2(X)=fft(fft(X).‘).’。

在MATLAB软件中，通过函数ifft2()进行二维快速傅立叶反变换，该函数和函数fft2()互为反函数。

通过fft2得到的频谱坐标原点位于左上角，函数fftshift()将变换后的坐标原点移到频谱窗口中央，坐标原点是低频，向外是高频。函数fftshift()能够进行傅里叶平移，函数ifftshift()能够进行傅里叶反平移。简单理解，ifftshift就是撤消fftshift的结果。

1、矩阵和图像的二维离散傅里叶变换

% 矩阵的二维离散傅里叶变换(变换后得到的矩阵和原矩阵大小相等)
I1=ones(4)
I2=[2 2 2 2;1 1 1 1;3 3 0 0;0 0 0 0]
J1=fft2(I1)
J2=fft2(I2)

% 图像的二维离散傅里叶变换
I=imread('cameraman.tif');
J=fft2(I);
K=abs(J/256);% 计算频率(频谱图即频率分布图)

subplot(121),imshow(I);
title('原始图像');
subplot(122),imshow(uint8(K));% 显示频谱图(通过fft2得到的频谱图中,坐标原点在窗口的左上角,窗口的四角分布低频部分)
title('灰度图像的频谱图');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

命令行结果：

I1 =

     1     1     1     1
     1     1     1     1
     1     1     1     1
     1     1     1     1
I2 =

     2     2     2     2
     1     1     1     1
     3     3     0     0
     0     0     0     0
J1 =

    16     0     0     0
     0     0     0     0
     0     0     0     0
     0     0     0     0
J2 =

  18.0000 + 0.0000i   3.0000 - 3.0000i   0.0000 + 0.0000i   3.0000 + 3.0000i
   2.0000 - 4.0000i  -3.0000 + 3.0000i   0.0000 + 0.0000i  -3.0000 - 3.0000i
  10.0000 + 0.0000i   3.0000 - 3.0000i   0.0000 + 0.0000i   3.0000 + 3.0000i
   2.0000 + 4.0000i  -3.0000 + 3.0000i   0.0000 + 0.0000i  -3.0000 - 3.0000i
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

实现效果：

2、通过fftshift()对矩阵和图像进行平移

% 通过函数fftshift()进行平移
close all;
clear all;
clc;
N=0:4
X=fftshift(N)            % 平移
Y=fftshift(fftshift(N))  % 平移后再进行平移
Z=ifftshift(fftshift(N)) % 平移后再进行反平移(恢复原来结果)

% 图像进行傅里叶变换和平移

% 在以后进行的傅里叶变化都进行平移！！！
I=imread('peppers.png');
J=rgb2gray(I);
K=fft2(J);     % 傅里叶变换
K=fftshift(K); % 平移
L=abs(K/256);

subplot(121),imshow(J);
title('原始图像');
subplot(122),imshow(uint8(L));% 显示频谱图(图像的能量主要集中在低频部分即频谱图的中央,四个角的高频部分幅值非常小)
title('灰度图像的频谱图');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

命令行结果：

N =

     0     1     2     3     4
X =

     3     4     0     1     2
Y =

     1     2     3     4     0
Z =

     0     1     2     3     4
1
2
3
4
5
6
7
8
9
10
11
12

实现效果：

3、图像变亮后进行傅里叶变换

close all;
clear all;
clc;
I=imread('peppers.png');
J=rgb2gray(I);
% 将灰度图像的数据矩阵乘以e使灰度图像变亮
% 注意本例中乘以e后像素值最大为255,若某些例子中乘以e后可能会产生大于255的像素值越界,需要找到并重新赋值即J(find(J>255))=255;
J=J*exp(1);        
K=fft2(J);          % 傅里叶变换
K=fftshift(K);      % 平移
L=abs(K/256);

subplot(121),imshow(J);
title('变亮后的图像');
subplot(122),imshow(uint8(L));% 显示频谱图
title('变亮后图像对应的频谱图');
% 对比fourier2和fourier3的频谱图,图像变亮后中央低频部分变大
% 频谱图中央低频部分代表了灰度图像的平均亮度
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

实现效果：

4、图像旋转后进行傅里叶变换

% 图像旋转后进行傅里叶变换
close all;
clear all;
clc;
I=imread('peppers.png');
J=rgb2gray(I);
% 图像旋转B=imrotate(A,angle,method)中method包括nearest最临近插值-默认值、bilinear双线性的和bicubic双三次的
J=imrotate(J,45,'bilinear');
K=fft2(J);      % 傅里叶变换
K=fftshift(K);  % 平移
L=abs(K/256);

subplot(121),imshow(J);
title('旋转后的图像');
subplot(122),imshow(uint8(L));% 显示频谱图(注意图像发生旋转,频谱图相应也会发生旋转)
title('旋转后图像对应的频谱图');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

实现效果：

5、图像中添加高斯噪声后进行傅里叶变换

% 图像中添加高斯噪声后进行傅里叶变换
close all;
clear all;
clc;
I=imread('peppers.png');
J=rgb2gray(I);
J=imnoise(J,'gaussian',0,0.01);% 添加高斯噪声
K=fft2(J);     % 傅里叶变换
K=fftshift(K); % 平移
L=abs(K/256);

subplot(121),imshow(J);
title('添加高斯噪声后的图像');
subplot(122),imshow(uint8(L));% 显示频谱图
title('添加高斯噪声后图像对应的频谱图');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

实现效果：

6、灰度图像的傅里叶变换和反变换

% 灰度图像的傅里叶变换和反变换
close all;
clear all;
clc;
I=imread('onion.png');
J=rgb2gray(I);
K=fft2(J);    % 傅里叶变换
L=fftshift(K);% 平移
M=ifft2(K);   % 傅里叶反变换

subplot(131),imshow(J);
title('原始图像');
subplot(132),imshow(uint8(abs(L/256)));% 经过傅里叶变换后得到的频谱图
title('傅里叶频谱');
subplot(133),imshow(uint8(M));         % 对傅里叶变换系数进行傅里叶反变换后得到的灰度图像
title('傅里叶反变换后的图像');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

实现效果：

7、灰度图像的幅值谱和相位谱

% 灰度图像的幅值谱和相位谱
close all;
clear all;
clc;
I=imread('peppers.png');
J=rgb2gray(I);
K=fft2(J);    % 傅里叶变换
L=fftshift(K);% 平移
fftr=real(L); % real函数求复数的实部
ffti=imag(L); % imag函数求复数的虚部
% 令R(u,v)和I(u,v)表示傅里叶变换F(u,v)的实部和虚部,则原图像函数f(x,y)傅里叶变换的幅值谱定义为|F(u,v)|=[(R(u,v))^2+(I(u,v))^2]^(1/2)
A=sqrt(fftr.^2+ffti.^2);% 幅值谱
A=(A-min(min(A)))/(max(max(A))-min(min(A)))*255;% 归一化到0-255(y=(x-min)/(max-min)为归一到0-1之间）)
B=angle(K);% 求复数相角
subplot(121),imshow(A);      % 灰度图像傅里叶变换后的幅值谱
title('傅里叶变换后的幅值谱');
subplot(122),imshow(real(B));% 灰度图像经过傅里叶变换后得到的相位谱
title('傅里叶变换后的相位谱');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

实现效果：

8、编程实现二维离散傅里叶变换

% 编程实现二维离散傅里叶变换
close all;
clear all;
clc;
I=imread('onion.png');
J=rgb2gray(I);
J=double(J);
s=size(J);
M=s(1);N=s(2);% 获取图像的行数和列数
for u=0:M-1
    for v=0:N-1
        k=0;
        for x=0:M-1
            for y=0:N-1
                k=J(x+1,y+1)*exp(-j*2*pi*(u*x/M+v*y/N))+k;% 根据二维离散傅里叶变换公式所得(注意公式中j表示复数的虚部比如a+bj)
            end
        end
        F(u+1,v+1)=k;% 傅里叶变换结果
    end
end
K=fft2(J);
% K和F图的结果基本相同,但是编程实现的二维离散傅里叶变换计算速度比较慢,fft2函数计算速度比较快
subplot(121),imshow(K);
title('函数fft2得到的傅里叶频谱');
subplot(122),imshow(F);
title('编程实现的傅里叶频谱');
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

实现效果：

9、图像处理之理想低通滤波器、巴特沃斯低通滤波器和高斯低通滤波器的matlab实现去噪

10、图像处理之理想高通滤波器、巴特沃斯高通滤波器和高斯高通滤波器的matlab简单实现

11、图像处理之理想带阻滤波器、巴特沃斯带阻滤波器和高斯带阻滤波器的matlab实现去噪

由于刚刚开始学习图像处理，对于很多知识理解不到位。如有错误，恳请指正，任重而道远，慢慢加油！