频域滤波中默认的边界条件——补零与不补零（答作者问）

这个问题源自于Rafael Gonzalez的《数字图像处理》中的这幅图，为什么他频域滤波要将图像零延拓到二倍尺寸？

完全没有没要，既浪费计算，又浪费空间。

注：图(d)错误，见这里。
这个他给的解释，并没有从理论高度做解释。

三个基础知识

卷积边界延拓的方式

离散傅里叶变换的周期性

线性卷积和循环卷积等效的条件

离散傅里叶变换的卷积定理

循环卷积对应于离散傅里叶变换的乘积，即离散傅里叶变换的卷积定理。

如果 X ( k ) = D F T { x ( n ) } X(k) = DFT\{x(n)\} X(k)=DFT{x(n)}和 Y ( k ) = D F T { y ( n ) } Y(k) = DFT\{y(n)\} Y(k)=DFT{y(n)}，那么

X ( k ) Y ( k ) = D F T { { x ( n ) } ⊛ { y ( n ) } } X(k)Y(k) = DFT\{\{x(n)\} \circledast \{y(n)\}\} X(k)Y(k)=DFT{{x(n)}⊛{y(n)}}

这里，$⊛$表示由下式定义的循环卷积：

{ x ( n ) } ⊛ { y ( n ) } = ∑ m = 0 N − 1 x ( m ) y ( ( n − m ) m o d N ) . \{x(n)\} \circledast \{y(n)\} = \sum_{m=0}^{N-1} x(m)y((n-m) \mod N). {x(n)}⊛{y(n)}=m=0∑N−1​x(m)y((n−m)modN).

$y((n-m)\mathrm{mod}N)$表示循环移位。

Conventional Shift

Circular Shift

线性卷积是这样的，$h(-n)$移位，计算与$f(n)$对位相乘再求和，一直下去就是线性卷积的结果$g(n)$。

然而，不幸的是，离散傅里叶变换的周期性造成时域是循环卷积。某人优点很突出，缺点也很突出，我很欣赏他的优点，怎么办呢？那只能解决掉他的缺点。

零延拓，补零的边界条件

循环卷积没办法了，那就想办法让循环卷积与线性卷积等价，于是就有了零延拓。对信号补零到线性卷积结果的长度。图中虽然循环移位，但是由于信号$f(n)$补零了，对应的那两个点就是零，相乘也是零。换句话说就是把位置给留出来。这样仍然是线性卷积的结果。

由于是线性卷积，那默认的边界延拓方式就是零延拓。看图说话，相当于$f(n)$补零后周期延拓，$h(-n)$移位，那就是零延拓的线性卷积。

看图像的结果，这里看不出大小。

因为fft2本身可以做任意点的傅里叶变换，所以无需提前延拓，也就没显示。

空域滤波器，频域实现。

blurtype = 'gaussian';
sigma = 2;
hsize = 2 * ceil(3 * sigma) + 1;
psf_gauss = fspecial('gaussian', hsize, sigma);
offset = floor(size(psf_gauss)/2);
S = size(Img) + size(psf_gauss) - 1;
F = fft2(Img, S(1), S(2));
H = fft2(psf_gauss, S(1), S(2));G=F.* H;
g=im2uint8(real(ifft2(G)));

那Gonzalez为什么要补零为二倍呢？这是由于他在频域直接对模拟滤波器采样为$M\times N$点，对应的空域就是$M\times N$点（频域滤波器不这样设计，很离谱，暂且不论）。所以为了不产生混叠，延拓的图像尺寸至少是$(2M-1)\times (2N-1)$。他的大错不少，小错不断，他认为是$2M\times 2N$也就见怪不怪了。

频域滤波器

filename = 'cameraman';
suffix = '.tif';
Img = im2double(im2gray(imread([filename, suffix])));
[M,N]=size(Img);
D0 = 0.15;
[U,V] = freqspace([2*M 2*N],'meshgrid');
d = sqrt(U.^2+V.^2);
H = exp(-(d/(2*D0.^2)));F = fftshift(fft2(Img,2*M,2*N));
G = F.* H;
g=im2uint8(real(ifft2(ifftshift(G))));

不补零，周期延拓的边界条件

回头再看不零延拓的话，会怎么样？头部数据混叠，尾部数据丢失。

但是图像处理的滤波通常是考虑same size的结果，可以先做循环移位。

与线性卷积的结果比较，中间的五个值是相同的，这是因为边界受周期延拓影响了。

>> gcgc =1.4641    1.5738    0.9545    0.5790    0.3512    0.2130    0.8149
>> gg =0.7071    1.4289    1.5738    0.9545    0.5790    0.3512    0.2130    0.1078    0.0352

所以MATLAB提出了psf2otf函数，先做循环移位，再计算离散傅里叶变换。如果有空域的卷积核，通过这个函数实现频域滤波。相当于周期延拓的边界条件。

空域滤波器，频域实现。

[M,N]=size(Img);
F = fft2(Img,M,N);blurtype = 'gaussian';
sigma = 2;
hsize = 2 * ceil(3 * sigma) + 1;
offset = (hsize-1)/2;
psf_gauss = fspecial('gaussian', hsize, sigma);
H = psf2otf(psf_gauss, [M,N]); %这里使用fft2就会产生错误，很多刊物这里有误。
G=F.*H;
g=im2uint8(real(ifft2(G)));

注意psf2otf的地方。

频域滤波器

[M,N]=size(Img);
D0 = 0.3;
[U,V] = freqspace([M N],'meshgrid');
d = sqrt(U.^2+V.^2);
H = exp(-(d/(2*D0.^2)));F = fftshift(fft2(Img));
G = F.* H;g=im2uint8(real(ifft2(ifftshift(G))));

结论

完全没有必要延拓，可以对原图像尺寸的频率直接频域滤波或者频域实现，就是相当于周期延拓。本来空域也有这样的延拓方式，而且我觉得比零延拓的效果好。