matlab编程实现卷积定理,利用时域（空域）卷积定理进行图像滤波（Matlab 实现）（二）...

引言

在 利用时域(空域)卷积定理进行图像滤波(Matlab 实现)(一)中我们知道了傅里叶变换是满足时域(空域)卷积定理的，在 从离散傅里叶变换到离散余弦变换 —— 公式证明 中我们也看到了离散余弦变换实际上就是对经过对称平移后的图像进行的傅里叶变换，那么离散余弦变换是否也满足时域(空域)卷积定理呢，请看下文。

发现新大陆

我们以二维空间的图像为例，列出其傅里叶变换的时域(空域)卷积定理：

F [ f ( x , y ) ∗ g ( x , y ) ] = F [ f ( x , y ) ] ⋅ F [ g ( x , y ) ] \mathcal{F}[f(x,y)*g(x,y)]=\mathcal{F}[f(x,y)]\cdot\mathcal{F}[g(x,y)]F[f(x,y)∗g(x,y)]=F[f(x,y)]⋅F[g(x,y)]

那么对应到离散余弦变换公式应该是什么样的呢，是下面这样吗？

D [ f ( x , y ) ∗ g ( x , y ) ] = D [ f ( x , y ) ] ⋅ D [ g ( x , y ) ] \mathcal{D}[f(x,y)*g(x,y)]=\mathcal{D}[f(x,y)]\cdot\mathcal{D}[g(x,y)]D[f(x,y)∗g(x,y)]=D[f(x,y)]⋅D[g(x,y)]

Sorry，正确的写法是这样的：

D [ f ( x , y ) ∗ g ( x , y ) ] = D [ f ( x , y ) ] ⋅ F [ g ( x , y ) ] \mathcal{D}[f(x,y)*g(x,y)]=\mathcal{D}[f(x,y)]\cdot\mathfrak{F}[g(x,y)]D[f(x,y)∗g(x,y)]=D[f(x,y)]⋅F[g(x,y)]

同时， g ( x , y ) g(x,y)g(x,y) 要满足上下对称，左右对称。记住，必须满足对称条件，后面我会进行说明。

这里我使用 F \mathfrak{F}F 符号代表了一种特殊的傅里叶变换，它代表先将 g ( x , y ) g(x,y)g(x,y) 扩展到 2 * 2 倍的大小，其余位置都填充0，对扩展后的矩阵求傅里叶变换，之后再截取到原先大小。这里表达可能不是很清楚，见下图(假定 f ( x , y ) f(x,y)f(x,y) 为图像，g ( x , y ) g(x,y)g(x,y) 为 Laplacian 卷积核)：

图中的卷积核 g ( x , y ) g(x,y)g(x,y) 为了与图像 f ( x , y ) f(x,y)f(x,y) 大小匹配，事先进行了扩展，具体扩展方式见 利用时域(空域)卷积定理进行图像滤波(Matlab 实现)(一)。

我们还需注意，当原始图像矩阵为上下且左右对称时，傅里叶变换后的结果也是上下且左右对称的。上图中第 2 步到第 3 步，尽管扩展后的矩阵并非是对称的，但是可以证明其傅里叶变换依然是上下且左右对称的(先将 3 * 3 的卷积核移到矩阵中心，再通过时移定理证明对称性)，不过只有当卷积核也是上下且左右对称时，结论才成立。

我们再仔细看看这个过程，有没有发现，整个操作跟离散余弦变换过程非常像，只是我们在扩展时，没有对图像做对称，仅仅是填充了 0。离散余弦变换因为开始对图像做了上下和左右对称，所以将对称图像做傅里叶变换后的结果也是上下且左右对称的，因此我们才能把离散余弦变换 D ( u , v ) D(u,v)D(u,v) 的 u , v u,vu,v 范围限制在 [ 0 , N − 1 ] [0, N-1][0,N−1] 而不损失信息。

离散余弦变换过程如下图所示：

综上，最终要证明的公式

D [ f ( x , y ) ∗ g ( x , y ) ] = D [ f ( x , y ) ] ⋅ F [ g ( x , y ) ] \mathcal{D}[f(x,y)*g(x,y)]=\mathcal{D}[f(x,y)]\cdot\mathfrak{F}[g(x,y)]D[f(x,y)∗g(x,y)]=D[f(x,y)]⋅F[g(x,y)]

可以简单的用下图表达(傅里叶变换后都要进行裁剪才能完全等价)：

图 A 为 f ( x , y ) ∗ g ( x , y ) f(x,y)*g(x,y)f(x,y)∗g(x,y) 做上下、左右对称后的图。图 B、C 为上两张变换示意图中 f ( x , y ) ,   g ( x , y ) f(x,y),\ g(x,y)f(x,y),g(x,y) 经过第一步变换后的图。

仔细对照一下就可以发现，这个等价的表达式与傅里叶变换的时域卷积定理能够完全匹配。当然，前提是卷积核必须是上下、左右对称的，否则不能将图 A 还原为图 B 与图 C 的卷积。

Matlab实现

代码：

clear, clc, close all

img_raw = imread('cameraman.tif');

img = double(img_raw);

freq_img = dct2(img); % 原图像的离散余弦变换

gaussian = [1, 2, 1; 2, 4, 2; 1, 2, 1] / 16; % Gaussian 模板

laplacian = [0, 1, 0; 1, -4, 1; 0, 1, 0]; % Laplacian 模板

filter = zeros(size(img) * 2); % 滤波器大小先扩展为图像 2 * 2 大小

filter_gau = filter; filter_gau([end,1,2], [end,1,2]) = gaussian;

filter_lap = filter; filter_lap([end,1,2], [end,1,2]) = laplacian;

freq_gau = fft2(filter_gau); % Gaussian 滤波器的傅里叶变换

freq_lap = fft2(filter_lap); % Laplacian 滤波器的傅里叶变换

freq_gau = freq_gau(1:256, 1:256); % 滤波器裁剪

freq_lap = freq_lap(1:256, 1:256); % 滤波器裁剪

% 利用新证明的定理得到滤波后的图像

img_gau = idct2(freq_gau .* freq_img); % Gaussian 滤波图像

img_gau = uint8(real(img_gau));

img_lap = idct2(freq_lap .* freq_img); % Laplacian 滤波图像

img_lap = uint8(real(img_lap));

% 使用自带卷积函数

img_gau_conv = uint8(conv2(img, gaussian, 'same'));

img_lap_conv = uint8(conv2(img, laplacian, 'same'));

% 画图

subplot(2, 3, 1); imshow(img_raw); title('原图');

subplot(2, 3, 2); imshow(img_gau); title('Gaussian 滤波图(新定理)');

subplot(2, 3, 3); imshow(img_gau_conv); title('Gaussian 滤波图(conv2)');

subplot(2, 3, 5); imshow(img_lap); title('Laplacian 滤波图(新定理)');

subplot(2, 3, 6); imshow(img_lap_conv); title('Laplacian 滤波图(conv2)');

再看图像矩阵，以 Laplacian 滤波图为例：

定理计算图像矩阵：

conv2 函数计算图像矩阵：

除外围像素，其他像素灰度完全一致，因此定理成立。