Matlab - 高斯的 FFT - 等价

Question

简单的问题：

我在 Matlab 中绘制了一个具有一定分辨率的二维高斯函数。我用方差或 sigma = 1.0 进行测试。我想将它与 FFT(Gaussian) 的结果进行比较，这应该会产生另一个方差为 (1./sigma) 的高斯。由于我使用 sigma = 1.0 进行测试，我认为我应该获得两个等效的 2D 内核。

即

g1FFT = buildKernel(rows, cols, mu, sigma)    % uses normpdf over arbitrary resolution (rows, cols, 3) with the peak in the center

构建内核：

function result = buildKernel(rows, cols, mu, sigma)  

result = zeros(rows, cols, 3);

center_w = floor(cols / 2);
center_h = floor(rows / 2);

for i = 1:rows
    for j = 1:cols
        distance = sqrt((center_w - j).^2 + (center_h - i).^2);
        g_val = normpdf(distance, mu, sigma);
        result(i, j, :) = g_val;
    end
end

% normalize so that kernel sums to 1
sumKernel = sum(result(:));
result = result ./ sumKernel;    

end

我正在使用 mu = 0.0（始终）和方差或 sigma = 1.0 进行测试。我想将它与 FFT(Gaussian) 的结果进行比较，这应该会产生另一个方差为 (1./sigma) 的高斯。

即

g1FFT = circshift(g1FFT, [rows/2, cols/2, 0]);       % fft2 expects center to be in corners 
freq_G1 = fft2(g1FFT);
freq_G1 = circshift(freq_G1, [-rows/2, -cols/2, 0]); % shift back to center, for comparison's sake

由于我使用 sigma = 1.0 进行测试，我认为我应该得到两个等效的 2D 内核，因为如果 sigma = 1.0，那么 1.0/sigma = 1.0。因此，g1FFT 将等于 freq_G1。

但是，我没有。即使在归一化之后，它们也有不同的大小。有什么我遗漏的吗？

Answer 1

为简单起见，我将首先介绍一维信号的情况。对于多维情况，也可以进行类似的观察。

连续时间高斯信号的傅立叶变换本身就是一个高斯函数，如this table 所示。可以注意到，时域中的高斯越宽，频域中变换后的高斯越窄，对于 mu=0 和 sigma=1/sqrt(2π)（对应于 a=1/(2*sigma^ 2)=上面变换表中的π)，连续时间函数的傅里叶变换

将是类似的功能（仅发生变量变化）：

这很好，但这是针对连续时间信号的，我们对离散时间信号非常感兴趣。 不幸的是，正如wikipedia 中所指出的，通过对连续时间高斯函数进行采样而获得的核的离散傅里叶变换本身并不是采样的高斯函数。 幸运的是，这种关系通常仍然是近似正确的（无需过多详细说明，它需要时域内核足够宽但又不能太宽，以至于频域近似也足够宽以使关系也近似逆变换为真）。在这种情况下，离散时间信号的周期扩展（周期为N）的离散傅里叶变换

其中 mu=0 和 sigma=sqrt(N/2π) 可以用类似的函数来近似（直到一个比例因子和变量的变化）：

然后您可以修改 buildKernel 以分别支持沿行和列的不同标准差 sqrt(rows/2π) 和 sqrt(cols/2π)：

function result = buildKernel(rows, cols, mu, sigma)  

  if (length(mu)>1)
    mu_h = mu(1);
    mu_w = mu(2);
  else
    mu_h = mu;
    mu_w = mu;
  endif
  if (length(sigma)>1)
    sigma_h = sigma(1);
    sigma_w = sigma(2);
  else
    sigma_h = sigma;
    sigma_w = sigma;
  endif

  center_w = mu_w + floor(cols / 2);
  center_h = mu_h + floor(rows / 2);

  r = transpose(normpdf([0:rows-1],center_h,sigma_h));
  c = normpdf([0:cols-1],center_w,sigma_w);
  result = repmat(r * c, [1 1 3]);

  % normalize so that kernel sums to 1
  sumKernel = sum(result(:));
  result = result ./ sumKernel;    

end

您可以使用它来获取其 FFT 是其自身的缩放版本的内核。换句话说，使用获得的内核

g1FFTin = buildKernel(rows, cols, mu, [sqrt(rows/2/pi) sqrt(cols/2/pi)]);

这样freq_G1（在您发布的代码中计算得出）几乎等于g1FFTin * sqrt(rows*cols)。

最后考虑到您的意图实际上只是为了测试内核的 FFT 也是（近似）高斯分布，您可能希望将具有标准偏差 sigma 的任意内核的 FFT 与 另一个 进行比较> 直接在频域中计算的适当缩放的高斯核。换句话说，假设通过以下方式获得空间域内核：

g1FFTin = buildKernel(rows, cols, mu, sigma);

用相应的频域表示获得：

g1FFT   = circshift(g1FFTin, [rows/2, cols/2, 0]);
freq_G1 = fft2(g1FFT);
freq_G1 = circshift(freq_G1, [-rows/2, -cols/2, 0]);

然后freq_G1 可以与另一个在频域中直接计算的适当缩放的高斯核进行比较：

freq_G1_approx = (rows*cols/(2*pi*sigma^2))*buildKernel(rows, cols, mu, [rows cols]/(2*pi*sigma));