数字图像处理与MATLAB 第三章学习笔记

第三章 频率域滤波

本章大部分内容与第二章讨论的滤波主题并行，只是所有滤波都是通过傅里叶变换在频率域中实现。

3.1 二维离散傅里叶变换

令f(x, y)表示一幅大小为M×N像素的数字图像，x∈[0,M-1]，y∈[0,N-1]，由F(u,v)表示的f(x, y)的二维离散傅里叶变换(DFT)由下式给出：$F(u,v)={\sum }_{x=0}^{M=1}{\sum }_{y=0}^{N=1}f(x,y)e^{-j2\pi (\frac{ux}{M}+\frac{vy}{N})}$，式中u∈[0,M-1]，v∈[0,N-1]。使用确定频率的变量u、v可将指数项展开为正弦和余弦函数。由u∈[0,M-1]、v∈[0,N-1]定义的大小为M×N的矩形区域通常称为频率矩形。

离散傅里叶反变换(IDFT)形式为$f(x,y)=\frac{1}{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,M-1]，y∈[0,N-1]；因此，给定F(u,v)即可借助IDFT得到f(x,y)。该式中F(u,v)值有时称为该展开式的傅里叶系数。

频率域原点处变换的值（如F(0,0)）称为傅里叶变换的直流分量，F(0,0)等于f(x, y)平均值的MN倍。

直观的分析一个变换的主要方法为计算它的频谱（即F(u, v)的幅度，是一个实函数），并将其显示为一幅图像。令R(u, v)和I(u, v)分别表示F(u, v)的实部和虚部，则傅里叶谱定义为$|F(u,v)|=[R^2(u,v)+I^2(u,V){]}^{\frac{1}{2}}$，变换的相角定义为$\varphi (u,v)=arctan[\frac{I(u,v)}{R(u,v)}]$。这两个函数可在极坐标形式下表示复函数F(u, v)：$F(u,v)=|F(u,v)|e^{-j\varphi (u,v)}$。

功率谱定义为幅度的平方：$P(u,v)=|F(u,v){|}^2=R^2(u,v)+I^2(u,v)$。为直观起见，用|F(u, v)|还是P(u, v)并不重要。

若f(x, y)是实函数，则其傅里叶变换关于原点共轭对称，即$F(u,v)=F^{\ast }(-u,-v)$；这意味着傅里叶谱也关于原点对称：$|F(u,v)|=|F(-u,-v)|$。直接代入F(u, v)的公式中即可证明$F(u,v)=F(u+k_{1}M,v)=F(u,v+k_{2}N)=F(u+k_{1}M,v+k_{2}N)$。式中$k_1$和$k_2$是整数，即DFT在u和v方向上是无穷周期的，周期由M和N决定。IDFT同理：$f(x,y)=f(x+k_{1}M,y)=f(x,y+k_{2}N)=f(x+k_{1}M,y+k_{2}N)$。

3.2 在MATLAB中计算和观察二维DFT

在实践中，DFT和IDFT可以用快速傅里叶变换(FFT)算法实现。使用函数fft2可以得到一个图像数组的FFT，语法为F=fft2(f)；此函数返回的傅里叶变换大小仍为M×N。使用傅里叶变换滤波时需要对数据进行零填充，此时语法变为F=fft2(f, P, Q)。

傅里叶谱可以使用函数abs获得：S=abs(F)，该函数计算数组中每个元素的幅度（实部和虚部平方和的平方根）。

使用函数fftshift将变换的原点移动到频率矩形的中心，语法Fc=fftshift(F)，其中F是fft2计算的变换，Fc是居中后变换的结果。函数fftshift是通过变换F的象限来运算的。

函数ifftshift使居中结果反转，语法F=iffshift(Fc)，也可以用来将对最初在矩形内居中的函数转换为中心点在矩形左上角的函数。

