如何获得图像中使用的一维和二维空间滤波器的光谱?

Tae*_*oon 4 matlab signal-processing image-processing filter

我想知道如何在图像处理中获得空间1维和2维锐化滤波器的光谱.

锐化滤镜是:

[-1, 3, -1];
[-1, -1, -1; -1, 8, -1; -1, -1, -1];
Run Code Online (Sandbox Code Playgroud)

在MATLAB中,我该怎么做才能获得滤波器的频谱,或者如何获得这些滤波器的频率成分?

ray*_*ica 7

您可以使用Luis Mendo的帖子立即确定滤镜的频谱.二者必选其一freqz的一维过滤器,或freqz2为2D滤波器.请注意,使用freqz和的图表freqz2是跨越的标准化频率[-1,1].

但是,我想写这篇文章来赞美Luis Mendo,并向您展示频谱是如何得出的.


让我们从您的第一个过滤器开始:

h1 = [-1,3,-1];
Run Code Online (Sandbox Code Playgroud)

如果您还记得,1D 离散时间傅里叶变换(DTFT)定义为:


旁白:为什么DTFT而不是离散傅里叶变换(DFT)?

请特别注意,这是不同的离散傅立叶变换(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]是1D信号,并且omega是以弧度定义的角频率.因此,您可以将滤波器视为此1D信号,当您在空间域中进行滤波时,它与获取此信号,将其转换为频域并与频域中的其他输入信号进行乘法相同.

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

h = [-1    3  -1]
      ^    ^   ^
k =   -1   0   1
Run Code Online (Sandbox Code Playgroud)

因此,频域表示仅仅是具有复指数项的滤波器系数的加权和.代h[k]入傅里叶变换公式,我们得到:

如果你还记得欧拉的公式,我们有:

同理:

如果我们将两个方程式加在一起,我们可以重新排列并求解cos(x):

回到我们的1D滤波器的变换,我们可以做一些因子:

注意,域之间是[-pi,pi],因为1D傅立叶变换是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]);
Run Code Online (Sandbox Code Playgroud)

linspace产生从线性增加的阵列-pipi100分在范围之间.您可以通过指定手动linspace指示要生成多少点的第三个参数来覆盖此行为,但如果省略此参数,则默认值为100.

请注意,我已经使y-axis从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]
Run Code Online (Sandbox Code Playgroud)

由于形状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');
Run Code Online (Sandbox Code Playgroud)

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

我们得到:

在此输入图像描述

请注意,这两个维度都是从[-pi,pi].如果检查频谱,您将看到DC分量被置为0,这实际上是一个高通滤波器,而不是锐化滤波器.请参阅下面的说明.


小调

顺便说一句,您的2D滤镜定义不是 2D锐化滤镜.它只是一个2D拉普拉斯算子,它是一个边缘检测,因此是一个高通滤波器.它只检测边缘.如果要正确锐化图像,请确保将1添加到内核的中心,因此您真正想要的是:

h2 = [-1 -1 -1; -1 9 -1; -1 -1 -1];
Run Code Online (Sandbox Code Playgroud)

偏移量1将确保您保持原始信号不变,但也会在原始信号的顶部添加高通结果,从而锐化您的输入.

  • 这是......全面!! (2认同)