我谈冈萨雷斯对频域滤波的误解——快速卷积与频域滤波之间的关系

在Rafael Gonzalez和Richard Woods所著的《数字图像处理》中，Gonzalez对频域滤波是有误解的，在频域设计滤波器不是非得图像和滤波器的尺寸相同，不是非得在频域通过乘积实现。相反，FIR滤波器设计都是构造空域脉冲响应。一般的原则是，小尺寸的滤波器在空域通过卷积实现更快，大尺寸的滤波器在频域通过频域滤波实现更快。

这里就引出一个快速卷积的概念，实际上原理上很简单，基础就是卷积定理。只是很多人没想过这个问题。

直接计算两个长度分别为N和M的离列序列的线性卷积需要 $O (N * M)$ 次乘法运算，当N和M较大时，这将是一个相当耗时的过程。

快速卷积是指使用快速傅里叶变换（Fast Fourier Transform, FFT）来加速两个序列的线性卷积的过程。

利用FFT可以将时间域中的卷积转换为频率域中的点乘运算，从而减少所需的计算量。具体步骤如下：

补零：首先，为了防止循环卷积带来的边界效应，通常需要将两个序列都扩展到至少N+M-1的长度，通常是选择2的幂次大小，这样更有利于FFT算法的实现。
执行FFT：对补零后的两个序列分别进行快速傅里叶变换。
频域相乘：将两个变换后的序列逐点相乘。
逆FFT：对相乘结果执行逆快速傅里叶变换（Inverse Fast Fourier Transform, IFFT），得到最终的卷积结果。
调整结果：由于使用了IFFT，可能需要对输出结果做一些调整，比如去除因补零而引入的额外部分。

这种方法的时间复杂度主要由FFT/IFFT决定，对于长度为 $L$ 的序列，其时间复杂度大约是 $\log L)$ 。

Using the typical convolution formula to compute the one-dimensional convolution of a P-element sequence A with Q-element sequence B has a computational complexity of $O (PQ)$ . However, the discrete Fourier transform (DFT)
can be used to implement convolution as follows:

Compute the L-point DFT of A, where $\geq P + Q - 1$ .
Compute the L-point DFT of B.
Multiply the two DFTs.
Compute the inverse DFT to get the convolution.

The overall computational complexity of these steps is $O(L\log L)$ . For P and
Q sufficiently large, then, using the DFT to implement convolution is a
computational win.

See the function fftfilt in the Signal Processing Toolbox.

在图像处理中，滤波通常用卷积实现。使用傅里叶变换对两个矩阵执行快速卷积。傅里叶变换的一个关键特性是两个傅里叶变换相乘对应于相关联的空间函数的卷积。此特性与快速傅里叶变换一起构成了快速卷积算法的基础。

基于 FFT 的卷积方法最常用于大型输入。对于小型输入，使用 imfilter 函数通常更快。

图中不同尺寸模糊核，两种方法的时间开销。卷积计算随尺寸持续增加，而基于FFT的运算基本稳定。
在这里插入图片描述

clear
close all
%% Read Image
name = 'Tsinghua Gate 512';
suffix = '.bmp';
Img = im2double(imread( ['images/', name, suffix] )) ;
sigmas = 21:31;
for k = 1:numel(sigmas)sigma = sigmas(k);hsize(k) = 2*ceil(3*sigma) + 1;psf   = fspecial( 'gaussian' , hsize(k), sigma );t_fft(k) = timeit(@() imfilter_fft(Img, psf)); %full convolutiont_normal(k) = timeit(@() imfilter(Img,psf,'circular'));
endfigure;
plot(hsize, t_normal, 'x-', hsize, t_fft, '*-', 'LineWidth',1)
grid on
xlabel('Kernel Size')
ylabel('Time/s')
box off
legend({'imfilter computation', 'FFT-based computation'},'Location','northwest')
title('Gaussian Filtering')%%
function y = imfilter_fft(x, psf)X = fft2(x);
OTF = psf2otf(psf, size(x));
y = ifft2(X.*OTF);end