二维傅里叶变换的实部R(u, v)和虚部I(u, v)使与F(u, v)大小相同的数组，因为R和I均可正可负，所以需要能在整个区间$[-\pi ,\pi ]$内计算反正切（具有此性质的函数称为四象限反正切），使用函数atan2完成计算，语法phi=atan2(I, R)，其中phi是I和R大小相同的数组，phi的元素是在区间$[-\pi ,\pi ]$内以弧度表示的角度，这些角度是关于实轴度量的。函数*real(arg)和imag(arg)*分别提供arg的实部和虚部，*P=angle(z)*返回复数组z的每个元素的相角，单位弧度，值域$\pm \pi$。

给定谱及相应的相角，使用下列语句可以得到DFT：F=S.×exp(i×phi)。

谱的成分决定了通过组合形成一幅图像的正弦的幅度，相位携带着不同正弦关于其原点的位移信息。

通过使用函数floor，同时记住MATLAB原点为(1, 1)，即可得MATLAB计算所用的频率矩形中心为[floor(M/2)+1, floor(N/2)+1]，不管M、N值是奇是偶，该公式给出的中心点都是正确的。

函数ifft2用于计算傅里叶反变换，语法f=ifft2(F)，F是傅里叶变换，f是结果图像。fft2会把输入图像毫无缩放地变为double类。

3.3 频率域滤波

3.3.1 基础

空间域和频率域地线性滤波的基础都是卷积定理，可写为$f(x,y)★h(h,y)\iff H(u,v)F(u,v)$，或写为$f(x,y)h(h,y)\iff H(u,v)\ast G(u,v)$，其中符号★表示两个函数的卷积，双箭头两边的表达式组成一个傅里叶变换对。第一个表达式表明两个空间函数的卷积（表达式左边的项），可以通过计算这两个函数的傅里叶变换的乘积（右侧）的反变换得到。相反，两个空间函数的卷积的傅里叶变换给出了两个函数的傅里叶变换的乘积。

低通滤波器：当乘以一个居中处理后的函数F(u, v)后会衰减F(u, v)的高频分量，而低频分量相对不变的滤波器。会模糊（平滑）图像。

空间滤波是指滤波器模板h(x, y)与一幅图像f(x, y)卷积，函数相对移位，直到一个函数全部滑过另一个函数为止。根据卷积定理，在频率域中让F(u, v)乘以空间滤波器的傅里叶变换H(u, v)，可得到相同的结果；但在处理离散量的时候，由于F、H具有周期性，即在离散频率域中执行的卷积是周期的，因此用DFT执行的卷积称为循环卷积。保证空间和循环卷积给出相同结果的唯一方法是使用适当的零填充。

基于卷积定理，要在空间域中得到相应的滤波后的图像，就需要计算成绩H(u, v)和F(u, v)的傅里叶反变换。在周期接近函数非零部分的持续周期时，对周期函数进行卷积会引起相邻周期的串扰，称为折叠误差，可以通过补零法避免。

卷积定理的离散版本要求参与卷积运算的两个函数必须在空间域内进行填充，为了避免折叠误差。进行滤波时卷积处理涉及的两个函数中的一个是滤波器；但在频率域中使用DFT进行滤波处理时会在域中指定一个滤波器，其尺寸等于填充后的图像尺寸。

3.3.2 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=ifftshift(H)。

5.用滤波器乘以该变换：G=H.×F;

6.获得G的IFFT：g=ifft2(G);

7.将左上部举行裁剪为原始大小：g=g(1:size(f, 1), 1:size(f, 2));

8.需要时，将滤波后的图像转换为输入图像的类：g=revertclass(g);

由线性系统理论可知，在某种合适的条件下，向线性系统输入一个冲激动，可以完全表征该系统。若该系统是一个滤波器，则可以通过观察它对冲激的响应来完全确定该滤波器，以这种方式确定的滤波器称为有限冲激响应(FIR)滤波器。本书中所有线性滤波器都是FIR滤波器。

3.3.3 用于频率域滤波的M函数

使用函数dftfilt接受输入图像和一个滤波函数，处理掉所有滤波细节，输出滤波后的结果并裁剪图像。

