1. 图像傅里叶变换的物理意义
图1 图像的傅里叶变换
图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。如：大面积的沙漠在图像中是一片灰度变化缓慢的区域，对应的频率值很低；而对于地表属性变换剧烈的边缘区域在图像中是一片灰度变化剧烈的区域，对应的频率值较高。傅里叶变换在实际中有非常明显的物理意义，设f是一个能量有限的模拟信号，则其傅里叶变换就表示f的频谱。从纯粹的数学意义上看，傅里叶变换是将一个函数转换为一系列周期函数来处理的。从物理效果看，傅里叶变换是将图像从空间域转换到频率域，其逆变换是将图像从频率域转换到空间域。换句话说，傅里叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数。
傅里叶逆变换是将图像的频率分布函数变换为灰度分布函数傅里叶变换以前，图像（未压缩的位图）是由对在连续空间（现实空间）上的采样得到一系列点的集合，通常用一个二维矩阵表示空间上各点，记为z=f(x,y)。又因空间是三维的，图像是二维的，因此空间中物体在另一个维度上的关系就必须由梯度来表示，这样我们才能通过观察图像得知物体在三维空间中的对应关系。
傅里叶频谱图上我们看到的明暗不一的亮点，其意义是指图像上某一点与邻域点差异的强弱，即梯度的大小，也即该点的频率的大小（可以这么理解，图像中的低频部分指低梯度的点，高频部分相反）。一般来讲，梯度大则该点的亮度强，否则该点亮度弱。这样通过观察傅里叶变换后的频谱图，也叫功率图，我们就可以直观地看出图像的能量分布：如果频谱图中暗的点数更多，那么实际图像是比较柔和的（因为各点与邻域差异都不大，梯度相对较小）；反之，如果频谱图中亮的点数多，那么实际图像一定是尖锐的、边界分明且边界两边像素差异较大的。
对频谱移频到原点以后，可以看出图像的频率分布是以原点为圆心，对称分布的。将频谱移频到圆心除了可以清晰地看出图像频率分布以外，还有一个好处，它可以分离出有周期性规律的干扰信号，比如正弦干扰。一幅频谱图如果带有正弦干扰，移频到原点上就可以看出，除了中心以外还存在以另一点为中心、对称分布的亮点集合，这个集合就是干扰噪音产生的。这时可以很直观的通过在该位置放置带阻滤波器消除干扰。
1 %CSDN:by J27 copyright! 2 I = imresize(im2double(imread('cameraman.tif')),2); 3 figure; 4 imshowpair(I,log(abs(fftshift(fft2(I)))+1),'montage');
2. 图像的傅里叶变换matlab实现（幅度谱及相位谱）
img=imread('cameraman.tif'); %img=double(img); f=fft2(img); %傅里叶变换 f=fftshift(f); %使图像对称 r=real(f); %图像频域实部 i=imag(f); %图像频域虚部 margin=log(abs(f)); %图像幅度谱，加log便于显示 phase=log(angle(f)*180/pi); %图像相位谱 l=log(f); subplot(2,2,1),imshow(img),title('源图像'); subplot(2,2,2),imshow(l,[]),title('图像频谱'); subplot(2,2,3),imshow(margin,[]),title('图像幅度谱'); subplot(2,2,4),imshow(phase,[]),title('图像相位谱');
图1 图像的频域图、幅度谱及相位谱
3. 图像的傅里叶变换的频谱特征 一（周期性，能量分布，fftshift，交错性）
3.1 周期性
DFT的周期性：对于DFT而言，空域和频域都沿着X和Y方向无限周期拓展的；
图2-1 DFT的周期性
Matlab代码：
%CSDN:by J27 copyright! % Periodic I = imresize(im2double(imread('cameraman.tif')),2); Iperiodic = repmat(I,3,3); figure; imshowpair(Iperiodic,repmat((log(abs(fftshift(fft2(I)))+1)), 3, 3),'montage');
如果只取其中的一个周期，则我们会得到如下的结果（即，频谱未中心化）。
图2-2 频谱未居中
为了便于频域的滤波和频谱的分析，常常在变换之前进行频谱的中心化。
频谱的中心化：从数学上说是在变换之前用指数项乘以原始函数，又因为e^jπ = 1，所以实际上是把原始矩阵乘以（-1）^(x+y)达到频谱居中的目的。如下图所示：1<----->3 对调，2<----->4 对调，matlab中的fftshit命令就是这样实现的。
图2-3 对调频谱的四个象限（中心化）
Matlab代码：
%CSDN:by J27 copyright! I = imresize(im2double(imread('cameraman.tif')),2); figure; imshowpair(log(abs((fft2(I)))+1),log(abs(fftshift(fft2(I)))+1),'montage');
经过中心化后的频谱
截取了其中的一个周期，作为图像的频谱：
3.2 高低频率的分布
除了周期性之外，还应该知道的就是哪里是高频哪里是低频。在经过频谱居中后的频谱中，中间最亮的点是最低频率，属于直流分量（DC分量）。越往边外走，频率越高。所以，频谱图中的四个角和X,Y轴的尽头都是高频。
没有经过频谱居中处理的频谱图则正好相反，中间区域是高频，而四个角则是DC低频分量。
如下图是正弦波的例子来展示频谱图的高低频的分布。(下图中，中高频的图像出现了混叠)
频谱中心化以后，正弦波的频点靠中心越近，频率越低，离中心越远，频率越高。
Matlab代码：
1 %CSDN:by J27 copyright! 2 % Length of signal 3 L = 512; 4 % Sampling frequency 5 FsLow = 150; FsMid = 500; FsHigh = 1000; 6 % Form sampling vectors 7 IndexLow = linspace(0,FsLow,L); 8 IndexMid = linspace(0,FsMid,L); 9 IndexHigh = linspace(0,FsHigh,L); 10 11 % build 1D sinewave 12 SineLow = sin(IndexLow); 13 SineMid = sin(IndexMid); 14 SineHigh = sin(IndexHigh); 15 16 % translate 1D wave into 2D sinewave 17 SinewaveL = repmat(SineLow,[L,1]); 18 SinewaveM = repmat(SineMid,[L,1]); 19 SinewaveH = repmat(SineHigh,[L,1]); 20 21 % Window function 22 Window = hanning(L)*hanning(L)'; 23 SinewaveLwindowed = SinewaveL.*Window; 24 SinewaveMwindowed = SinewaveM.*Window; 25 SinewaveHwindowed = SinewaveH.*Window; 26 figure; 27 subplot(2,3,1),imshow(SinewaveL,[]),title('Low-freq. sinewave'), 28 subplot(2,3,2),imshow(SinewaveM,[]),title('Med-freq. sinewave'), 29 subplot(2,3,3),imshow(SinewaveH,[]),title('Hig-freq. sinewave'), 30 subplot(2,3,4),imshow(log(abs(fftshift(fft2(SinewaveLwindowed)))+1),[]),title('Spectrum of Low-freq. sinewave'), 31 subplot(2,3,5),imshow(log(abs(fftshift(fft2(SinewaveMwindowed)))+1),[]),title('Spectrum of Med-freq. sinewave'), 32 subplot(2,3,6),imshow(log(abs(fftshift(fft2(SinewaveHwindowed)))+1),[]),title('Spectrum of Hig-freq. sinewave');
3.3 频谱图的能量分布
频谱中的能级分布，则如下图所示：DC分量所占能量最大最多，不论是二维还是一维都应该是这样。频率越高的部分，能量越少。如下图所示，图示画的不好，勉强能够理解就好。中间最小的那个圆圈内包含了大约85%的能量，中间那个圈包含了大约93%的能量，而最外面那个圈则包含了几乎99%的能量。
3.4 方向一致性
在二维傅里叶变换中，空间域中横向的周期变化会反应在频谱图中的横轴上，而空间域中纵向的周期变化会反应在频谱图中的纵轴上。空间域中东南方向的周期变化会反应在频谱图中的东北方向，反之亦然。说明见下图。
Matlab代码：
%CSDN:by J27 copyright! % black&white Isize = 512; Iblackwhite = zeros(Isize); Iblackwhite(1:Isize/2,:) = 1; figure; imshowpair(Iblackwhite,log(abs(fftshift(fft2(Iblackwhite)))+1),'montage'); Iwhiteblack = zeros(Isize); Iwhiteblack(:,Isize/2:end) = 1; figure; imshowpair(Iwhiteblack,log(abs(fftshift(fft2(Iwhiteblack)))+1),'montage'); % black&white triangle Hanning = hanning(Isize)*hanning(Isize)'; ItriangleUpper = triu(ones(Isize),0); ItriangleUpper = ItriangleUpper.*Hanning; figure; imshowpair(ItriangleUpper,log(abs(fftshift(fft2(ItriangleUpper)))+1),'montage'); % ItriangleLower = tril(ones(Isize),0); % ItriangleLower = ItriangleLower.*Hanning; % figure; % imshowpair(ItriangleLower,log(abs(fftshift(fft2(ItriangleLower)))+1),'montage'); % flip the upper triangle matrix ItriangleUpperMirror = fliplr(triu(ones(Isize),0)); ItriangleUpperMirror = ItriangleUpperMirror.*Hanning; figure; imshowpair(ItriangleUpperMirror,log(abs(fftshift(fft2(ItriangleUpperMirror)))+1),'montage');
正弦波灰度值分别在横向、纵向变化的例子.(可对比一维正弦信号的频域分析来理解，如下图)
Matlab代码：
1 %CSDN:by J27 copyright! 2 % form horizontal and vertical sinewaves 3 SineHorz = SinewaveM; 4 SineVert = rot90(SineHorz,1); 5 SineHorzwindowed = SineHorz.*Window; 6 SineVertwindowed = SineVert.*Window; 7 8 figure; 9 subplot(2,2,1),imshow(SineHorz,[]),title('Vertical sinewave'), 10 subplot(2,2,2),imshow(SineVert,[]),title('Horizontal sinewave'), 11 subplot(2,2,3),imshow(log(abs(fftshift(fft2(SineHorzwindowed)))+1),'InitialMagnification','fit'),title('Spectrum of Vertical sinewave'), 12 subplot(2,2,4),imshow(log(abs(fftshift(fft2(SineVertwindowed)))+1),'InitialMagnification','fit'),title('Spectrum of Horizontal sinewave'),
参考原文：https://blog.csdn.net/daduzimama/article/details/80109139
Matlab代码：
1 %**************************************************************************************** 2 % 正弦信号及频域分析 3 %*************************************************************************************** 4 % 5 fs = 1000; % 采样频率1K 6 N = 256 % 采样点数 7 n = 0:N-1; 8 A1 = 100; % 幅值 9 W1 = 50; % 频率 10 P1 = 45; % 初始相位 11 t = n/fs; % 时间序列 12 f = n*fs/N; % 频率序列 13 14 x = A1*sin(2*pi*W1*t + P1/180*pi); %50Hz和200Hz正弦波合成 15 16 subplot(211); 17 plot(t, x); %绘制信号Mix_Signal的波形 18 xlabel('时间'); 19 ylabel('幅值'); 20 title('正弦信号'); 21 grid on; 22 23 subplot(212); 24 y = fft(x, N); %对信号 Mix_Signal做FFT 25 plot(f, abs(y)); 26 xlabel('频率/Hz'); 27 ylabel("振幅"); 28 title('正弦信号的频域图'); 29 grid on;
4. 图像的噪音与频域滤波
4.1 添加噪声
读入一幅图像，分别为图像添加椒盐、高斯噪声，做FFT变换，并显示其原图及对应的频谱图。
椒盐噪声：对原图像模拟叠加密度为0.04的椒盐噪声；
高斯噪声：模拟均值为0方差为0.02的高斯噪声。
Matlab代码：
1 I=imread('D:\1\lena512.bmp'); %读入图像 2 F=fft2(im2double(I)); %FFT 3 F=fftshift(F); %FFT频谱平移 4 F=real(F); 5 T=log(F+1); %频谱对数变换 6 subplot(3,2,1),imshow(I),title('原始图像'); 7 subplot(3,2,2),imshow(T,[]),title('原始图像其频谱图'); 8 9 S=imnoise(I,'salt & pepper', 0.04); %模拟叠加密度为0.04的椒盐噪声 10 K=fft2(im2double(S)); %FFT 11 K=fftshift(K); %FFT频谱平移 12 K=real(K); 13 T=log(K+1); %频谱对数变换 14 subplot(3,2,3),imshow(S),title('添加椒盐噪声后的图像'); 15 subplot(3,2,4),imshow(T,[]),title('椒盐噪声频谱图'); 16 17 G=imnoise(I,'gaussian',0,0.02);%模拟均值为0方差为0.02的高斯噪声， 18 H=fft2(im2double(G)); %FFT 19 H=fftshift(H); %FFT频谱平移 20 H=real(H); 21 T=log(H+1); %频谱对数变换 22 subplot(3,2,5),imshow(G),title('添加高斯噪声后的图像'); 23 subplot(3,2,6),imshow(T,[]),title('高斯噪声频谱图');
注：lena测试图像的下载链接：https://www.ece.rice.edu/~wakin/images/
4.2 图像频域滤波
读入一幅图像，对图像分别进行高斯低通、巴特沃兹低通、高斯高通和巴特沃兹高通频域滤波，比较其锐化和平滑效果；
对于高斯低通频域滤波，首先求原图像的频谱图，然后根据二维高斯低通滤波器(GLPF)定义，对频谱图做高斯低通滤波，使低频通过而使高频衰减，最后做快速傅里叶逆变换，结果发现滤波后的图像变模糊，比原始图像减少尖锐的细节部分而突出平滑过渡部分；对于巴特沃兹低通频域滤波，然后根据二级巴特沃思低通滤波器(BLPF)定义，对频谱图做巴特沃兹低通滤波，使低频通过而使高频衰减，最后做快速傅里叶逆变换，结果发现滤波后的图像变模糊，比原始图像减少尖锐的细节部分而突出平滑过渡部分；经对比图像后发现，经过巴特沃思低通滤波的图像比经过高斯低通频域滤波的图像更平滑。
对于高斯高通频域滤波，首先求原图像的频谱图，然后根据截频距原点为D0的高斯高通滤波器(GHPF)定义，对频谱图做高斯高通滤波，使高频通过而使低频衰减，最后做快速傅里叶逆变换，结果发现滤波后的图像变锐化，比原始图像减少平滑过渡而突出边缘等细节部分；对于巴特沃兹高通频域滤波，然后根据二阶且截至频率距原点的距离为D0的巴特沃思高通滤波器(BHPF)定义，对频谱图做巴特沃兹高通滤波，使高频通过而使低频衰减，最后做快速傅里叶逆变换，结果发现滤波后的图像变锐化，比原始图像减少平滑过渡而突出边缘等细节部分；经对比图像后发现，经过高斯低通滤波的图像比经过高斯低通巴特沃思低通频域滤波的图像更平滑。
高斯低通频域滤波Matlab代码：
1 %高斯低通频域滤波 2 I=imread('D:\1\lena512.bmp'); %读入图像 3 subplot(1,2,1),imshow(I),title('原始图像'); 4 I=double(I); 5 S=fftshift(fft2(I)); 6 [M,N]=size(S); 7 n=2; 8 d0=30; %GLPF滤波，d0=5，15，30（程序中以d0=30为例） 9 n1=floor(M/2); 10 n2=floor(N/2); 11 for i=1:M 12 for j=1:N 13 d=sqrt((i-n1)^2+(j-n2)^2); 14 h=1*exp(-1/2*(d^2/d0^2)); 15 S(i,j)=h*S(i,j); 16 end 17 end 18 S=ifftshift(S); 19 S=uint8(real(ifft2(S))); 20 subplot(1,2,2),imshow(S),title('高斯低通滤波图像');
巴特沃斯低通频域滤波Matlab代码：
%巴特沃斯低通频域滤波 I=imread('D:\1\lena512.bmp'); %读入图像 subplot(1,2,1),imshow(I),title('原始图像'); F=double(I); % 数据类型转换，MATLAB不支持图像的无符号整型的计算 G=fft2(F); % 傅立叶变换 G=fftshift(G); % 转换数据矩阵 [M,N]=size(G); nn=2; % 二阶巴特沃斯(Butterworth)低通滤波器 d0=30; %截止频率为30 m=fix(M/2); n=fix(N/2); for i=1:M for j=1:N d=sqrt((i-m)^2+(j-n)^2); h=1/(1+0.414*(d/d0)^(2*nn)); % 计算低通滤波器传递函数 result(i,j)=h*G(i,j); end end result=ifftshift(result); Y2=ifft2(result); Y3=uint8(real(Y2)); subplot(1,2,2),imshow(Y3),title('巴特沃斯低通滤波')
巴特沃斯高通频域滤波Matlab代码：
%巴特沃斯高通频域滤波 I=imread('D:\1\lena512.bmp'); %读入图像 subplot(121),imshow(I); title('原始图像'); F=double(I);% 数据类型转换，MATLAB不支持图像的无符号整型的计算 G=fft2(F);% 傅立叶变换 G=fftshift(G);%转换数据矩阵 [M,N]=size(G); nn=2;% 二阶巴特沃斯(Butterworth)高通滤波器 d0=30; m=fix(M/2);n=fix(N/2); for i=1:M for j=1:N d=sqrt((i-m)^2+(j-n)^2); if (d==0) h=0; else h=1/(1+0.414*(d0/d)^(2*nn));% 计算传递函数 end result(i,j)=h*G(i,j); end end result=ifftshift(result); J2=ifft2(result); J3=uint8(real(J2)); subplot(122),imshow(J3); title('巴特沃斯高通滤波后图像'); % 滤波后图像显示