matlab - 使用FFT属性查找2D函数的导数

Question

我正在尝试使用FFT属性来获取2D函数(尤其是2D高斯函数)的i阶导数。我也想在MATLAB中以数字方式执行此操作。

实际上，我对正在做的事情一无所知，但是我在Internet上已经读了很多书，并且在MATLAB的所有帮助下都读了很多书，似乎没有什么对我有帮助，所以我将在这里提出问题。

这是我到目前为止的内容:

h = fspecial('gaussian', size(imr), 10);
imshow(h, []);
G = fft2(h);

这是 imr或512 x 512大小的高斯的0阶导数。我知道我必须进行逆运算并乘以某物，但我不知道该怎么做。

有人能给我一个有关使用FFT属性获得2D高斯函数的i阶导数的线索吗？

提前致谢

Answer 1

免责声明-2015年9月16日编辑

由于我对频域微分运算的理解存在一个小而根本的缺陷，因此该帖子已发生了显着变化。自上次迭代以来，许多帖子已更改。

我要感谢Luis Mendo帮助我调试了为什么这对其他图片不起作用，但该帖子的原始版本中没有提供其他图片。

当涉及离散时间数据时，不存在任何派生的东西。导数是仅在连续时间中存在的分析工具。在离散时间内，我们只能通过使用差异来近似。因此，我将假设您是指差运算而不是导数。导数最常见的近似之一是使用前向差分运算。具体而言，对于一维离散序列x[n]，将应用前向差分运算的输出序列y[n]定义为:

y[n] = x[n+1] - x[n], for n = 1, 2, ..., M-1

M是离散序列中的样本总数。自然地，由于前向差操作的性质，值会少一个。

要在频域中完成相同的操作，您必须使用离散傅立叶变换属性中的 shift theorem。具体地说，给定信号 x[n]的DFTed版本 X(k)，以下属性将时域和频域联系在一起:

在此， i表示复数或-1的平方根，并且 k是频域中的角频率。 F^{-1}是信号 X(k)的 逆傅立叶变换，它是时域信号 x[n]的傅立叶变换。 M也是信号 x[n]的长度，而 m是您要移动信号的偏移量。因此，这就是说，如果要通过按 m位置移位来计算移位操作，则只需进行傅立叶变换，将每个分量逐个乘以 exp(-i*2*pi*m*k/M)，其中 k是指向傅立叶点的索引变换并对该中间​​结果进行傅立叶逆变换。

特别要注意的是，傅立叶变换的最后一个元素是没有意义的，因为差分运算导致信号中的元素减少了一个，因此从该结果中删除最后一个元素非常重要。

有了上面提到的内容，为了计算导数的近似值，我们需要找到 y[n]或 Y(k)的傅立叶变换，并且可以像这样进行计算:

请注意，移位是 -1，即 x[n+1] = x[n-(-1)]。因此，您将计算傅立叶变换，将每个对应的值乘以 exp(i*2*pi*k/M) - 1并求逆。

进入2D并不是一件容易的事。在这里，我假设我们正在水平和垂直两个方向上执行差异操作。记住，我们有两个代表水平和垂直空间坐标的变量。我们也将其称为 空间域，而不是时域，因为变量是根据我们在2D空间中的位置而定的。按照上述公式，进入2D非常简单:

在此， u和 v表示2D离散差操作 y[u,v]的空间坐标。 U和 V是2D频率分量，而 M和 N是2D信号的列数和行数。请特别注意您实际上是在进行水平差分操作，然后再进行垂直差分操作。具体来说，您将对每一行分别进行水平差分操作，然后对每一列分别进行垂直差分操作。实际上，顺序无关紧要。您可以按照选择的顺序依次执行一个。需要注意的重要一点是，信号的2D傅立叶变换保持不变，而您只是将每个值乘以某个复数常数。但是，您需要确保删除结果的最后一行和最后一列，因为由于我们谈到的差异操作观察，每次操作所产生的信号比每行和每列的原始长度小一个点。

