matlab - 如何获得图像中使用的1D和2D空间滤波器的光谱？

Question

我想知道如何在图像处理中获得空间一维和二维锐化滤镜的光谱。

锐化滤镜为：

[-1, 3, -1];
[-1, -1, -1; -1, 8, -1; -1, -1, -1];

在MATLAB中，应该怎么做才能获得滤波器的频谱，或者如何获得这些滤波器的频率分量？

Answer 1

您可以使用Luis Mendo的帖子立即确定滤波器的频谱。将freqz用于1D过滤器，或freqz2用于2D过滤器。请注意，使用freqz和freqz2的图是根据跨越[-1,1]的归一化频率表示的。

但是，我想写这篇文章来赞扬路易斯·门多（Luis Mendo），并向您展示频谱的推导方式。



让我们从第一个过滤器开始：

h1 = [-1,3,-1];

如果您还记得，一维Discrete-Time Fourier Transform (DTFT)定义为：





除了：为什么使用DTFT而不使用离散傅立叶变换（DFT）？

请特别注意这与Discrete Fourier Transform (DFT)不同。它们之间的区别在于DTFT是应用于离散信号的传统傅里叶变换。从连续的角度来看，我们应用正常的傅立叶变换，其中频谱在频率上是连续的。这里对于离散信号本质上是相同的，但是该变换被应用于离散信号。输出在频率上也是连续的，另外还有频谱是2*pi周期性的约束。对于DFT，它实质上是频域中DTFT的采样版本。因此，它们之间的核心区别在于，在频域中，一个是连续的，而另一个是连续副本的采样版本。我们需要对DTFT进行采样的原因是，您可以将光谱实际存储在计算机上，并让程序处理光谱（即MATLAB）。

您的问题是确定频谱是什么，并且从理论角度来讲，我将向您介绍DTFT。当我们实际绘制频谱时，无论如何我们实际上都在采样DTFT，因此当我们看到频谱时，您将可视化DFT（或多或少！）。

可以在DSP StackExchange上找到有关它们之间差异的很好解释：https://dsp.stackexchange.com/questions/16586/difference-between-discrete-time-fourier-transform-and-discrete-fourier-transfor



回到您的问题

对于DTFT的定义，h[k]是一维信号，omega是以弧度定义的角频率。因此，您可以将滤波器视为此一维信号，并且在空间域中进行滤波时，与获取该信号，将其转换为频域并与频域中的另一个输入信号相乘相同。

因此，如果我们认为您的滤波器是对称的，则将值3定义为中心点。因此，您可以将h[k]视为：

h = [-1    3  -1]
      ^    ^   ^
k =   -1   0   1

因此，频域表示只是带复指数项的滤波器系数的加权和。将h[k]代入傅立叶变换公式，我们得到：



如果您回想起Euler's formula，我们有：



类似地：



如果将两个方程式相加，则可以重新排列并求解cos(x)：



回到我们对一维滤波器的转换，我们可以做一些分解：



请注意，域在[-pi,pi]之间，因为一维傅立叶变换是2*pi对称的。因此，如果要显示光谱，只需使用上述域并使用我刚刚创建的方程式绘制绝对值即可，因为光谱的幅值始终为正：

w = linspace(-pi,pi);
h = abs(3 - 2*cos(w));
plot(w,h);
title('Magnitude of 1D spectrum');
axis([-pi, pi, 0, 5]);

linspace生成一个从-pi到pi的线性增加的数组，范围之间有100个点。您可以通过指定第三个参数来覆盖此行为，该参数手动告诉linspace您要生成多少个点，但是如果省略此参数，则默认值为100。

请注意，我已经使y轴从0开始延伸，因此您可以看到曲线的起始位置，然后上升到5，因为这是h可能的最大值。这是我们得到的：



实际上，上述频谱看起来确实像是一个高通滤波器，并且由于DC频率（w = 0）的幅度为1，因此您实际上是在高通滤波结果的基础上加上原始信号，从而“ ”您的信号。



您可以对2D情况进行相同的处理，尽管会更加复杂。在2D情况下，离散傅立叶变换定义为：



我们需要考虑两个独立变量，并且我们将拥有2D空间频率。 w1和k1将沿2D信号的行进行操作，而w2和k2将沿2D信号的行进行操作。

考虑到您的面具：

h2 = [-1 -1 -1; -1 8 -1; -1 -1 -1]

由于h2的形状是对称的，因此8的值将是位置(w1,w2) = (0,0)。因此，当我们使用以上公式计算频谱时，我们得到：



我将通过简化为您省去麻烦，但是通过重新安排并利用欧拉公式，我们得到：



注意：



我使用上面的属性来简化。对于2D情况，我们将定义两个维度的坐标meshgrid，均在[-pi,pi]之间，并使用surf将其绘制在表面图上。切记取滤波器的绝对值来显示幅度：

[w1,w2] = meshgrid(linspace(-pi,pi));
h = abs(8 - 2*cos(w1+w2) - 2*cos(w1) - 2*cos(w2) - 2*cos(w1-w2));
surf(w1,w2,h);
title('2D spectrum of filter');

w1和w2提供在阵列的每个相应空间位置处水平和垂直定义的频率的唯一组合。 w1和w2实际上将是二维的。对于w1和w2的每个唯一组合，我们确定光谱的幅度。一旦计算出这些，就可以将所有这些信息放到一个漂亮的三维表面图中。

因此，我们得到：



请注意，这两个维度均跨越[-pi,pi]。如果检查频谱，将会看到DC分量归零，这实际上是一个高通滤波器，而不是锐化滤波器。请参阅下面的注释。



小音符

顺便说一句，您的2D滤镜定义不是2D锐化滤镜。它只是2D拉普拉斯算子，它是边缘检测，因此是高通滤波器。它仅检测边缘。如果要适当锐化图像，请确保将1加到内核的中心，这样，您真正需要的是：

h2 = [-1 -1 -1; -1 9 -1; -1 -1 -1];

偏移量1将确保您保持原始信号不变，而且还将高通结果添加到原始信号之上，从而使您的输入更加清晰。