使用 FFT 属性求二维函数的导数
Finding the derivative of a 2D function using FFT properties
我正在尝试使用 FFT 属性来获取二维函数的第 i 个导数 - 特别是二维高斯函数。我也想在 MATLAB 中以数字方式执行此操作。
实际上，我不知道我在做什么，但我在网上看了很多，所有的 MATLAB 帮助，似乎没有什么帮助我，所以我打算在这里问
这是我目前的情况：
h = fspecial('gaussian', size(imr), 10);
imshow(h, []);
G = fft2(h);
这是大小为 imr 或 512 x 512 的高斯函数的 0 阶导数。我知道我必须进行逆运算并乘以某些东西，但我不知道是什么。
任何人都可以给我提供有关使用 FFT 属性获取二维高斯函数的第 i 个导数的线索吗？
提前致谢
免责声明 - 2015 年 9 月 16 日编辑
这个 post 由于我对频域微分运算的理解存在一个小而根本的缺陷，因此发生了重大变化。自上次迭代以来，其中很多 post 已更改。
我要感谢 Luis Mendo 帮助我调试为什么除了 post 的原始版本中出现的这张图片之外，这对其他一些图像不起作用。
对于离散时间数据，没有导数这样的东西。导数是一种只存在于连续时间中的分析工具。在离散时间，我们只能使用差异来近似。因此，我假设您指的是差分运算而不是导数。导数最常见的近似之一是使用前向差分运算。具体地，对于一维离散序列x[n]，应用前向差分运算的输出序列y[n]定义为：
y[n] = x[n+1] - x[n], for n = 1, 2, ..., M-1
M 是离散序列中的样本总数。由于前向差分运算的性质，自然会少一个值。
要在频域中完成同样的事情，您必须使用离散傅立叶变换属性中的 shift theorem。具体来说，给定信号 x[n] 的 DFT 版本，即 X(k)，以下 属性 将时域和频域联系在一起：
这里，i表示复数，也就是-1的平方根，k是频域中的angular频率。 F^{-1}是信号X(k)的inverse傅里叶变换，也就是时域信号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 space 中的位置一致。按照上面的公式，进入 2D 非常简单：
这里，u和v表示二维离散差分运算的空间坐标y[u,v]。 U和V是二维频率分量，M和N是二维信号的列数和行数。请特别注意，您实际上是在进行水平差分运算，然后再进行垂直差分运算。具体来说，就是对每一行分别进行水平差分运算，然后对每一列分别进行垂直差分运算。事实上，顺序并不重要。您可以按照您选择的任何顺序先做一个再做另一个。需要注意的重要一点是，信号的二维傅里叶变换保持不变，您只是将每个值乘以某个复杂常数。但是，您需要确保删除结果的最后一行和最后一列，因为由于我们谈到的差异操作观察，每次操作都会导致一个信号比每行和每列的原始长度少一个点.
假设你已经有了函数的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] 之间。在显示两个域中导数运算的等价性时，我们需要考虑到这一点。此外，为了使事情变得更简单，如果您 移动 频率分量，使它们出现在图像的 centre 中，这也是明智的.执行 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 值，看看会得到什么。这个 属性 我在实践中没有看到经常使用，但知道它肯定很好。
我正在尝试使用 FFT 属性来获取二维函数的第 i 个导数 - 特别是二维高斯函数。我也想在 MATLAB 中以数字方式执行此操作。
实际上，我不知道我在做什么，但我在网上看了很多，所有的 MATLAB 帮助，似乎没有什么帮助我，所以我打算在这里问
这是我目前的情况：
h = fspecial('gaussian', size(imr), 10);
imshow(h, []);
G = fft2(h);
这是大小为 imr 或 512 x 512 的高斯函数的 0 阶导数。我知道我必须进行逆运算并乘以某些东西，但我不知道是什么。
任何人都可以给我提供有关使用 FFT 属性获取二维高斯函数的第 i 个导数的线索吗？
提前致谢
免责声明 - 2015 年 9 月 16 日编辑
这个 post 由于我对频域微分运算的理解存在一个小而根本的缺陷，因此发生了重大变化。自上次迭代以来，其中很多 post 已更改。
我要感谢 Luis Mendo 帮助我调试为什么除了 post 的原始版本中出现的这张图片之外，这对其他一些图像不起作用。
对于离散时间数据，没有导数这样的东西。导数是一种只存在于连续时间中的分析工具。在离散时间，我们只能使用差异来近似。因此，我假设您指的是差分运算而不是导数。导数最常见的近似之一是使用前向差分运算。具体地，对于一维离散序列x[n]，应用前向差分运算的输出序列y[n]定义为：
y[n] = x[n+1] - x[n], for n = 1, 2, ..., M-1
M 是离散序列中的样本总数。由于前向差分运算的性质，自然会少一个值。
要在频域中完成同样的事情，您必须使用离散傅立叶变换属性中的 shift theorem。具体来说，给定信号 x[n] 的 DFT 版本，即 X(k)，以下 属性 将时域和频域联系在一起：
这里，i表示复数，也就是-1的平方根，k是频域中的angular频率。 F^{-1}是信号X(k)的inverse傅里叶变换，也就是时域信号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 space 中的位置一致。按照上面的公式，进入 2D 非常简单：
这里，u和v表示二维离散差分运算的空间坐标y[u,v]。 U和V是二维频率分量，M和N是二维信号的列数和行数。请特别注意，您实际上是在进行水平差分运算，然后再进行垂直差分运算。具体来说，就是对每一行分别进行水平差分运算，然后对每一列分别进行垂直差分运算。事实上，顺序并不重要。您可以按照您选择的任何顺序先做一个再做另一个。需要注意的重要一点是，信号的二维傅里叶变换保持不变，您只是将每个值乘以某个复杂常数。但是，您需要确保删除结果的最后一行和最后一列，因为由于我们谈到的差异操作观察，每次操作都会导致一个信号比每行和每列的原始长度少一个点.
假设你已经有了函数的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] 之间。在显示两个域中导数运算的等价性时，我们需要考虑到这一点。此外，为了使事情变得更简单，如果您 移动 频率分量，使它们出现在图像的 centre 中，这也是明智的.执行 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 值，看看会得到什么。这个 属性 我在实践中没有看到经常使用，但知道它肯定很好。