第四章：频率域滤波

Joyrich.

已于 2024-07-03 08:12:32 修改

阅读量963

点赞数 26

文章标签：算法

于 2024-06-28 23:34:00 首次发布

本文链接：https://blog.youkuaiyun.com/qq_63029071/article/details/139904895

版权

4.3.2 空间域滤波和频率域滤波之间的对应关系

4.1一维傅里叶变换及其反变换

任何周期函数都可以表示为不同频率的正弦和或余弦和的形式，每个正弦和/或余弦和乘以不同的系数(现在称这个和为傅里叶级数)。

非周期的函数(曲线是有限的情况)也可以用正弦/或余弦乘以加权函数的积分来表示。在这种情况下的公式就是傅里叶变换。

单变量连续函数fx)的傅里叶变换F(u)定义为等式:

$F(\mu )=\int_{-\infty }^{+\infty }f(x)^{}e^{^{-j2\pi\mu x}}dx$

反变换定义为

$f(x)=\int_{-\infty }^{+\infty }F(\mu )^{}e^{^{j2\pi\mu x}}d\mu$

单变量离散函数f(r)(其中x=0,1,2...,m-1)的傅里叶变换由以下等式给出:

$F(\mu )=\frac{1}{M}\sum_{x=0}^{M-1}f(x)e^{-j2\pi \mu x/M}$ $for$ $u=0,1, ,M-1$

同样。给出F(u)，能用反DFT来获得

$f(x)=\sum_{u=0}^{M-1}F(u)e^{j2\pi ux/M}$

4.2 二维DFT及其反变换

二维连续傅里叶变换与反变换

$F(u,v)=\int_{-\infty }^{+\infty}\int_{-\infty }^{+\infty}f(x,y)e^{-j2\pi(ux+vy)}dxdy$

$f(x,y)=\int_{-\infty }^{+\infty}\int_{-\infty }^{+\infty}F(u,v)e^{j2\pi(ux+vy)}dudv$

二维离散傅里叶变换与反变换

$F(u,v)=\frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)e^{-j2\pi(ux/M+vy/N)}$

$f(x,y)=\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}F(u,v)e^{j2\pi(ux/M+vy/N)}$

对于图像尺寸为M*N的函数f(x, y)的二维离散傅里叶变换为：

$F(0,0)=\frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)$

给出F(u, v)，可通过反DFT得到f(x, y）如下：

$f(x,y)=\sum_{u=0}^{M-1}\sum_{v=0}^{N-1}F(u,v)e^{j2\pi(ux/M+vy/N)}$

二维离散的性质：

平移，用指数项乘以f(x,y)将使DFT的原点移到点 $(u_0,v_0)$ ;反之，用负指数项乘以F(x,y)将使f(x,y)的原点移到点 $(x_0,y_0)$
旋转, 若f(x,y)旋转 $\theta _0$ 角度，则F(u,v)也旋转相同的角度，反之若F(u,v)旋转一个角度,f（x,y)也旋转相同角度
周期性

4.3 频率域滤波

傅里叶变换的结果，与图像中的强度变化模式具有一定的联系。例如，变化最慢的频率成分(u=v=0)对应一幅图像的平均灰度级

$F(0,0)=\frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)$

反映的是对图像每个像素点的灰度级求和然后除以MN，得到平均灰度级。低频对应着图像的总体灰度级的显示，高频对应图像中变化快的分量(图像的细节)。

频率域滤波步骤：

