matlab 查看函数属性,matlab – 使用FFT属性查找2D函数的导数

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

由于我对频域中微分运算的理解存在一个小但基本的缺陷,这篇文章发生了重大变化.自上次迭代以来,很多帖子都发生了变化.

我要感谢Luis Mendo帮助我调试为什么这不适用于除了这个帖子的原始版本中提供的其他图像之外的其他图像.

当涉及离散时间数据时,没有衍生物这样的东西.衍生物是一种仅在连续时间内存在的分析工具.在离散时间,我们只能通过使用差异来近似.因此,我将假设您的意思是差异操作而不是衍生操作.导数的最常见近似之一是使用前向差分运算.具体地,对于1D离散序列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信号的列数和行数.特别注意您实际上正在进行水平差异操作,然后进行垂直差分操作.具体来说,您将独立地为每一行执行水平差异操作,然后分别对每列执行垂直差异操作.实际上,顺序并不重要.你可以按照你选择的顺序做一个跟随另一个.需要注意的是,信号的二维傅立叶变换保持不变,您只需将每个值乘以一些复数常数即可.但是,您需要确保删除结果的最后一行和最后一列,因为每次操作都会产生比每行和每列的原始长度小一个点的信号,这是由于我们所讨论的差异操作观察.

假设你已经有了函数的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将所有内容转移到中心.

现在,让我们从一些小事做起.为了这个过程,我们假设每个维度的大小是奇数,以允许通过ffthift移动频率更容易和明确.假设我们有一个101 x 101高斯滤波器,标准差为10,我们采用FFT,然后进行移位.所以：

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');

我们得到：

凉！这与我所知的高斯导数一致.

奖金

具体来说,这就是高斯看起来如果你在3D中看它,我们有水平和垂直坐标,Z值是3D中高斯的高度：

所以,如果我们独立地采用两个方向的导数就会发生这种情况.顶部显示了空间发生的情况,底部显示了如果我们直接在鸟瞰图上看到这一点会发生什么：

表面的最负面部分是用黑色绘制的,表面的最正面部分是用白色绘制的,其他所有部分都在强度方面进行缩放……所以,如果我们采用一个方向的空间导数,然后使用这个中间结果,并采取另一个方向的空间导数,你会看到我们会得到我们上面看到的.玩弄不同的m和n值,看看你得到了什么.我在实践中经常使用的这个属性,但它确实很有用.