傅里叶变换的相关实验——matlab实现

图像变换是数字图像处理中常用的技术，在图像增强、图像复原、图像压缩编码等数字图像处理中，都会用到图像变换技术，傅里叶变换是数字图像处理中应用最广的一种变换。

1. 从空间域到频率的唯一途径是进行傅里叶变换；

2. 傅里叶变换有明确的物理意义，即任何时域连续的复杂波形，通过傅里叶变换都可以变成一系列简谐波之和；

3. 傅里叶变换的频谱分析是全局性的；

4. 傅里叶变换在频率域里是具有实部和虚部的复函数；

5. 一维离散傅里叶频谱：$|F(u)|=\sqrt{R^2(u)+I^2(u)}$ ；能量谱：$E(u)=|F(u){|}^2$ ；

   相位谱：$\varphi (u)=\arctan (I(u)/R(u))$；

一、实验目的

（1）通过实验进一步加深对图像傅立叶变换的理解；
（2）计算离散图像的傅里叶变换；
（3）掌握图像的傅里叶频谱图及离散傅里叶变换性质；
（4）掌握 MATLAB中的傅立叶变换函数；
（5）实现数字图像的傅立叶变换与反变换；
（6）了解相位谱和幅值谱的意义，并且分别重构；
（7）实现双谱重构实验

二、实验原理

对数字图像$f(x,y)$，其二维离散傅里叶变换定义为：
$$F(u,v)=\frac{1}{\sqrt{MN}}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y)e^{-j2\pi (\frac{ux}{M}+\frac{vy}{N})}$$
其中，u=0,1,…,M-1，v=0,1,…,N-1。

二维离散傅里叶反变换定义为：
$$f(x,y)=\frac{1}{\sqrt{MN}}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}F(u,v)e^{j2\pi (\frac{ux}{M}+\frac{vy}{N})}$$
其中，x=0,1,…,M-1，y=0,1,…,N-1。

傅里叶频谱：
$$|F(u,v)|=\sqrt{R^2(u,v)+I^2(u,v)}$$
能量谱：
$$E(u,v)=|F(u,v){|}^2$$
相位谱：
$$\varphi (u,v)=\arctan \frac{I(u,v)}{R(u,v)}$$
注意：正反傅里叶变换的唯一区别是幂的符号。

图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。如：大面积的沙漠在图像中是一片灰度变化缓慢的区域，对应的频率值很低；而对于地表属性变换剧烈的边缘区域在图像中是一片灰度变化剧烈的区域，对应的频率值较高。傅里叶变换在实际中有非常明显的物理意义，设 f 是一 个能量有限的模拟信号，则其傅里叶变换就表示 f 的谱。
从纯粹的数学意义上看，傅里叶变换是将一个函数转换为一系列周期函数来处理的。从物理效果看，傅里叶变换是将图像从空间域转换到频率域，其逆变换是将图像从频率域转换到空间域。换句话说，傅里叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，傅里叶逆变换是将图像的频率分布函数变换为灰度分布函数。

不同频率信息在图像结构中有不同的作用：

• 图像的主要成分是低频信息，它形成了图像的基本灰度等级，对图像结构的决定作用较小；
• 中频信息决定了图像的基本结构，形成了图像的主要边缘结构；
• 高频信息形成了图像的边缘和细节，是在中频信息上对图像内容的进一步强化。

简短概括：

• 图像的高频，意味着灰度变化剧烈
• 图像的低频，意味着灰度变化平坦

另外说明以下几点：
图像经过二维傅里叶变换后，其变换系数矩阵表明：

1. 若变换矩阵 Fn原点设在中心，其频谱能量集中分布在变换系数矩阵的中心附近。若所用的二维傅里叶变换矩阵 Fn 的原点设在左上角，那么图像信号能量将集中在系数矩阵的四个角上。这是由二维傅里叶 变换本身性质决定的。同时也表明一般图像能量集中低频区域。
2. 变换之后的图像在原点平移之前四角是低频，最亮，平移之后中间部分是低频，最亮，亮度大说明
   低频的能量大（幅角比较大）。

三、实验准备

了解matlab实现傅里叶变换的相关函数

函数名	作用
fft2	二维快速傅里叶正变换
ifft2	二维快速傅里叶反变换
fftshift	将傅里叶变换原点移到中心
ifftshift	fftshift的反变换
real	复数的实部
imag	复数的虚部

注意：傅里叶变换一般都将利用平移性质，将中心移到零点。

1. fft2
   
2. ifft2
   
3. fftshif
   
4. ifftshift
   

四、实验内容

1.图像傅里叶正反变换实验

利用 MATLAB 对图像进行离散傅里叶正变换及反变换，分析图像的傅里叶频谱

