六、数字图像处理之频域处理（1）

（一）在MATLAB中计算及观察二维DFT

常用的DFT函数总结：

F = fft2(f)           %获得一幅图像的傅里叶变换
F = fft2(f,P,Q)       %对图像填充所需数目的0，结果大小变为P×Q
S = abs(F)            %获得傅里叶谱
Fc = fftshift(F)      %将变换原点移动到频域矩形的中心
S2 = log(1+abs(Fc))   %用log变换增强的频谱
F = ifftshift(Fc)     %反居中变换函数
phi = atan2(I,R)      %表示相位
例如：phi = atan2(1,1)  %表示45°  
phi = atan2(1,-1)      %表示135° 
f = ifft2(F)          %傅里叶逆变换

实验代码：

f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0303(a).tif');
subplot(331),imshow(f), title('原图像图像');
F = fft2(f);  %获得一幅图像的傅里叶变换
subplot(332),imshow(uint8(real(F))), title('傅里叶变换图像');
Fc = fftshift(F);      %将变换原点移动到频域矩形的中心
subplot(333),imshow(Fc,[]), title('变换移到频谱中心的图像');
S = abs(F);            %获得傅里叶谱
subplot(334),imshow(S,[]), title('傅里叶频谱图像');
phi = angle(F);
F = S.*exp(i*phi);
subplot(335),imshow(F,[]), title('傅里叶相位谱图像');
S2 = log(abs(Fc));   %用log变换增强的频谱
subplot(336),imshow(S2,[]), title('变用log变换增强的频谱的图像');
F = ifftshift(Fc);     %反居中变换函数
subplot(337),imshow(uint8(real(F))), title('反居中变换函数的图像');
f = ifft2(F);          %傅里叶逆变换
subplot(338),imshow(f,[]), title('傅里叶逆变换的图像');

实验效果：

实验说明：
仔细一下，其实右上图中间有一个亮点，中左和正中间图像的四个角有亮点，这四个亮点集中在右上图的图像上（即将傅里叶变换移动到频谱中心），中右图显示的是右上图用Log增强以后的效果图。

（二）频域滤波

在这里，符号"星号"表示两个函数的卷积，双箭头两边的表达式组成傅立叶变换对。例如, 第一个表达式表明两个空间函数(表达式左侧的项)的卷积可以通过计算两个傅立叶变换函数(表达式的右侧)的乘积的反变换得到。相反，两个空间函数的卷积的正傅立叶变换给出了两个函数傅立叶变换的乘积。同样的解释也出现在第二个表达式中。

若使用DFT进行滤波操作，则图像及其变换都将自动地看成是周期性的。若周期关于函数的非零部分的持续时间很靠近，则对周期函数执行卷积运算会导致相邻周期之间的干扰，称之为折叠误差的干扰。可通过使用零来填充函数的方法避免。若处理的函数大小均是M×N，则填充值可设置为：P≥2M-1 和 Q ≥2N-1。

paddedsize函数用于计算满足P和Q的最小偶数值，函数也提供填充输入图像，从而形成的方形图像的大小等于最接近的2的整数次幂。

实验代码：

参数 paddedsize 函数：

function PQ = paddedsize(AB, CD, PARAM)
%PADDEDSIZE Computes padded sizes useful for FFT-based filtering. 
%   PQ = PADDEDSIZE(AB), where AB is a two-element size vector,
%   computes the two-element size vector PQ = 2*AB.
%
%   PQ = PADDEDSIZE(AB, 'PWR2') computes the vector PQ such that
%   PQ(1) = PQ(2) = 2^nextpow2(2*m), where m is MAX(AB).
%
%   PQ = PADDEDSIZE(AB, CD), where AB and CD are two-element size
%   vectors, computes the two-element size vector PQ.  The elements
%   of PQ are the smallest even integers greater than or equal to 
%   AB + CD - 1.
%
%   PQ = PADDEDSIZE(AB, CD, 'PWR2') computes the vector PQ such that
%   PQ(1) = PQ(2) = 2^nextpow2(2*m), where m is MAX([AB CD]). 

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.5 $  $Date: 2003/08/25 14:28:22 $

if nargin == 1
   PQ  = 2*AB;
elseif nargin == 2 & ~ischar(CD)
   PQ = AB + CD - 1;
   PQ = 2 * ceil(PQ / 2);
elseif nargin == 2
   m = max(AB); % Maximum dimension.

   % Find power-of-2 at least twice m.
   P = 2^nextpow2(2*m);
   PQ = [P, P];
elseif nargin == 3
   m = max([AB CD]); % Maximum dimension.
   P = 2^nextpow2(2*m);
   PQ = [P, P];
else 
   error('Wrong number of inputs.')
end

低通滤波器的传递函数lpfilter函数：

function H = lpfilter(type, M, N, D0, n)
[U, V] = dftuv(M, N);
D = hypot(U,V);
switch type
case 'ideal'
   H = double(D <= D0);
