【数字图像处理笔记】Matlab实现离散傅立叶变换 （二）

• 💌 所属专栏：【数字图像处理笔记】

• 😀 作  者：我是夜阑的狗🐶

• 🚀 个人简介：一个正在努力学技术的CV工程师，专注基础和实战分享 ，欢迎咨询！

• 💖 欢迎大家：这里是CSDN，我总结知识的地方，喜欢的话请三连，有问题请私信 😘 😘 😘
  
  

您的点赞、关注、收藏、评论，是对我最大的激励和支持！！！🤩 🤩 🤩

---

前言

  大家好，又见面了，我是夜阑的狗🐶，本文是专栏【数字图像处理笔记】专栏的第2篇文章;
  这是今天学习到数字图像处理笔记 – Matlab实现离散傅立叶变换 💖💖💖，开启新的征程，记录最美好的时刻🎉，每天进步一点点。
  专栏地址：【数字图像处理笔记】 ， 此专栏是我是夜阑的狗对数字图像处理学习过程的总结，希望能够加深自己的印象，以及帮助到其他的小伙伴😉😉。
  如果文章有什么需要改进的地方还请大佬不吝赐教👏👏。

一、原理

  要在数字图像处理中应用傅立叶变换， 还需要解决两个问题：

• 问题一：在数学中进行傅立叶变换的f(x)为连续（模拟）信号， 而计算机处理的是数字信号（图像数据）；
• 问题二：数学上采用无穷大概念，而计算机只能进行有限次计算。

  通常， 将受这种限制的傅立叶变换称为离散傅立叶变换（Discrete Fourier Transform，DFT)。离散傅立叶变换的定义如下。

• 离散傅立叶正变换：

F ( u ) = ∑ x = 0 N − 1 f ( x ) e − j 2 π u x N F(u) = \sum_{

}^{N-1}f(x)e^{-j{\frac {2 \pi ux} N}} F(u)=x=0​∑N−1​f(x)e−jN2πux​
  其中 $u=0,1,2,\Lambda ,N-1$；

• 离散傅立叶逆变换：

f ( u ) = 1 N ∑ u = 0 N − 1 F ( x ) e − j 2 π u x N f(u) = {\frac 1 N}\sum_{

}^{N-1}F(x)e^{-j{\frac {2 \pi ux} N}} f(u)=N1​u=0​∑N−1​F(x)e−jN2πux​
  其中 $x=0,1,2,\Lambda ,N-1$；

二、实现过程

1、离散傅立叶变换（直接调用Matlab内置函数）：

  代码：

I=imread('cameraman.tif');
% I=rgb2gray(I);	%将图像转换为灰度图
I=im2double(I);     %将图像的数据格式转换为double型的 
F=fft2(I);               %傅立叶变换
%不进行象限转换
F1=abs(F);             %傅里叶变换的模
T1=log(F1+1);		 %数据范围压缩
%进行象限转换
F=fftshift(F);          %对傅里叶变换后的图像进行象限转换
F=abs(F);               %傅里叶变换的模
T=log(F+1);           %数据范围压缩
figure;                    %显示
subplot(2,2,1);imshow(I);
subplot(2,2,2);imshow(I);
subplot(2,2,3);imshow(T,[]); 
subplot(2,2,4);imshow(T1,[]);

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

  运行结果如下图所示：

  从图中可以看出，原图像经过快速傅立叶变换 fft 后可以得到频谱图；以及通过 ffshift 函数将零频点移到频谱的中间，可以更直观观察傅立叶变换，傅立叶变换能使从空间域与频率域两个不同角度来分析图像问题。还可以通过 ifftshift 函数从频率域反映射到空间域，傅里叶逆变换基本还原图像。
  在进行傅立叶变换的具体实现中，要注意对fft变换后的频域矩阵数值比较大，一般采用log(abs[F]+1) 来对缩放频域，保证映射都在正数范围内，abs 函数可以求出复数矩阵的模缩放处理后进行显示。
  频谱图表示信号变化的快慢剧烈程度，通过分析不同图像和频谱图的关系时发现，高频信号对应图像的边缘细节，低频信号对应图像的大致概貌和轮廓。因此，对于细节边缘清晰或者噪声较多的图像，经过傅立叶变换后，靠外边的高频信号更多；而对于平滑、具备大体轮廓的图像，经过傅立叶变换后，靠中间的低频信号更多。
  所以，频谱图中的四个角和X,Y轴的尽头都是高频。可以以发现频谱图中的较亮白线应该是对应相机三脚架的边缘以及衣服的边缘。

2、离散傅立叶变换（自己编写）：

  代码：

clc; clear;

data = imread('cameraman.tif');  % 数据——最好比卷积核的尺寸大
data = im2double(data); 
% data = rgb2gray(data);     % rgb转为灰度图像
subplot(1,3,1);
imshow(data);
title('原始图像')

zidai = fft2(data);   % matlab自带函数，来用对比
T1= abs(zidai);    %傅里叶变换的模
T1= fftshift(T1);   %对傅里叶变换后的图像进行象限转换
T1 = log(T1+1);	  %数据范围压缩
subplot(1,3,2);
% imshow(real(zidai));  % 一般只要实部
imshow(T1,[]); 
title('fft2');

size_data = size(data);
M = size_data(1);  % 图(原始数据矩阵)的长
N = size_data(2);  % 图(原始数据矩阵)的宽