（1）首先， 生成一幅大小为 512*512 的黑色背景的中间叠加一个尺寸为 40*40 的白色矩阵的图像实现傅里叶正反变换；

matlab代码如下：

%% 生成一幅大小为 512*512 的黑色背景的中间叠加一个尺寸为 40*40 的白色矩阵的图像实现傅里叶正反变换；
clc;clear;close all;    % 初始化
i=zeros(512,512);   % 生成512*512为0的矩阵
i(236:276,236:276)=1;   % 中间40*40赋值为1
i=im2double(i); % 将数据类型转换为double，数值范围变为[0,1]
i1=fft2(i); % 快速傅里叶变换
i2=fftshift(i1); % 将傅里叶变换进行象限转换，即将中心移到零点
i3=abs(i2); % 求傅里叶变换的模，将复数的double变成实数的double
i4=log(i3+1);   % 将傅里叶变换结果的幅值映射到较小的正数
subplot(121);imshow(i);title('原图');
subplot(122);imshow(i4,[]);title('傅里叶变换频谱'); %显示傅里叶变换结果图，映射到[0,1]

运行结果如下：

（2）将相应的数字图像进行傅立叶变换并且反变换

matlab代码如下：

%% 将相应的数字图像进行傅立叶变换并且反变换
clc;clear;close all;    % 初始化
i=imread('C:\Users\xcz\Desktop\matlab_example\Images\lena.bmp');
i=rgb2gray(i);  % 将真彩色图像转换为灰度图像
i=im2double(i); % 将数据类型转换为double，数值范围变为[0,1]
i1=fft2(i); % 进行快速傅里叶变换
i2=fftshift(i1);  % 将傅里叶变换进行象限转换，即将中心移到零点
i3=abs(i2); % 求傅里叶变换结果的模，即傅里叶频谱
i4=log(i3+1);   % 将傅里叶变换结果映射到较小的正数
i5=ifftshift(i2);   % 将变换的象限转换回去，即将中心复位
i6=ifft2(i5);   % 快速傅里叶反变换
subplot(131);imshow(i);title('原图');
subplot(132);imshow(i4,[]);title('傅里叶变换频谱');%显示傅里叶变换结果图，映射到[0,1]
subplot(133);imshow(i6);title('傅里叶逆变换图像');

运行结果如下：

2. 傅里叶变换性质实验

通过实验加深对傅立叶变换几个性质的理解，包括平移性质、尺度变换性质、旋转特性。

（1）平移性质：当空域中f(x,y)产生移动时，在频域中只发生相移，并不影响它的傅立叶变换的幅值。

%% 傅里叶变换平移性质：
%% 当空域中f(x,y)产生移动时，在频域中只发生相移，并不影响它的傅立叶变换的幅值。
clc;clear;close all;    % 初始化
i1=zeros(512,512);   % 生成512*512为0的矩阵
i1(236:276,236:276)=1;   % 中间40*40赋值为1
i1=im2double(i1); % 将数据类型转换为double，数值范围变为[0,1]
i2=zeros(512,512);   % 生成512*512为0的矩阵
i2(36:76,36:76)=1;   % 左上角40*40赋值为1 xy皆发生了移动
i2=im2double(i2); % 将数据类型转换为double，数值范围变为[0,1]
% i1为图1，i2为图2，相对图1 中间白色部分向左上角移动了
% 对图1进行傅里叶变换
i1_fft=fftshift(fft2(i1));  % 进行傅里叶变换，并且象限转换
i1_margian=abs(i1_fft); % 求模，即幅值谱
i1_phase=angle(i1_fft); % 相位谱
% 对图2进行傅里叶变换
i2_fft=fftshift(fft2(i2));  % 进行傅里叶变换，并且象限转换
i2_margian=abs(i2_fft); % 求模，即幅值谱
i2_phase=angle(i2_fft); % 相位谱
subplot(231);imshow(i1);title('图1');
subplot(232);imshow(i2);title('图2');
subplot(233);imshow(log(i1_margian+1),[]);title('图1幅值谱');  % 参数'[]'是为了将其值线性拉伸
subplot(234);imshow(i1_phase,[]);title('图1相位谱');   
subplot(235);imshow(log(i2_margian+1),[]);title('图2幅值谱');
subplot(236);imshow(i2_phase,[]);title('图2相位谱');

运行结果如下：

由结果可知，图像中某一部位的位置改变并不会改变图像的幅度谱，因为幅度谱表示各像素点的亮度信息，但具体点是被打乱的，所以幅度谱并不会改变。而相位谱因为记录了图像的所有点的相位信息，所以当图像中部分的位置发生改变时，相位谱会随之变化。

