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('巴特沃斯高通濾波后圖像'); % 濾波后圖像顯示