% 下面是傅里叶正变换必备的一些矩阵:
Wm = exp(-j*2*pi/M);
Wn = exp(-j*2*pi/N); % 不同G中用不同的W
Em = zeros(M);
En = zeros(N);     % E是辅助计算矩阵
Gm = zeros(M)+Wm;
Gn = zeros(N)+Wn;  % G是计算时要用的矩阵
F = zeros(M,N);    % F是转换到频域的结果

% 对Gm的计算: 循环长度为M
fprintf('二维离散傅里叶变换开始:\n');
for row = 0:M-1
    for col = 0:M-1
        Em(row+1,col+1) = row * col;
        Gm(row+1,col+1) = Gm(row+1,col+1)^Em(row+1,col+1);
    end
end
% 对Gn的计算: 循环长度为N
for row = 0:N-1
    for col = 0:N-1
        En(row+1,col+1) = row * col;
        Gn(row+1,col+1) = Gn(row+1,col+1)^En(row+1,col+1);
    end
end

% F = real(Gm*data*Gn);  % F = Gm*f*Gn是计算公式，一般只要实部
F = Gm*data*Gn;
subplot(1,3,3);
T=fftshift(F);   %对傅里叶变换后的图像进行象限转换
T=abs(T);        %傅里叶变换的模
T=log(T+0.7);  %数据范围压缩
% imshow(F);
imshow(T,[]);
title('myfft2');

error = sum(sum((real(F)-real(zidai)).^2));
if error < 10^(-10)
    fprintf('自带与手写结果一致!\n');
else
    fprintf('不一致!\n');
end

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

  运行结果如下图所示：

  从图中的对比结果可以看出，通过编程实现傅里叶变换与 Matlab 自带的傅里叶变换函数得到的频谱图基本一致，说明了该程序的正确性。

3、离散傅里叶逆变换（自己编写）：

  代码：

% 鉴于正向fft2手写与自带结果一致；
% ifft2的输入就直接用自带的fft2的结果。
clc; clear;

data = imread('cameraman.tif');  % 数据——最好比卷积核的尺寸大
data = im2double(data); 
% data = rgb2gray(data);     % rgb转为灰度图像
subplot(1,3,1);
imshow(data);
title('原始图像')

F = fft2(data);
T1= abs(F);
T1= fftshift(T1);
T1 = log(T1+1);	
subplot(1,3,2);
imshow(T1,[]);
% imshow(real(F));  % 一般画图只要实部, 作为输入时实虚都要！！
title('fft2');

% s = ifft2(F);
% subplot(1,3,3);
% imshow(s);
% return;

size_data = size(F);
M = size_data(1);  % 图(原始数据矩阵)的长
N = size_data(2);  % 图(原始数据矩阵)的宽

% 下面是傅里叶逆变换必备的一些矩阵:
Wm = exp(-j*2*pi/M);
Wn = exp(-j*2*pi/N);  % 不同G中用不同的W
Em = zeros(M);
En = zeros(N);        % E是辅助计算矩阵
Gm = zeros(M)+Wm;
Gn = zeros(N)+Wn;  % G是计算时要用的矩阵
f = zeros(M,N);    % F是转换到频域的结果

% 对Gm的计算: 循环长度为M
fprintf('二维离散反傅里叶变换开始:\n');
for row = 0:M-1
    for col = 0:M-1
        Em(row+1,col+1) = -row * col;
        Gm(row+1,col+1) = Gm(row+1,col+1)^Em(row+1,col+1);
    end
end
Gm = Gm/M;
% 对Gn的计算: 循环长度为N
for row = 0:N-1
    for col = 0:N-1
        En(row+1,col+1) = -row * col;
        Gn(row+1,col+1) = Gn(row+1,col+1)^En(row+1,col+1);
    end
end
Gn = Gn/N;   % 注意：这个/N和上面的/M都是算完G之后才除以的！因为上面计算的时候是幂项变化！

f = real(Gm*F*Gn);  % f = Gm*F*Gn是计算公式，一般只要实部
subplot(1,3,3);
imshow(f);
title('逆变换myidft2');

error = sum(sum((real(f)-real(data)).^2));
if error < 10^(-10)
    fprintf('反变换后与原图一致!\n');
else
    fprintf('不一致!\n');
end

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

  运行结果如下图所示：

  从图中的傅里叶逆变换处理结果可以看出，通过编程实现的傅里叶逆变换得到图像与原图基本一致，说明了该程序的正确性。
  通过本次操作基本掌握傅立叶变换的基本原理和实现算法。在实际过程中遇到以下几个问题及解决方法：

• （1）要傅里叶变换频谱图进行分析，频谱没有中心化。解决方法：调用 fftshift() 函数对频谱图进行中心化。

---

总结

  感谢观看，如果觉得有帮助，请给文章点个赞吧，让更多的人看到。🌹 🌹 🌹

  也欢迎你，关注我。👍 👍 👍

  原创不易，还希望各位大佬支持一下，你们的点赞、收藏和留言对我真的很重要！！！💕 💕 💕 最后，本文仍有许多不足之处，欢迎各位认真读完文章的小伙伴们随时私信交流、批评指正！下期再见。🎉

更多专栏订阅:

• 😀 【LeetCode题解（持续更新中）】
• 🌼 【鸿蒙系统】
• 👑 【Python脚本笔记】
• 🚝 【Java Web项目构建过程】
• 💛 【微信小程序开发教程】
• ⚽ 【JavaScript随手笔记】
• 🤩 【大数据学习笔记（华为云）】
• 🦄 【程序错误解决方法（建议收藏）】
• 🚀 【软件安装教程】

订阅更多，你们将会看到更多的优质内容！！