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 ，那么理想低通频域滤波器可形式化地描述为：

$LPF(i,j)=\left\{\begin{matrix} 1, & \sqrt{(i-row/2)^{2}+(j-col/2)^{2}} \leq D_{0}\\ 0, & \sqrt{(i-row/2)^{2}+(j-col/2)^{2}} > D_{0} \end{matrix}\right.$

其中 $D_{0}$ 表示理想低通滤波器的截止频率，滤波器的频率域原点在频谱图像的中心处，在以截止频率为半径的圆形区域之内的滤镜元素值全部为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('高通滤波')