Matlab对灰度图像的频域进行高通滤波和低通滤波

Matlab对灰度图像的频域进行高通滤波和低通滤波,第1张

1. 要求

 对灰度图像进行离散傅里叶变换(Discrete Fourier Transfom, DFT)变换,在频域上分别使用理想的高通和低通滤波器进行滤波,显示滤波后的频域图像,以及逆离散傅里叶变换(Inverse Discrete Fourier Transfom, IDFT)变换后的空域图像,观察振铃现象。

2. 读取灰度图像

这里读取matlab自带的“摄影师”灰度图像

%% 读取图像
x = imread('cameraman.tif');
figure, imshow(x), title('原图')

 

3. 离散傅里叶变换

关于傅里叶变换的理解可以参考这篇文章

通俗讲解:图像傅里叶变换 - 知乎 (zhihu.com)https://zhuanlan.zhihu.com/p/99605178结合这篇文章,使用matlab实现离散傅里叶变化可以参考这篇文章

使用matlab对图像进行傅里叶变换 - 三山音 - 博客园 (cnblogs.com)https://www.cnblogs.com/sanshanyin/p/5426822.html这里使用的fft2是二维傅里叶变换,fftshift可以将将零频分量移到频谱中心,abs是绝对值和复数的模,如abs(3+4i)=5,imshow(I, [])不同于imshow(I),它是显示灰度图像 I,根据 I 中的像素值范围对显示进行转换。imshow 使用 [min(I(:)) max(I(:))] 作为显示范围。imshow 将 I 中的最小值显示为黑色,将最大值显示为白色。

%% 傅里叶变换
y = fft2(x);       % 二维傅里叶变换
y = fftshift(y);   % 将零频分量移到频谱中心
y_abs = abs(y);    % 对复数进行模运算 
z = log(y_abs+1);  % 映射到小范围区间
% 根据z中的像素值范围对显示进行转换,将z中的最小值显示为黑色,最大值显示为白色。
figure, imshow(z, []), title('频域');

  

4. 低通滤波

如果图像的宽度为col ,高度为row ,那么理想低通频域滤波器可形式化地描述为:

D_{0} \end{matrix}\right." src="https://latex.codecogs.com/gif.latex?LPF%28i%2Cj%29%3D%5Cleft%5C%7B%5Cbegin%7Bmatrix%7D%201%2C%20%26%20%5Csqrt%7B%28i-row/2%29%5E%7B2%7D+%28j-col/2%29%5E%7B2%7D%7D%20%5Cleq%20D_%7B0%7D%5C%5C%200%2C%20%26%20%5Csqrt%7B%28i-row/2%29%5E%7B2%7D+%28j-col/2%29%5E%7B2%7D%7D%20%3E%20D_%7B0%7D%20%5Cend%7Bmatrix%7D%5Cright." />

其中表示理想低通滤波器的截止频率,滤波器的频率域原点在频谱图像的中心处,在以截止频率为半径的圆形区域之内的滤镜元素值全部为1,而该圆之外的滤镜元素值全部为0。理想低通滤波器的频率特性在截止频率处十分陡峭,无法用硬件实现,这也是使用者称之为理想的原因,但其软件编程的模拟实现较为简单。

由此,我们可以制作一个低通滤波器LPF,LPF初始化为一个全部为1的矩阵,接着遍历y,对其中的每个点(i,j),判断其频率是否高于截止频率,若高于截止频率,则将低通滤波器中对应的点置为0,若小于等于截止频率,则不做处理。

这里y.*LPF,是将y中的所有点分别乘以LPF中的所有点,对于y(i,j),如果其频率高于截止频率,则其对应位置的LPF(i,j)=0,两者想相结果为0,即过滤了高频;如果其频率小于等于截止频率,则其对应位置的LPF(i,j)=1,两者相乘结果仍为y(i,j),即保留了低频。

%% 低通滤波器
LPF = ones(size(y));        % 初始化低通滤波器
threshold_low = 50;         % 设置截至频率
[row, col] = size(y);       % 得到图像y的高度和宽度
for i = 1: row
    for j = 1: col
        % 计算频率
        if sqrt((i-row/2)^2+(j-col/2)^2) > threshold_low
            % 将低通滤波器中高于截止频率的频率变为0
            LPF(i,j) = 0;
        end
    end
end
y_LPF = y.*LPF;             % 低通滤波
y_LPF = ifftshift(y_LPF);   % 逆零频平移
x_LPF = ifft2(y_LPF);       % 傅里叶逆变换
x_LPF = uint8(real(x_LPF)); % 转换成实数,并将double转成uint8
figure, imshow(x_LPF), title('低通滤波')

更一般地,低通滤波器仅仅是一个形式,每次为了得到指定频率的滤波器(低通也好,高通也罢),如果都需要遍历row*col次,效率也是很低的。

本质上,理想低通滤波器和理想高通滤波器都是基于频率的,我们可以先计算出每个点对应的频率,将其存放在一个频率矩阵(频率本质上也是距离)distance中,接下来,对于任意一个给定的截止频率,我们只需要对distance矩阵进行处理即可。

(上面这段的代码和下面的两段代码在效果上是相同的)

%% 计算频率,频率 = 当前点(i,j) 到 中心点(row/2, col/2)的距离
[row, col] = size(y);
distance = zeros(row, col);
for i = 1: row
    for j = 1: col
        % 计算向量的模长,即两点间距离
        distance(i,j) = norm([i-row/2 j-col/2]);
    end
end

ifftshift是逆零频平移,是fftshift的逆变换,ifft2是二维快速傅里叶逆变换,是fft2的逆变换,

real用来得到复数的实部,返回的数据类型是double,一般灰度图像的数据类型是uint8,所以需要使用uint8进行转换。

%% 低通滤波 (Low Pass Filtering, LPF)
y_LPF = y;                  % 复制一个正常的y
threshold_low = 50;         % 设置阈值
% 过滤频率高于阈值的点
y_LPF(distance > threshold_low) = 0;
y_LPF = ifftshift(y_LPF);   % 逆零频平移
x_LPF = ifft2(y_LPF);       % 傅里叶逆变换
x_LPF = uint8(real(x_LPF)); % 转换成实数,并将double转成uint8
figure, imshow(x_LPF), title('低通滤波')

 

5. 高通滤波

高通滤波同理

%% 高通滤波 (HPF, High Pass Filtering)
y_HPF = y;
threshold_high = 30;
y_HPF(distance < threshold_high) = 0;
y_HPF = ifftshift(y_HPF);
x_HPF = ifft2(y_HPF);
x_HPF = uint8(real(x_HPF));
figure, imshow(x_HPF), title('高通滤波')

6. 完整代码

关注公众号TobleroneWind,回复“数字图像处理”即可获取完整代码。

欢迎分享,转载请注明来源:内存溢出

原文地址: http://outofmemory.cn/web/1295096.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2022-06-10
下一篇 2022-06-10

发表评论

登录后才能评论

评论列表(0条)

保存