通过 FFT 执行 2D 卷积时的奇怪行为

Question

我正在尝试构建二维盒函数的 FFT 和二维高斯函数的 FFT 的乘积。之后，我找到了逆 FFT，因此我期待这两个函数的卷积。但是，如下所示，我得到了一个奇怪的片面结果。结果显示在subplot 的右下角。

我为重现上述subplot 而编写的 Octave 代码以及我为构建卷积而执行的计算如下所示。谁能告诉我我做错了什么？

clear all;
clc;
close all;
% domain on each side is 0-9
L = 10;
% num subdivisions
N = 32;
delta=L/N;
sigma = 0.5;
% get the domain ready
[x,y] = meshgrid((0:N-1)*delta);
% since domain ranges from 0-(N-1) on both sdes
% we need to take the average
xAvg = sum(x(1, :))/length(x(1,:));
yAvg =  sum(y(:, 1))/length(x(:,1));
% gaussian
gssn =  exp(- ((x - xAvg) .^ 2 + (y - yAvg) .^ 2) ./ (2*sigma^2));


function ret = boxImpulse(a,b)
    n = 32;
    L=10;
    delta = L/n;
    nL = ((n-1)/2-3)*delta;
    nU = ((n-1)/2+3)*delta;
    if ((a >= nL) && (a <= nU) && ( b >= nL) && (b <= nU) )
        ret=1;
    else
        ret=0;
    end
    ret;
endfunction

boxResponse = arrayfun(@boxImpulse, x, y);

subplot(2,2,1);mesh(x,y,gssn); title("gaussian fun");
subplot(2,2,2);mesh(x,y, abs(fft2(gssn)) .^2); title("fft of gaussian");
subplot(2,2,3);mesh(x,y,boxResponse); title("box fun");
inv_of_product_of_ffts =  abs(ifft2(fft2(boxResponse) * fft2(gssn))) .^2 ;
subplot(2,2,4);mesh(x,y,inv_of_product_of_ffts); title("inv of product of fft");

Answer 1

让我们首先解决您最明显的错误。这是您在频域中计算卷积的方式：

inv_of_product_of_ffts =  abs(ifft2(fft2(boxResponse) * fft2(gssn))) .^2 ;

第一个问题是您使用的是*，即矩阵乘法。在频域中，element-wise乘法相当于频域中的卷积，所以需要使用.*来代替。第二个问题是您也不需要.^2 术语和abs 操作。您正在执行卷积，然后找到绝对值并对每个项进行平方。这不是必需的。删除abs 和.^2 操作。

因此，你需要的是：

inv_of_product_of_ffts = ifft2(fft2(boxResponse).*fft2(gssn)));

---

但是，您将得到的是这个结果。让我们把它放在一个新的图中，而不是 subplot^*:

figure;
mesh(x,y,ifft2(fft2(boxResponse).*fft2(gssn))); 
title('Convolution... not right though');

您可以看到它是正确的结果...但它没有居中...这是为什么呢？ 这实际上是在频域计算卷积时最常见的问题之一。事实上，即使是最有经验的人也会遇到这个问题，因为他们不了解 FFT 的工作原理。

这是 MATLAB 和 Octave 计算 2D FFT 的结果。具体来说，MATLAB 和 Octave 定义了一个 meshgrid 的坐标，从 0,1,...M-1 的水平和 0,1,...N-1 的垂直，给我们一个 M x N 结果。这意味着原点 / DC 分量位于矩阵的 左上角 角，而不是我们通常定义的 中心。具体来说，传统的 2D FFT 定义了从-(M-1)/2, ..., (M-1)/2 的水平坐标和-(N-1)/2, ..., (N-1)/2 的垂直坐标。

参考上述内容，您使用中心定义信号，假设它是原点而不是左上角。为了弥补这一点，您需要添加一个fftshift（MATLAB doc，Octave doc），以便 FFT 的输出现在以原点而不是左上角为中心，以使事情恢复原状是，因此你真的需要：

figure;
mesh(x,y,fftshift(ifft2(fft2(boxResponse).*fft2(gssn)))); 
title('Convolution... now it is right');

如果您想仔细检查我们的结果是否正确，您可以在空间域执行直接卷积，我们可以在频域比较方法之间的结果。

这是您在空间域中直接计算结果的方式：

figure;
mesh(x,y,conv2(gssn,boxResponse,'same'));
title('Convolution... spatial domain');

conv2 (MATLAB doc, Octave doc) 在两个信号之间执行 2D 卷积，'same' 标志确保输出大小是两个信号中最大的，即高斯滤波器或箱形滤波器。你会看到它是同一条曲线，为了简洁，我不会在这里展示它。

但是，我们可以比较两个结果，看看它们在元素方面是否相同。一种方法是确定减去结果中的每个元素是否小于某个阈值......比如说......1e-10：

>> out1 = conv2(boxResponse, gssn, 'same');
>> out2 = fftshift(ifft2(fft2(boxResponse).*fft2(gssn)));
>> all(abs(out1(:)-out2(:)) < 1e-10)

ans =

     1

这意味着这两者确实是相同的。

---

希望这是有道理的！请记住，由于您将信号定义在原点位于中心的位置，因此一旦找到逆 FFT，您必须将事物移向后移动，以便原点现在位于中心，而不是顶部-左角。

---

^{*：次要注意 - 此答案中生成的所有图都是使用 MATLAB R2015a 生成的。但是，该代码在 Octave 中完全可以工作 - 使用 Octave 4.0.0 进行测试。}