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

我正在尝试构建2D盒函数的FFT和2D高斯函数的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");

rayryeng - R.. 5

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

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

第一个问题是你正在使用*,这是矩阵乘法.在频域中,逐元素乘法相当于频域中的卷积,因此需要使用.*.第二个问题是你也不需要.^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'标志确保输出大小是两个信号中的最大值,即高斯滤波器或Box滤波器.你会看到它是相同的曲线,为简洁起见,我不会在这里展示它.

但是,我们可以比较两个结果,看看它们是否是相同的元素.执行此操作的方法是确定减去结果中的每个元素是否小于某个阈值...比如.. 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进行测试.}

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

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

第一个问题是你正在使用*,这是矩阵乘法.在频域中,逐元素乘法相当于频域中的卷积,因此需要使用.*.第二个问题是你也不需要.^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'标志确保输出大小是两个信号中的最大值,即高斯滤波器或Box滤波器.你会看到它是相同的曲线,为简洁起见,我不会在这里展示它.

但是,我们可以比较两个结果,看看它们是否是相同的元素.执行此操作的方法是确定减去结果中的每个元素是否小于某个阈值...比如.. 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进行测试.}