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

Question

我正在尝试构造二维框函数的 FFT 和二维高斯函数的 FFT 的乘积。之后，我找到了逆 FFT，结果我期待这两个函数的卷积。但是，我得到了一个奇怪的单方面结果，如下所示。结果显示在 subplot.

的右下角

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

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 结果。这意味着原点/直流分量位于矩阵的 top-left 角，而不是我们通常定义的 centre 。具体来说，传统的 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，就必须 shift 使原点现在位于中心，而不是左上角。

^{*：小注 - 此答案中生成的所有图都是使用 MATLAB R2015a 生成的。但是，该代码在 Octave 中完全有效 - 使用 Octave 4.0.0 进行了测试。}

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

Weird behavior when performing 2D convolution by the FFT

matlab

signal-processing

fft

gaussian

octave