由实验可知，因为傅里叶变换的平移性质，当空域中f(x,y)产生移动时，在频域中只发生相移，并不影响它的傅立叶变换的幅值，即相位谱会发生改变，而幅值谱不会改变。

（2）尺度变换：当空域中f(x,y)产生缩放时，频域中傅立叶变换的幅值也相应的缩放。

%% 尺度变换：
%% 当空域中f(x,y)产生缩放时，频域中傅立叶变换的幅值也相应的缩放。
clc;clear;close all;    % 初始化
I=imread('C:\Users\xcz\Desktop\matlab_example\Images\lena.bmp');
I=im2double(I);     % 将数据类型转换为double，数值范围改为[0,1]
I=rgb2gray(I);  %将真彩色图像转换为灰度图像
I1=imresize(I,0.5,'bicubic');   % 将图像缩小0.1倍，使用三次内插法
I2=imresize(I,5,'bicubic');   % 将图像放大5倍，使用三次内插法

% 分别求三幅图的幅值谱
I_fft=fftshift(fft2(I));    %进行傅里叶变换，并且进行象限转换
I_margian=abs(I_fft);   % 幅值谱

I1_fft=fftshift(fft2(I1));   
I1_margian=abs(I1_fft);  

I2_fft=fftshift(fft2(I2));   
I2_margian=abs(I2_fft);   
subplot(231);imshow(I);title('原图');axis on;
subplot(232);imshow(I1);title('缩放0.5倍');axis on;
subplot(233);imshow(I2);title('放大5倍');axis on;
 % 参数'[]'是为了将其值线性拉伸,log是为了更好显示图，加1是为了让图更亮
subplot(234);imshow(log(I_margian+1),[]);title('原图幅值谱');axis on;
subplot(235);imshow(log(I1_margian+1),[]);title('缩放0.5倍后幅值谱');axis on;
subplot(236);imshow(log(I2_margian+1),[]);title('放大5倍后幅值谱');axis on;

运行结果如下：

由运行结果可知，随着图像的缩放，图像傅里叶变换后的幅值谱也相应的缩放

（3）旋转特性：如果f(x,y)旋转了一个角度，那么对应的傅立叶变换也旋转了相同的角度。

%% 旋转特性：
%% 如果f(x,y)旋转了一个角度，那么对应的傅立叶变换也旋转了相同的角度。
clc;clear;close all;    % 初始化
I=imread('C:\Users\xcz\Desktop\matlab_example\Images\lena.bmp');
I=im2double(I);     % 将数据类型转换为double，数值范围改为[0,1]
I=rgb2gray(I);  %将真彩色图像转换为灰度图像
I1=imrotate(I,45);  %将图片逆时针旋转45° 
I2=imrotate(I,90);  %将图片逆时针旋转90° 

% 分别求三幅图的幅值谱
I_fft=fftshift(fft2(I));    %进行傅里叶变换，并且进行象限转换
I_margian=abs(I_fft);   % 幅值谱

I1_fft=fftshift(fft2(I1));   
I1_margian=abs(I1_fft);  

I2_fft=fftshift(fft2(I2));   
I2_margian=abs(I2_fft);   
subplot(231);imshow(I);title('原图');axis on;
subplot(232);imshow(I1);title('逆时针旋转45°');axis on;
subplot(233);imshow(I2);title('逆时针旋转90°');axis on;
 % 参数'[]'是为了将其值线性拉伸,log是为了更好显示图，加1是为了让图更亮
subplot(234);imshow(log(I_margian+1),[]);title('原图幅值谱');axis on;
subplot(235);imshow(log(I1_margian+1),[]);title('旋转45°后幅值谱');axis on;
subplot(236);imshow(log(I2_margian+1),[]);title('旋转90°后幅值谱');axis on;

运行结果如下：

由结果可知，随着图像的旋转，其傅里叶变换结果也旋转相同的角度。

3. 图像的幅值谱和相位谱及双谱重构实验

傅里叶变换的频谱一般分为幅值谱和相位谱，频谱是一个以频率为自变量的函数。频谱在每一个频率点的取值是一个复数。一个复数由模和辐角唯一地确定，所以可将频谱分解为幅度谱（即复数的模关于频率的函数）和相位谱（即复数的辐角关于频率的函数）。一般将频谱取模即是幅值谱

利用matlab求图像的频谱图、幅值谱和相位谱。并根据相位谱、幅值谱重建图像。

