我需要在MATLAB中实现一个高通Butterworth滤波器，用于图像过滤。我已经实现了一个，但是看起来它不起作用。以下是我编写的代码，请问有人知道哪里出了问题吗？
n=1;
d=50;
A=1.5;
im=imread('imagex.jpg');
h=size(im,1);
w=size(im,2);
[x y]=meshgrid(-floor(w/2):floor(w-1/2),-floor(h/2):floor(h-1/2));
hhp=(1./(d./(x.^2+y.^2).^0.5).^(2*n));
image_2Dfilter=fftshift(fft2(im));
Image_butterworth=image_2Dfilter;
imshow(Image_butterworth);
ifftshow(Image_butterworth);
首先，没有名为ifftshow的命令。其次，你没有进行任何滤波处理。你所做的只是可视化图像的频谱。
在可视化频谱方面，你当前的方法非常危险。首先，你正在可视化每个空间频率分量上的系数，其本质上是复值。如果你想以对大多数人有意义的方式可视化频谱，最好查看幅度或相位。然而，由于这是一个巴特沃斯滤波器，所以最好将其应用于滤波器的幅度。
你可以使用abs函数来获取频谱的幅度。即使这样做，如果你直接在幅度上使用imshow，你会得到一个几乎除了中心外的所有位置都为零的可视化图像。这是因为DC分量非常大，其余频谱与之相比很小。
让我举个例子。这是图像处理工具箱中的摄影师图像：
im = imread('cameraman.tif');
figure;
imshow(im);
现在，让我们可视化频谱，并确保直流分量位于图像中心 - 您已经使用fftshift完成了此操作。 还应该将图像转换为double以确保数据的最佳精度。 此外，请确保应用abs来查找幅度：
fftim = fftshift(fft2(double(im)));
mag = abs(fftim);
figure;
imshow(mag, []);
正如您所看到的，由于我提到的原因，它并不是非常有用。更好地可视化图像的谱通常是将对数转换应用于谱。如果您想去除均值或平均值以便动态范围更适合显示，这也很有用。换句话说，您将添加1到幅度，然后对幅度应用对数，使得较高值可以逐渐减小。使用哪个基数并不重要，因此我将使用自然对数来封装log命令：
figure;
imshow(log(1 + mag), []);
现在好多了。现在我们将进入您的过滤机制。 您的Butterworth过滤器略有不正确。 坐标的meshgrid略有不对。 结束间隔处的-1操作需要移到外面：
[x y]=meshgrid(-floor(w/2):floor(w/2)-1,-floor(h/2):floor(h/2)-1);
记住，您正在定义关于图像中心的对称区间，而您最初的结果是不正确的。我还要提到这看起来像是一个高通滤波器，所以输出应该类似于边缘检测。另外，Butterworth高通滤波器的定义是不正确的。在频域中，滤波器的正确定义为:
D(u,v) 是图像在频域中心的距离，Do 是截止距离，而 B 是控制比例因子，控制截止距离处所需增益是多少。n 是滤波器的阶数。在您的情况下，Do 为 d = 50。在实践中，B = sqrt(2) - 1 ，以便在截止距离处，D(u,v) = 1 / sqrt(2) = 0.707，这是电子电路滤波器中经常看到的3 dB截止频率。有时你会看到 B 被设置为1以简化计算，但将其设置为 B = sqrt(2) - 1 是很常见的。
然而，您当前的代码并没有进行任何过滤。要在频域中进行滤波，只需将图像的频谱与滤波器本身的频谱相乘即可。这等效于空间域中的卷积。完成后，只需撤消在图像上执行的fftshift，取逆FFT，然后消除由于数值不精确导致的任何虚部。此外，让我们转换为uint8以确保我们尊重原始图像类型。
可以这样做:
%// Your code with meshgrid fix
n=1;
d=50;
h=size(im,1);
w=size(im,2);
fftim = fftshift(fft2(double(im)));
[x y]=meshgrid(-floor(w/2):floor(w/2)-1,-floor(h/2):floor(h/2)-1);
%hhp=(1./(d./(x.^2+y.^2).^0.5).^(2*n));
%%%%%%// New code
B = sqrt(2) - 1; %// Define B
D = sqrt(x.^2 + y.^2); %// Define distance to centre
hhp = 1 ./ (1 + B * ((d ./ D).^(2 * n)));
out_spec_centre = fftim .* hhp;
%// Uncentre spectrum
out_spec = ifftshift(out_spec_centre);
%// Inverse FFT, get real components, and cast
out = uint8(real(ifft2(out_spec)));
%// Show image
imshow(out);
如果您想查看滤波后的频谱是什么样子，只需执行以下操作：
figure;
imshow(log(1 + abs(out_spec_centre)), []);
我们得到：
这很有道理。您可以看到，在频谱的中心，与频谱的外缘相比稍微更暗一些。这是因为使用高通Butterworth滤波器时，您正在放大更高频率的项，并且它被可视化为更高的强度。
现在， out 包含您的过滤图像，我们最终得到：
看起来结果不错！然而，天真地将图像转换为 uint8 将负值截断为0，将大于255的正值截断为255。因为这是边缘检测，您想要检测负和正的变化...所以一个好主意是标准化输出，使其范围为[0,1]，然后在乘以255之后再进行类型转换为uint8。这样，图像中没有变化被可视化为灰色，负变化被可视化为黑色，正变化被可视化为白色....所以您可以这样做：
%// Your code with meshgrid fix
n=1;
d=50;
h=size(im,1);
w=size(im,2);
fftim = fftshift(fft2(double(im)));
[x y]=meshgrid(-floor(w/2):floor(w/2)-1,-floor(h/2):floor(h/2)-1);
%hhp=(1./(d./(x.^2+y.^2).^0.5).^(2*n));
%%%%%%// New code
B = sqrt(2) - 1; %// Define B
D = sqrt(x.^2 + y.^2); %// Define distance to centre
hhp = 1 ./ (1 + B * ((d ./ D).^(2 * n)));
out_spec_centre = fftim .* hhp;
%// Uncentre spectrum
out_spec = ifftshift(out_spec_centre);
%// Inverse FFT, get real components
out = real(ifft2(out_spec));
%// Normalize and cast
out = (out - min(out(:))) / (max(out(:)) - min(out(:)));
out = uint8(255*out);
%// Show image
imshow(out);
我们得到了这个：
我认为你应该以稍微不同的方式工作
n=1;
D0=50; % change the name for d0, d is usuaally the (u²+v²)⁽1/2)
A=1.5; % normally the amplitude is 1
im=imread('cameraman.jpg');
[M,N]=size(im); % is easy to get the h and w like this
% compute the 2d fourier transform in order to multiply
F=fft2(double(im));
% compute your filter and do the meshgrid for your matrix but it is M*n, and get only the real part
u=0:(M-1);
v=0:(N-1);
idx=find(u>M/2);
u(idx)=u(idx)-M;
idy=find(v>N/2);
v(idy)=v(idy)-N;
[V,U]=meshgrid(v,u);
D=sqrt(U.^2+V.^2);
H =A * (1./(1 + (D0./D).^(2*n)));
% multiply element by element
G=H.*F;
g=real(ifft2(double(G)));
subplot(1,2,1); imshow(im); title('Input image');
subplot(1,2,2); imshow(g,[ ]); title('filtered image');
D = (x.^2 + y.^2).^(0.5);，而是使用函数sqrt。 - Jonathan Hsqrt()函数。否则，这将调用通用的幂函数，并且不如使用sqrt精确。 - rayryeng