低通滤波器的问题

Question

低通滤波器的问题

matlabimage-processingsignal-processingdftlowpass-filter

4

我需要在一张图片上应用一个简单的低通滤波器，但是这需要在频域中完成。然而最终结果是一张非常嘈杂的图片，这表明我的处理过程出现了问题。

我的思路是：

1）将图像（Testpatt）居中并应用2D FFT

tpa1 = size(Testpatt, 1); tpa2 = size(Testpatt, 2);
dTestpatt = double(Testpatt);


for i = 1:tpa1
    for j = 1:tpa2

        dTestpatt(i,j) = (dTestpatt(i,j))*(-1)^(i + j);
    end
end

fftdTestpatt = fft2(dTestpatt);

2) 生成低通滤波器并将其与傅里叶变换后的图像相乘（注意：低通滤波器需要是半径为Do的1的圆）。

lowpfilter = zeros(tpa1, tpa2);
Do = 120;

for i = 1:tpa1
    for j = 1:tpa2


       if sqrt((i - tpa1/2)^2 + (j - tpa2/2)^2) <= Do
           lowpfilter(i,j) = 1;
       end

   end
end


filteredmultip = lowpfilter*fftdTestpatt;

3）提取逆二维FFT的实部并还原中心化。

filteredimagecent = real(ifft2(filteredmultip));


for i = 1:tpa1
    for j = 1:tpa2

        filteredimage(i,j) = (filteredimagecent(i,j))*(-1)^(i + j);
    end
end

任何帮助或意见都将不胜感激。

- Rallad

请问您能否添加图片（或类似的）？ - Finn

任何图像都可以，效果应该是一种容易识别的模糊效果。 - Rallad

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- rayryeng · Accepted Answer

我很惊讶这没有给你一个错误，但你在频域执行的是矩阵乘法而不是按元素相乘。记住，在频域中进行过滤需要对两个变换信号进行逐元素相乘，以执行等效的空间域卷积操作。你只需要更改乘法语句，使其成为按元素的等效语句即可：

filteredmultip = lowpfilter.*fftdTestpatt;

此外，在显示图像之前，确保将图像转换回与原始图像相同的数据类型。例如，如果您的图像是uint8，则需要使用im2uint8进行转换。这主要是为了显示和写入图像到文件方便。如果您像在代码中所做的那样将其保留为double，则显示图像将被视为大部分白色，因为显示类型为double的图像假定值的范围从[0,1]。

顺便说一句，我怀疑您没有收到错误的原因是因为您的图像和过滤器是正方形的，因此在进行矩阵乘法时维度是有效的。如果您决定对非正方形图像执行此操作，则肯定会出现错误。

警告 - 使用理想低通滤波器

您正在实施的是使用理想低通滤波器进行滤波，当您进行低通滤波时，会出现振铃效应。原因是这直接来自信号处理理论。实现频域硬边缘需要大量（或者说无限...）的空间域正弦波。请记住，FFT是一个正弦波分解的信号。当您使用此低通滤波器来过滤图像时，这被视为重建图像中的振铃，因为要创建硬边缘需要大量正弦波求和（因此呈现出振铃）。

为了演示这些效果，让我们以示例图像为例。我将使用作为图像处理工具箱一部分的摄影师图像。如果我使用半径Do = 40并运行您的代码（已更正），则这是原始图像和过滤后的结果：

如您所见，那非常糟糕。振铃来自频域滤波器的硬边缘。人们倾向于随着距离中心的增加而逐渐降低幅度。一个很好的例子是高斯模糊。您需要做的是定义标准差，然后创建一个与您的图像大小相同的掩模，该掩模代表一个2D高斯函数。

请记住，2D高斯可以表示为：

^{来源：维基百科}

你只需循环所有像素坐标并计算上述方程。我没有在前面乘以缩放因子，因为你实际上不需要这样做。因此，你可以使用这段代码生成高斯掩模，这等效于使用两个for循环，但更加向量化。我定义了一个跨越与你的图像相同大小的2D坐标网格，然后在图像中的每个位置运行该方程。当然，我们需要将核心放在中心，因此你必须通过图像中心偏移，就像你之前的低通滤波代码一样。我还决定将你的Do变量用作标准差。

Do = 40;
[X,Y] = meshgrid(1:tpa2, 1:tpa1);
lowpfilter = exp(-((X - (tpa2/2)).^2 + (Y - (tpa1/2)).^2) / (2*Do*Do));

因此，我们现在得到了这个滤波器响应：

... 当我进行滤波时，我获得：

与原始图像进行比较：

正如您所看到的，模糊效果更好。没有振铃。因此，在过滤图像时，请确保避免使用过滤器中的硬边缘，因为这将在空间域中产生振铃伪影。

一些更多建议

我有几个建议可以让您使此代码快速运行。

避免使用循环来居中图像

正如您从信号处理理论中已经知道的那样，如果您通过(-1)^(i+j)对图像中的像素值进行预先乘法，其中(i,j)是要过滤的像素的空间位置，这将在频率域中居中您的图像。我建议避免使用循环来执行此操作并首先对FFT进行中心化。您可以使用一个名为fftshift的函数为您执行中心化。首先对图像进行FFT，然后调用此函数：

fftdTestpatt = fft2(dTestpatt);
fftdTestpatt = fftshift(fftdTestpatt); % Centre the FFT

避免使用循环来生成过滤器

我也建议避免使用循环来生成您的过滤器。特别是对于具有理想过滤器的代码，请将您的代码替换为以下代码，该代码遵循与高斯过滤器相同的逻辑。我们还可以删除sqrt操作，并确保您正在与半径的平方进行比较，以避免任何不必要的计算：

[X,Y] = meshgrid(1:tpa2, 1:tpa1);
lowpfilter = (X - tpa1/2).^2 + (Y - tpa2/2).^2 <= Do^2;
lowpfilter = double(lowpfilter);

在以上代码中特别注意逐元素的幂运算（.^）。最后一条语句很重要，因为我们需要将结果转换回double类型，因为首先生成过滤器实际上会给你一个logical类型的结果。我们需要double类型以确保FFT的精度得到尊重。

避免使用循环来取消图像滤波后的中心化

完成滤波后，再次乘以(-1)^(i+j)来取消图像的中心化。使用相关的ifftshift函数在使用FFT进行滤波后取消中心化图像：

filteredmultip = ifftshift(filteredmultip); % Uncentre first
filteredimage = real(ifft2(filteredmultip));