%% 利用matlab求图像的频谱图、幅值谱和相位谱
clc;clear;close all;    % 初始化
I=imread('cameraman.tif');
I=im2double(I);
I1=fft2(I); % 傅里叶变换
I2=fftshift(I1);    %象限转换
r=real(I2); i=imag(I2);     % 求图像频域的实部和虚部
margin=log(abs(I2+1));      %图像幅度谱，加log便于显示
%phase=log(angle(I2)*180/pi);     %图像相位谱
phase=angle(I2);
l=log(I2+1);   % 图像频谱 加1是为了让频谱图更亮
p1_r=ifft2(abs(I1));    %幅值谱重构 使用未进行象限转换的幅值谱
p2_r=ifft2(exp(1i*(angle(I1)))); %相位谱重构 使用未进行象限转换的相位谱
% 此处1i表示复数
subplot(231);imshow(I);title('原图');
subplot(232);imshow(l,[]);title('频谱图');
subplot(233);imshow(margin,[]);title('幅值谱');
subplot(234);imshow(phase,[]);title('相位谱');
subplot(235);imshow(p1_r,[]);title('幅值谱重构');
subplot(236);imshow(p2_r,[]);title('相位谱重构');

运行结果如下：

由结果可知图像的幅度谱代表的是图像各像素点的亮度信息，即该像素应该显示什么颜色，但是做出来的幅度谱却不知道每一点在原图像中具体是哪一点，即幅度谱虽然存储了各个像素点的幅值信息，但是原像素点的位置已经被打乱，所以仅凭幅度谱是没有办法重构原图像的。幅度谱的中心是低频部分，越亮的地方代表的幅度越大。幅度谱中“十”字形亮线表示原图像中水平和垂直方向的分量较其他方向要多，因为在人们周围的自然场景中水平和垂直的线条出现的可能性较大。

而==相位谱记录的是所有点的相位信息，==看起来相位谱是一团噪声，这也说明相位信息是以一种更为隐蔽的方式出现，但它非常重要，因为相位信息中携带者图像的位置信息，没有它将无法从频谱还原出原图像。

可以看到仅有相位谱重构出来的图像还是可以看到一些轮廓信息的，利用相位谱记录的位置信息和幅度谱记录的亮度信息，就可以用双谱重构的方法恢复出原图像。

当使用双谱重构时：

%% 双谱重构关键代码
p3_r = ifft2(abs(I1).*exp(1i*(angle(I1))));%双谱重构
figure(1);imshow(p3_r,[]);title('双谱重构');

运行结果如下：

4. 图像的相位谱和幅值谱重组实验

对两幅图像分别求其幅值谱和相位谱，将相位谱交换，分别进行图像重构。

matlab代码如下：

%% 对两幅图像分别求其幅值谱和相位谱，将相位谱交换，分别进行图像重构。
clc;clear;close all;    %初始化
I1=imread('baby.jpg');
I2=imread('car2.jpg');
% 将彩色图像转换为灰度图片并且调整为相同大小
I1=rgb2gray(I1);I2=rgb2gray(I2);
I2=imresize(I2,[300 300]);I1=imresize(I1,[300 300]);
% 求傅里叶变换
I1_fft=fft2(I1); I2_fft=fft2(I2);
% 分别求幅值谱和相位谱
I1_margain=abs(I1_fft); I1_phase=angle(I1_fft);
I2_margain=abs(I2_fft); I2_phase=angle(I2_fft);
% 交换相位谱并重建复数矩阵
I1_re=I1_margain.*cos(I2_phase)+I1_margain.*sin(I2_phase).*1i;
I2_re=I2_margain.*cos(I1_phase)+I2_margain.*sin(I1_phase).*1i;
% 傅里叶反变换
I1_re1=abs(ifft2(I1_re)); I2_re1=abs(ifft2(I2_re)); % 此处取模是为了取复试的实部
%% 法二 效果相同
I1_re2 = ifft2(I1_margain.*exp(1i*(I2_phase)));%双谱重构
I2_re2 = ifft2(I2_margain.*exp(1i*(I1_phase)));%双谱重构
% 显示图像
subplot(221);imshow(I1);title('男孩原图');
subplot(222);imshow(I2);title('汽车原图');
subplot(223);imshow(I1_re1,[]);title('男孩的幅值谱和汽车的相位谱组合');
subplot(224);imshow(I2_re1,[]);title('汽车的幅值谱和男孩的相位谱组合');

运行结果如下：

五、实验总结

经过此次的实验，

• 明白了傅里叶变换在matlab的具体实现，进一步理解了傅里叶变换的相关性质，如平移性质、尺度变换性质和旋转性质等；
• 理解了相位谱、幅值谱的相关意义;
• 掌握了几个傅里叶变换的函数，初步知道了空间域和频率域之间的区别和联系，以及图像的频率域的变化规律；

参考资料

图像傅里叶变换的幅度谱、相位谱以及双谱重构原图像_逸凌Time的博客-程序员宅基地_图像的幅度谱和相位谱

数字图像处理：实验二 数字图像的傅里叶变换_nochengzi的博客-程序员宅基地