Matlab

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）上面变换中的=πtable），连续时间函数的傅立叶变换

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

这很好，但这是针对连续时间信号的，我们对离散时间信号非常感兴趣。不幸的是，正如 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));

Matlab - 高斯的 FFT - 等效

Matlab - FFT of Gaussian - Equivalency

fft

gaussian