case 'btw'
   if nargin == 4
      n = 1; 
   end
   H = 1./(1 + (D./D0).^(2*n));
case 'gaussian'
   H = exp(-(D.^2)./(2*(D0^2)));
otherwise
   error('Unknown filter type.')
end

dftuv函数：

function [U,V] = dftuv(M,N)
u = double(0:(M-1));
v = double(0:(N-1));
idx = find(u >M/2);
u(idx) = u(idx) - M;
idy = find(v >N/2)
v(idy) = v(idy) - N;

[V,U] = meshgrid(v,u);
[M,N] = size(f);
F = fft2(f);
sig = 10;
H = lpfilter('gaussian',M,N,sig);
G = H.*F;
g = ifft2(G);
g = revertclass(g);
imshow(g);

end

（ dftfilt函数的代码见下面）

f1 = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0309(a).tif');
%未使用填充的频率低通滤波处理
[M,N] = size(f1);          %设置维数
F = fft2(f1);                  %计算傅里叶变换
sig = 10;                   % 指定低通滤波器的标准偏差
H = lpfilter('gaussian',M,N,sig);     %高斯低通滤波器
G = H.*F;                         %将变换乘以滤波函数
g =real( ifft2(G));              %获得傅里叶逆变换的实部
subplot(131), imshow(f1), title('原图像');
subplot(132), imshow(g,[]), title('未使用填充的频率低通滤波');

%使用填充的频率低通滤波处理
f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0217(a).tif');
PQ = paddedsize(size(f));                   %获取填充参数
Fp = fft2(f1,PQ(1),PQ(2));                   %得到使用填充的傅里叶变换
Hp = lpfilter('gaussian',PQ(1),PQ(2),2*sig);    %高斯低通滤波器
Gp = Hp.*Fp;                                %将变换乘以滤波函数
gp = real(ifft2(Gp));                       %获得傅里叶逆变换的实部
gpc = gp(1:size(f1,1),1:size(f,2));          %将左上部的矩形修剪为原始大小
subplot(133);imshow(gp,[ ]),title('使用填充的低通滤波');

实验结果：

（三）DFT滤波的基本步骤

（1）用函数tofloat把输入图像变换为浮点图像

[f,revertclass] = tofloat(f);

（2）用函数paddedsize获得填充参数

PQ = paddedsize(size(f));

（3）得到有填充的傅里叶变换

F = fft2(f,PQ(1),PQ(2));

（4）生成大小为PQ（1）×PQ（2）的滤波函数H，在这一步可以使用本节中提到的任何算法。滤波函数H必须是DFT滤波格式。如果是居中的，在使用滤波器之前令H=fftshift(H)。
（5）用滤波器乘以FFT变换

G = H.^F;

（6）获得G的逆FFT变换

g = ifft2(G);

（7）修剪左上部矩形为原始大小

g = g(1:size(f,1),size(f,2));

（8）把滤波过的图像变换为输入图像的类，如果希望的话

g = reverttclass(g);

（四）频域滤波的M-函数

function g = dftfilt(f, H, classout) 
[f, revertClass] = tofloat(f);               %把输入图像变换为浮点图像：
F = fft2(f, size(H, 1), size(H, 2));         %获得有填充的傅立叶变换
g = ifft2(H.*F);                             %获得H.*F的逆FFT变换
g = g(1:size(f,1), 1:size(f, 2));           %修剪图像为原始大小
if nargin == 2 || strcmp(classout,'priginal')        %将输出图像转化为与输入图像相同的类
    g = revertClass(g);                  
 elseif strcmp(classout,'fltpoint')         
   return 
 else
   error('Undefined class for the output image.’)
end

（五）空域滤波和频域滤波器的比较

说明：通常，当滤波器较小时，使用空域滤波要比频域滤波在计算上更有效。在图像处理工具箱中函数freqz2可以输出相应的频域滤波器和计算FIR滤波器的频率响应。
语法表示为：

H = freqz2(h,R,C); %h是二维空域滤波器，H是相应的二维频域滤波器,R是行数，C是列数

实验代码如下：

f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0309(a).tif');
F = fft2(f);
S = fftshift(log(1+abs(F)));
h = fspecial('sobel');
PQ = paddedsize(size(F));
H = freqz2(h,PQ(1),PQ(2));
H1 = ifftshift(H);
gs = imfilter(double(f),h);
gf = dftfilt(f,H1);
subplot(231);imshow(f);title('灰度级图像');
subplot(232);imshow(S,[]);title('原图像傅里叶谱');
subplot(233);imshow(abs(H),[ ]);title('对应垂直sobel空域滤波器的频域滤波 ');
subplot(234);imshow(abs(H1),[ ]);title('用ifftshift处理后的滤波器');
subplot(235);imshow(gs,[ ]);title('以图形显示的空域滤波器');
subplot(236);imshow(gf,[ ]);title('用ifftshift处理后的滤波器的图像');

实验结果：