用 $(-1)^{x+y}$ 乘以输入图像来进行中心变换
由(1)计算图像的DFT，即F(u, v)
用滤波器函数H(u,v)乘以F(u,v)
计算(⑶中结果的反DFT
得到(4)中结果的实部用 $(-1)^{x+y}$
乘以(5)中的结果
H(u,v)被称为滤波器(滤波器传递函数)。

4.3.1 陷波滤波器

图像的平均值由F(0,0)给出。如果在频率域中设置此项为零,并进行度变换,那么结果图像的平均值将为零。滤波函数可选为

$H(u,v)=\begin{cases} 0 & \text (u,v)= (M/2,N/2)\\ 1 & \text other \end{cases}$

此滤波器可以设置F(0,0)为零(正如所希望的那样)，而保留其他傅里叶变换的频率成分不变。处理后的图像(具有0平均值)可以通过对H(u,v)F(u,v)进行傅里叶反变换来获得

4.3.2 空间域滤波和频率域滤波之间的对应关系

空间域和频率域之间最基本的联系是由卷积定理的有关结论建立的。

在空间域中将滤波的模板在图像中逐像素移动，并对每个像素进行指定数量的计算的过程就是卷积过程。形式上,大小为M×N的两个函数f(x, y)和h(x, y)的离散卷积表示为,f (x, y) * h(x, y)，并定义如下:

$f(x,y)\times h(x,y)=\frac{1}{MN}\sum_{m=0}^{M-1}\sum_{n=0}^{N-1}f(m,n)h(x-m,y-n)$

除了前面的常数、负号以及求和的上下限之外，整个表达式与空间域的线性滤波相似。特别是负号只说明函数h关于原点镜像对称。这是卷积定义中自带的。

上式是一种实现：

关于原点翻转函数
通过改变(x,y)的值相对于一个函数移动到另外一个函数
对每一个(x,y)的位移值，计算所有m和n值成绩的和。（(x,y)的位移值是以整数增加的，当函数不再重叠时停止）

用F(u,v)和H(u,v)分别表示f(x,y)和h(x,y)的傅里叶变换，卷积定理说明f(x,y)*h(x,y)和F(u,v)H(u,v)组成傅里叶变换对。类似的结果是频率域的卷积简化为空间域的乘法，形式上表示如下：

$f(x,y) *h(x,y) \Leftrightarrow F(u,v)H(u,v);f(x,y)h(x,y)\Leftrightarrow F(u,v)*H(u,v)$

注意，前述的所有函数均为相同尺寸M*N。因此，在实际中，指定一个频率域的滤波器，然后进行反变换以计算相同尺寸的空间域的相应滤波器，但这种方法从计算的角度来看并不能解决太大问题。如果两个滤波器是相同尺寸，那么通常在频率域进行滤波计算更为有效。但是，在空间域更适用于更小的滤波器。这正是我们所感兴趣的联系。滤波器在频率域中更为直观，但在空间域使用更小的滤波器模板更为明智。

由此，可以在频率域指定滤波器，做反变换，然后在空间域使用结果滤波器作为在空间域构建更小的空间滤波模板的指导。

4.4 平滑的频率域滤波器

边缘和其他尖锐变化(如噪声)在图像的灰度级中主要处于傅里叶变换的高频部分。因此,平滑(模糊)可以通过衰减指定图像傅里叶变换中高频成分的范围来实现。

目标是选择一个滤波器变换函数H(u, v)以通过衰减F(u, v)的高频成分产生G,(u, v)。

主要使用低通滤波器，（只通过低频部分，过滤掉高频部分）

理想低通滤波器
巴特沃思低通滤波器
高斯低通滤波器

4.4.1 理想低通滤波器

最简单的低通滤波器是"截断"傅里叶变换中所有高频成分，这些成分处在距变换原点的距离比指定距离 $D_0$ 要远得多的位置。这种滤波器称为二维理想低通滤渡器,其变换函数为:

$H(u,v)=\begin{cases} 1 & \text D(u,v)\leqslant D_{0}\\ 0 & \text D(u,v)> D_{0} \end{cases}$

其中， $D_0$ 是指定的非负数值，D(u, v)是(u, v)点距频率矩形原点的距离。由于变换被中心化了，如果要研究的图像尺寸为M×N,从点(u, v)到傅里叶变换中心(原点)的距离如下所示:
$D(u,v)=[(u-M/2)^{2}+(v-N/2)^{2}] ^{\frac{1}{2}}$

理想低通滤波器模糊和振铃特性可参考卷积定理来解释。

卷积定理指出在空间域的相应过程:g( x ,y) = h( z, y)* f(x , y)

其中h(x, y)是滤波器变换函数H(u,v)的反变换。

4.4.2 巴特沃思低通滤波器

n阶布特沃斯低通滤波器(BLPF)的传递函数（且截止频率距原点的距离为 $D_0$ ）的定义如下：

$H(u,v)=\frac{1}{1+[ D(u,v)/D_0] ^{2n}}$

其透视图、图像显示以及BLPF的径向横截面如图

4.4.3 高斯低通滤波器

$H(u,v)=e^{-D^2(u,v)/2D_0^2}$

其中， $D_0$ 是截止频率。当当D(u,v)= $D_0$ 时，滤波器下降到它最大值的0.607处。高斯滤波器有一个重要的特性，那就是在高斯滤波器中是没有振铃的。其平滑效果常不如布特沃斯低通滤波器。

4.4 频率域锐化滤波器

灰度级的边缘和其他地方的急剧变化与高频成分有关，图像的锐化能够在频率域用高通滤波处理实现，而衰减低频成分并不会扰乱傅里叶变换的高频信息。

高通滤波器，可以通过对低通滤波器进行精确的反操作来得到。

4.5 频率域的拉普拉斯算子

频率域的拉普拉斯算子可由如下滤波器实现：

$H(u,v)=-4\pi^2(u^2+v^2)$

之前的公式在处理数字图像时，还需进行比例的变换。例如对于大小为M*N的数字图像f(x, y)，实际的频率域拉普拉斯算子为：

$H(u,v)=-4\pi^2[(u/M)^2+(v/M)^2]$

如果不做这个处理，u和v的取值范围从0至M-1和0至N-1的，u^2 + v^2其实是一个非常大的值，这样不适合在数字图像中做处理的。

计算F(u,v)时需要先中心化，也就是将中心移到M/2和N/2处，则此时为：

$H(u,v)-4\pi^2\left \{ [(u-M/2)/M]^2+[(v-N/2)/N]^2 \right \}$

4.6 同态滤波器

照度-反射模型可用来开发一种频域处理,通过同时进行灰度范围的压缩和对比度增强来改进一幅图像的外观。
$f(x,y)=i(x,y)r(x,y)$

采用上述概念的增强方法总结于下图。这个方法是基于称为同态系统的一类系统的特例。在这个特殊应用中，方法的关键是将照射分量和反射分量分开。同态滤波函数H(u,v)能分别对两个分量进行操作.

二、实验

1.理想低通滤波

I = imread('01.jpg');  % 读入图像
I = rgb2gray(I);
I = im2double(I);
M = 2 * size(I, 1);   % 滤波器的行数
N = 2 * size(I, 2);   % 滤波器的列数
u = -M/2:(M/2-1);
v = -N/2:(N/2-1);
[U, V] = meshgrid(u, v);
D = sqrt(U.^2 + V.^2);
D0 = 100;               % 截止频率
H = double(D <= D0);   % 理想低通滤波器
J = fftshift(fft2(I, size(H, 1), size(H, 2)));  % 空域图像转换到频域
K = J .* H;                                     % 滤波处理
L = ifft2(ifftshift(K));                        % 傅里叶反变换 
L = L(1:size(I, 1), 1:size(I, 2));
% 显示原始图像和滤波后的图像
figure;
% 显示原始图像
subplot(221);
imshow(I);
title('原始图像');
% 显示滤波后的图像
subplot(222);
imshow(L);
title('滤波后的图像');
% 显示原始图像的频域
subplot(223);
imagesc(log(1 + abs(J)));
title('原始图像的频域');
% 显示滤波后图像的频域
subplot(224);
imagesc(log(1 + abs(K)));
title('滤波后图像的频域');

2.理想高通滤波

I = imread('01.jpg');  % 读入图像
I = rgb2gray(I);
I = im2double(I);
M = 2 * size(I, 1);   % 滤波器的行数
N = 2 * size(I, 2);   % 滤波器的列数
u = -M/2:(M/2-1);
v = -N/2:(N/2-1);
[U, V] = meshgrid(u, v);
D = sqrt(U.^2 + V.^2);
D0 = 100;               % 截止频率
H = double(D >= D0);   % 理想高通滤波器
J = fftshift(fft2(I, size(H, 1), size(H, 2)));  % 空域图像转换到频域
K = J .* H;                                     % 滤波处理
L = ifft2(ifftshift(K));                        % 傅里叶反变换 
L = L(1:size(I, 1), 1:size(I, 2));
% 显示原始图像和滤波后的图像
figure;
% 显示原始图像
subplot(221);
imshow(I);
title('原始图像');
% 显示滤波后的图像
subplot(222);
imshow(L);
title('滤波后的图像');
% 显示原始图像的频域
subplot(223);
imagesc(log(1 + abs(J)));
title('原始图像的频域');
% 显示滤波后图像的频域
subplot(224);
imagesc(log(1 + abs(K)));
title('滤波后图像的频域');