假设您已经具有该功能的FFT，则只需获取每个2D空间坐标并乘以 (exp(i*2*pi*U/M) - 1)*(exp(i*2*pi*V/N) - 1)即可。对于存储在 (U,V)处的FFT的每个2D分量，我们使用这些相同的频率坐标，并将该位置乘以前面两个东西的乘积。之后，您将进行逆FFT，并将其与实际的空间域导数进行比较。我们应该看到它们是相同的。

请注意，执行2D FFT时，频率是 归一化的，因此它们在两个维度上都在 [-1,1]之间。当显示两个域中的导数运算的等效性时，我们需要考虑到这一点。另外，为简化起见，如果您要 移动频率分量，使其出现在图像的 中心中，这也是明智的做法。执行2D FFT时，将分量放置为使DC分量或原点位于图像的左上角。让我们使用 fftshift 将所有内容移到中心。

现在，让我们从一些小事情开始。为了这个过程，我们假设每个维度的大小都是奇数，以使通过 fftshift的频率移动变得更加容易和明确。假设我们有一个101 x 101高斯滤波器，标准偏差为10，然后对此进行FFT，然后进行 fftshift。所以:

h = fspecial('gaussian', [101 101], 10);
hf = fftshift(fft2(h));

如果要找到导数，则定义一个频点网格，该频点网格跨越频域中图像的大小，从 0作为原点开始，现在是 中心。我们可以通过 meshgrid 做到这一点，所以:

[U,V] = meshgrid(-50:50, -50:50);

通常，此水平和垂直坐标的范围是水平的 [-floor(cols/2),floor(cols/2)]和垂直的 [-floor(rows/2),floor(rows/2)]。请注意，上面的网格与图像具有相同的尺寸，但是中心在 (0,0)上，我们在两个尺寸上的范围都在-50到50之间。

现在，在两个方向上在频域中执行差分操作。让我们对 x和 y都对 这样做，所以这两个方向的一阶导数:

hf2 = (exp(1i*2*pi*U/101) - 1).*(exp(1i*2*pi*V/101) - 1).*hf;

请注意，由于FFT的频率已归一化，因此我们必须通过除以101来对频率进行归一化，以将频率归一化。 101 x 101是由于图像的尺寸。如果要在时域中返回，只需进行逆FFT。但是，在进行逆FFT之前，我们需要 撤消与 ifftshift 所做的转换。我们还使用 real消除了由于数值错误而导致的任何残留复数值分量。您会看到输出的值为复数，但是虚部的数量非常小，我们可以忽略它们。因此:

out_freq = real(ifft2(ifftshift(hf2)));

记住也要删除最后一行和最后一列:

out_freq = out_freq(1:end-1,1:end-1);

如果您想将此输出与过滤器 h的时域差值操作进行比较，我们可以通过在行上使用 diff 来精确地做到这一点，然后使用此结果遍历各列。因此:

out_spatial = diff(diff(h, 1, 1), 1, 2);

第二个参数表示差异运算的顺序。这是一阶差分，因此水平和垂直差分运算都设置为1。第三个参数告诉您要在哪个尺寸上进行差分。我们首先在列上执行此操作，然后在行上应用中间结果。然后，我们可以显示它们的外观:

figure;
subplot(1,2,1);
imshow(out_spatial, []);
title('Derivative from spatial domain');
subplot(1,2,2);
imshow(out_freq, []);
title('Derivative from frequency domain');

我们得到:



凉爽的!这与我对高斯导数的了解一致。

奖金

作为一点奖金，我想向您指出康奈尔大学的一些精彩幻灯片: http://www.cs.cornell.edu/courses/cs6670/2011sp/lectures/lec02_filter.pdf

具体来说，如果您要在3D模式下查看高斯图像，它就是水平和垂直坐标，并且 Z值是3D模式下高斯图像的高度:

因此，如果我们独立地采用两个方向的导数，就会发生这种情况。顶部显示您在空间上发生的事情，底部显示您在鸟瞰图上方直接观察到的情况:

表面的最负部分用黑色绘制，表面的最正部分用白色绘制，其他所有部分都在强度上按比例缩放...因此，如果我们采用一个方向的空间导数，然后使用此中间结果并采用另一个方向的空间导数，您会看到我们将得到上面所看到的结果。尝试不同的 m和 n值，看看会得到什么。我从未见过的该属性在实践中经常使用，但一定很高兴知道。