function g = dftfilt(f, H, classout)
%DFTFILT Performs frequency domain filtering.
% g = DFTFILT(f, H, CLASSOUT) filters f in the frequency domain
% using the filter transfer function H. The output, g, is the
% filtered image, which has the same size as f.
%
% Valid values of CLASSOUT are
%
% 'original' The ouput is of the same class as the input.
% This is the default if CLASSOUT is not included
% in the call.
% 'fltpoint' The output is floating point of class single, unless
% both f and H are of class double, in which case the
% output also is of class double.
%
% DFTFILT automatically pads f to be the same size as H. Both f
% and H must be real. In addition, H must be an uncentered,
% circularly-symmetric filter function.

% Convert the input to floating point.
[f, revertClass] = tofloat(f);

% Obtain the FFT of the padded input.
F = fft2(f, size(H, 1), size(H, 2));

% Perform filtering.
g = ifft2(H.*F);

% Crop to original size.
g = g(1:size(f, 1), 1:size(f, 2)); % g is of class single here.

% Convert the output to the same class as the input if so specified.
if nargin == 2 || strcmp(classout, 'original')
    g = revertClass(g);
elseif strcmp(classout, 'fltpoint')
    return
else
    error('Undefined class for the output image.')
end

3.4 从空间滤波器获得频率域滤波器

通常，滤波器较小时在计算上空间滤波要比频率域滤波更有效。

函数freqz2计算FIR滤波器的频率响应，结果是所希望的频率域滤波器，语法形式为H=freqz2(h, R, C)，h是一个二维空间滤波器，H是相应的二维频率域滤波器，R是行数，C是我们希望滤波器H所具有的列数。

3.5 在频率域中直接生成滤波器

3.5.1 创建用于实现频率域滤波器的网格数组

M函数dftuv提供了距离计算及其他类似应用所需的网格数组，其结果是按fft2和ifft2处理的要求排列的，不需要重排数据。

函数hypot与语句*D=sqrt(U.^2+V.2)*等效，但速度更快。

3.5.2 低通（平滑）频率域滤波器

理想低通滤波器(ILPF)具有如下传递函数：$H(u,v)=\{\begin{array}{l}1,D(u,v)\le D_0\\ 0,D(u,v)>D_0\end{array}$，式中$D_0$为正数，$D(u,v)$为点(u, v)到滤波器中心的距离。满足$D(u,v)=D_0$的点的轨迹是一个圆。若滤波器H(u, v)乘以一幅图像的傅里叶变换，则可以观察到一个理想滤波器会切断该圆以外的所有F(u, v)分量，而保留圆上和圆内的所有分量不变。

n阶巴特沃斯低通滤波器(BLPF)在距滤波中心$D_0$处有截止频率，其传递函数为$H(u,v)=\frac{1}{1+[\frac{D(u,v)}{D_0}{]}^{2n}}$，当$D(u,v)=D_0$时，H(u, v)=0.5（降为最大值1的50%）。与ILPF不同的是，BLPF的传递函数在$D_0$点处不存在尖锐的不连续。对于具有平滑传递函数的滤波器，通常将截止频率轨迹定义在H(u, v)降低为其最大值的一个指定的比例点处。

高斯低通滤波器(GLPF)的传递函数由下式给出：$H(u,v)=e^{-D^2(u,v)/2{\sigma }^2}$，σ为标准差；通过令$\sigma =D_0$，根据截止参数可得$H(u,v)=e^{-D^2(u,v)/2D_0^2}$。当$D(u,v)=D_0$时，该滤波器降到其最大值1的60.7%。

表3.1 低通滤波器。[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-5kI3CI8s-1647150015442)(D:/daily/Matlab%E7%94%A8/%E6%95%B0%E5%AD%97%E5%9B%BE%E5%83%8F%E5%A4%84%E7%90%86%E4%B8%8EMATLAB%EF%BC%88%E7%AC%94%E8%AE%B0%EF%BC%89.assets/clip_image002-164592440497019-164592440632520.png)]是截止频率，n是巴特沃斯滤波器的阶。

