Matlab - 高斯函数的FFT - 等效性

Question

Matlab - 高斯函数的FFT - 等效性

4

简单问题：

我在Matlab中以一定的分辨率绘制了一个二维高斯函数。我测试方差或 sigma = 1.0 的情况。我想将其与FFT（高斯函数）的结果进行比较，这应该得到另一个方差为(1./sigma) 的高斯函数。由于我正在测试 sigma = 1.0，我认为我应该得到两个相等的二维核。

例如：

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（高斯函数）的结果进行比较，这应该会产生另一个方差为(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。

但是，事实并非如此。即使在归一化之后，它们也具有不同的幅度。我错过了什么吗？

- kdottiemo

请提供所有能够重现错误的代码。您说的高斯傅里叶变换确实会得到另一个高斯函数，但是如果没有代码展示您的操作，我们无法判断您是否正确地进行了操作。 - rayryeng

非常简单，只是添加了代码。 - kdottiemo

为什么你的内核是一个三维矩阵？你说它是一个二维高斯函数，但有三个通道。 - rayryeng

1

你可能已经发现了，但是你不能直接将连续时间傅里叶变换对应到离散时间/DFT上。http://users.ece.gatech.edu/mrichard/Gaussian%20FT%20and%20random%20process.pdf 显示出离散情况下的等效性取决于标准差以样本为单位。 - Ahmed Fasih

你是否知道fftshift函数可以用来移动fft操作的结果？ - rubenvb

显示剩余3条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- SleuthEye · Answer 1

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

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

enter image description here

将会是类似的函数（只有变量的更改发生了）：

enter image description here

这都很好，但是这是针对连续时间信号的，我们真正感兴趣的是离散时间信号。不幸的是，正如wikipedia上也指出的那样，通过对连续时间高斯函数进行采样得到的核的离散傅里叶变换本身并不是一个采样的高斯函数。幸运的是，这种关系通常仍然近似成立（不过不详细讨论的话，需要时域核足够宽但又不能太宽，以使频域近似也足够宽，从而使逆变换的关系也近似成立）。在这种情况下，周期为N的离散时间信号的周期扩展的离散傅里叶变换

enter image description here

当mu=0且sigma=sqrt(N/2π)时，可以通过类似的函数（经过缩放因子和变量更改）进行近似：

enter image description here

你可以修改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));