matlab - 通过 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 结果。这意味着原点/直流分量位于矩阵的左上角，而不是我们通常定义的中心。具体来说，传统的 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 ) 在两个信号之间执行二维卷积，'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 进行了测试。}