理想	巴特沃斯	高斯
$H(u,v)=\{\begin{array}{l}1,D(u,v)\le D_0\\ 0,D(u,v)>D_0\end{array}$	$H(u,v)=\frac{1}{1+[\frac{D(u,v)}{D_0}{]}^{2n}}$	$H(u,v)=e^{-D^2(u,v)/2{\sigma }^2}$

函数lpfilter用于生成上表的几个低通滤波器的传递函数，也可用于生成高通滤波器。

3.5.3 绘制线框图和表面图

函数mesh可以绘制大小为M×N的二维函数，语法mesh(H)。该函数针对x=1:M和y=1:N绘出一幅线框图，M和N很大时通常会绘出密集的线框图；为避免如此可使用下列语法每隔k个点进行绘制：mesh(H(1:k:end, 1:k:end))。通常沿每个坐标轴40~60个点可在分辨率和外观之间得到较好的平衡。

函数colormap可以设置线框图的颜色，语法*colormap([0 0 0])*可将线框图设置为黑色。

MATLAB默认会将网格和坐标轴叠放到网格图上，使用命令grid off可关闭网格，命令axis off可以关闭坐标轴。语法view(az, el)可以控制观察点的位置，az是方位角，el是仰角，单位为度，默认值是-37.5和30。语句[az, el]=view可以确定当前的几何视图，语句*view(3)*可将观察点设置为默认值。

函数surf可以将函数绘制成表面图，基本语法surf(H)。该函数与mesh生成相同的图形，但surf生成的网格中四边形会被彩色填充（称为小面描影），使用命令colormap(gray)可将彩色转换为灰度；命令shading interp可平滑小面描影并去掉网格线。

当目标是绘制双变量的一个解析函数时，可使用meshgrid生成坐标值，坐标值可生成将在mesh或surf中使用的离散（采样）矩阵。

3.6 高通（锐化）频率域滤波器

高通滤波的方法是衰减傅里叶变换的低频部分而保持高频部分相对不变。

若给定低通滤波器的传递函数$H_{LP}(u,v)$，则相应的高通滤波器传输函数为$H_{HP}(u,v)=1-H_{LP}(u,v)$。下表显示了对应表3.1中低通滤波器的高通滤波器的传输函数。

表3.2 高通滤波器。$D_0$是截止频率，n是巴特沃斯滤波器的阶。

理想	巴特沃斯	高斯
$H(u,v)=\{\begin{array}{l}0,D(u,v)\le D_0\\ 1,D(u,v)>D_0\end{array}$	$H(u,v)=\frac{1}{1+[\frac{D_0}{D(u,v)}{]}^{2n}}$	$H(u,v)=1-e^{-D^2(u,v)/2D_0^2}$

3.6.1 一个用于高通滤波的函数

高通滤波器函数hpfilter可由低通滤波器函数lpfilter推导而来。使用高通滤波后图像中标远和其他灰度急剧过渡得到了增强，但由于图像平均值由F(0, 0)给出，且迄今为止讨论的高通滤波器偏离了傅里叶变换的原点，因此该图像失去了大部分原图像中呈现的灰色调。

3.6.2 高频强调滤波

高通滤波器偏离了直流项，因此会将图像的平均值降低为0，补偿方法之一是给高通滤波器加一个偏移量；若把偏移量与将滤波器乘以一个大于1的常数结合起来，则这种方法称为高频强调滤波（该常量乘数会突出高频部分、增大低频部分的幅度，但只要该偏移量与乘数相比较小，则低频增强的影响就弱于高频增强的影响）。高频强调滤波器的传递函数为$H_{HFE}(u,v)=a+b{H}_{HP}(u,v)$，a是偏移量，b是乘数，$H_{HFE}(u,v)$是高通滤波器的传递函数。

高频强调滤波优点是图像中由低频分量引起的灰色调得以保留。由较窄灰度级范围内的灰度表征的图像是直方图均衡化的理想选择。