巴特沃斯滤波原理及代码实现（matlab详细过程版）

一、算法原理

1、原理概述

巴特沃斯滤波器（Butterworth filter）是一种连续衰减的滤波器，所以也被称为最大平坦滤波器，在该滤波器不会出现太大陡峭的变化。其特点是在通频带内呈现出最大限度的平坦的频率响应曲线，没有纹波，同时在阻频带内则逐渐下降为0。其主要原理是先通过离散傅里叶变换把图像转换到频域，再进行巴特沃斯低通滤波，然后用傅里叶逆变换转换回空域，最终实现图像增强效果。其中离散傅里叶变换就是傅里叶变换在时间和频率域上均以离散的形式存在。其一维离散傅里叶变换定义为：
在这里插入图片描述
其中 $f (x)$ 为离散序列，长度为 $N$ ， $u = 0, 1, \dots, N - 1$ ，则 $F (u)$ 的一维离散傅里叶反变换可以定义为：

在这里插入图片描述
根据公式（3）和（4）可以推广出二维离散傅里叶变换，设离散函数 $f (x, y)$ ，其中 $x = 0, 1, \dots, M - 1, y = 0, 1, \dots, N - 1$ 。
其定义为：

在公式（5）中， $u = 0, 1, \dots, M - 1 ， v = 0, 1, \dots, N - 1$ 。 $F (u, v)$ 二维离散傅里叶反变换可以定义为：
在这里插入图片描述
通过傅里叶变换将图像从时域转换到频域之后，就可以对其进行下一步操作，巴特沃斯低通滤波的定义如公式(7)所示：

在这里插入图片描述
其中 $D (u, v)$ 为该点到中心点的距离， $D_0$ 为截止频率，也就是振幅下降为-3 分贝时的频率， $D_0$ 选取的值不同其对图像的处理效果也不同。滤波器阶数 $n$ 越大滤波器的形状越陡峭。通常取 $n = 2$ ，此时没有明显的振铃效应，数值更大时会有模糊效应。
巴特沃斯高通滤波的定义如公式(8)所示：
在这里插入图片描述
为了在高频率域来降噪时较好的保存边缘信息，首要工作就是要找到图像结构特征的频率带，在这些频率带上不能做过多的缩减，而在离这些频率带比较远的区域，可以增加缩减的程度。巴特沃斯响应函数在频率中心区域的半径小于 $D (u, v)$ 时，函数值很接近1，变化不大；而随着距离中心点越来越远，函数值迅速减少。这一特性可以被用来进行图像降噪。

2、参考文献

[1]贾亮,邢轶博,徐善博. 基于改进巴特沃斯滤波的红外图像增强算法 [J]. 电脑与电信, 2022, (07): 58-62+76. DOI:10.15966/j.cnki.dnydx.2022.07.017.

二、代码实现

clc;
clear;
close all;

%% ------------------------------读取图像----------------------------------
Img = imread('Y.png');
%% ---------------------------可视化原始图像-------------------------------
subplot(121),imshow(Img);
title('原始图像');
%% --------------------------巴特沃斯高通滤波------------------------------
% 图像的数据类型由uint8转化为double类型
Imgd = im2double(Img);
% 二维傅立叶变换
fuv = fft2(Imgd);
% 对傅里叶变换后得到的频谱进行平移,将变换后的图像频谱中心从矩阵的原点移到矩阵的中心
trans = fftshift(fuv);
[row,col] = size(trans); % 求二维傅里叶变换后图像大小
n = 2;                   % 巴特沃斯高通滤波器的阶数n
d0 = 30;                 % 巴特沃斯高通滤波器的截止频率D0
n1 = round(row/2);
n2 = round(col/2);
% 计算频率点(i,j)与频域中心的距离
for i = 1:row      
    for j = 1:col 
        distance = sqrt((i-n1)^2+(j-n2)^2);
        if distance==0 
            h = 0; 
        else
            h = 1/(1+(d0/distance)^(2*n));% 巴特沃斯高通滤波
        end
        trans(i,j) = h*trans(i,j);% 频域图像乘以滤波器的系数
    end
end
invTrans = ifftshift(trans);
% 二维傅里叶反变换转换为时域图像
fxy = ifft2(invTrans);
% real函数取元素的实部
trans = real(fxy);
%% ----------------------可视化滤波后的图像--------------------------------
subplot(122),imshow(trans,[]);
title('巴特沃斯高通滤波图像');