1、基础

卷积定理

折叠误差补零

当处理DFT时，图像及其变换是周期的。在周期接近函数非零部分的持续周期时，对周期函数进行卷积会引起相邻周期的串扰。

这种称为 折叠误差 的串扰可通过补零方法来避免。

假设函数 f(x,y) 和 h(x,y) 的大小分别为 A x B 和 C x D。通过对 f 和 g 补零，构造两个大小均为 P x Q 的扩展函数。

按如下方式进行选择可以避免折叠误差：

下面 paddedsize 的函数可以满足上面等式的 P 和 Q 的最小偶数值。

function PQ = paddedsize(AB, CD, PARAM)

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);

P = 2^nextpow2(2*m);

PQ = [P, P];

elseif (nargin == 3) && strcmp(PARAM, 'pwr2')

m = max([AB CD]);

P = 2^nextpow2(2*m);

PQ = [P, P];

else

error('Wrong number of inputs.')

end

示例

[f, revertcalss] = tofloat(f);

sig = 10;

PQ = paddedsize(size(f));

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

Hp = lpfilter('gaussian', PQ(1), PQ(2), 2*sig); %滤波器大小是未填充时的两倍

Gp = Hp .* Fp;

gp = ifft2(Gp);

gpc = gp(1:size(f, 1), 1:size(f, 2));

gpc = revertcalss(gpc);

imshow(gp)

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);

用于频率域滤波的M函数

function g = dftfilt(f, H, classout)

[f, revertcalss] = tofloat(f);

F = fft2(f, size(H, 1), size(H, 2))

g = ifft2(H .* F);

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

if nargin == 2 || strcmp(classout, 'original')

g = revertcalss(g);

elseif strcmp(classout, 'fltpoint')

return

else

error('Undefined class fot the output image.')

end

3、从空间滤波器获得频域滤波器

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

当滤波器大约有 32 个或更多元素时，使用 FFT 算法的滤波处理要比空间实现快。

函数 freqz2 可以将空间滤波器转换为频率域滤波器

H = freqz2(h, R, C)

示例

h = fspecial('sobel')'

freqz2(h)

% 频率域滤波器

PQ = paddedsize(size(f));

H = freqz2(h, PQ(1), PQ(2));

H1 = ifftshift(H)

函数 dftfilt 用于处理频率域滤波

gf = dftfilt(f, H1);

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

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

function [U, V] = dftuv(M, N)

u = single(0:(M - 1));

v = single(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);

% ----------------------- %

[U, V] = dftuv(8, 5);

DSQ = hypot(U, V)

低通(平滑)频率域滤波器

理想低通滤波器(ILPF)

巴特沃斯低通滤波器(BLPF)

高斯低通滤波器

示例，对图像 f 应用一个高斯低通滤波器，其中 D0 为所填充图像宽度的 5%。

[f, revertclass] = tofloat(f);

PQ = paddedsize(size(f));

[U, V] = dftuv(PQ(1), PQ(2));

D = hypot(U, V);

D0 = 0.05 * PQ(2);

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

H = exp(-(D .^ 2) / (2 * (D0^2)));

g = dftfilt(f, H);

g = revertclass(g)

低通滤波器整合函数

function H = lpfilter(type, M, N, D0, n)

[U, V] = dftuv(M, N);

D = hypot(U, V);

switch type

case 'ideal'

H = single(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

5、绘制线框图和表面图

绘制二维函数的最简单方法是使用函数 mesh

mesh(H(1:k:end, 1:k:end))

默认绘出彩色网格图命令为

colormap([0 0 0])

关闭网格和坐标轴命令

grid off

axis off

观察点由函数 view 控制

view(az, el) %az,el 代表方位角和仰角

view(3) %默认观察点

示例，绘制一个高斯低通滤波器

H = fftshift(test2('gaussian', 500, 500, 50));

mesh(double(H(1:10:500, 1:10:500)))

axis tight

colormap([0 0 0])

untitled.jpg

表面图

surf(H)

6、高通(锐化)频率域滤波器

理想高通滤波器(IHPF)

巴特沃斯高通滤波器(BHPF)

高斯高通滤波器

基于 lpfilter 构建一个高通滤波器函数

function H = hpfilter(type, M, N, D0, n)

if nargin == 4

n = 1;

end

H1p = lpfilter(type, M, N, D0, n);

H = 1 - H1p;

% ---------------------------%

H = fftshift(hpfilter('ideal', 500, 500, 50));

高频强调滤波

高通滤波器偏离了直流项，将图像的平均值降低为0.补偿方法之一是给高通滤波器加上一个偏移量。把偏移量与将滤波器乘以一个大于1的常数结合起来的方法就称为 高频强调滤波。

示例，联合使用高频强调滤波和直方图均衡

PQ = paddedsize(size(f));

D0 = 0.05 * PQ(1);

HBW = hpfilter('btw', PQ(1), PQ(2), D0, 2);

H = 0.5 + 2*HBW;

gbw = dftfilt(f, HBW, 'fltpoint');

gbw = gscale(gbw);

ghf = dftfilt(f, H, 'fltpoint');

ghf = gscale(ghf);

ghe = histeq(ghf, 256);