数字图像处理实验（5）：PROJECT 04-01 [Multiple Uses]，Two-Dimensional Fast Fourier Transform 标签： 图像处理MATLAB数字图像处理

实验要求：

Objective:
To further understand the well-known algorithm Fast Fourier Transform (FFT) and verify its effectiveness to calculate the discrete Fourier transform (DFT).
Main requirements:
Ability of programming with C, C++, or Matlab.
Instruction manual:
The purpose of this project is to develop a 2-D FFT program “package” that will be used in several other projects that follow. Your implementation must have the capabilities to:
(a) Multiply the input image by (-1)x+y to center the transform for filtering.
(b) Multiply the resulting (complex) array by a real function (in the sense that the the real coefficients multiply both the real and imaginary parts of the transforms). Recall that multiplication of two images is done on pairs of corresponding elements.
(c) Compute the inverse Fourier transform.
(d) Multiply the result by (-1)x+y and take the real part.
(e) Compute the spectrum.
Basically, this project implements Fig. 4.5. If you are using MATLAB, then your Fourier transform program will not be limited to images whose size are integer powers of 2. If you are implementing the program yourself, then the FFT routine you are using may be limited to integer powers of 2. In this case, you may need to zoom or shrink an image to the proper size by using the program you developed in Project 02-04.
An approximation: To simplify this and the following projects (with the exception of Project 04-05), you may ignore image padding (Section 4.6.3). Although your results will not be strictly correct, significant simplifications will be gained not only in image sizes, but also in the need for cropping the final result. The principles will not be affected by this approximation.

实验要求我们编写程序实现2维图像的快速傅里叶变换，并运用于频率域滤波。
简要介绍一下频率域滤波的步骤：

1. 给定一幅大小为M*N的输入图像f(x,y)，计算得到填充参数P和Q。通常，取P=2M，Q=2N；
2. 添加必要数量的0，形成大小为P*Q的填充后的图像fp(x,y)；
3. 用(-1)^(x+y)乘以fp(x,y)，将图像移到其傅里叶变换的中心处；
4. 计算（3）中的图像DFT，得到F(u,v)；
5. 生成一个实的、对称的滤波器函数H(u,v)，其大小与fp(x,y)相同，为P*Q，中心在(P/2, Q/2)。用阵列相乘形成乘积G(u,v)；
6. 计算（5）中结果的逆DFT，取出实部，得到gp(x,y)；
7. 从gp(x,y)的左上象限提取M*N的图像，得到最终处理结果g(x,y)。

上程序：

%%
close all;
clc;
clear all;

%%
% 读取图像
img = imread('gray_image.jpg');
figure(1)
imshow(img);
title('original A');

% 得到填充参数P和Q
[M, N] = size(img);
P = 2 * M;
Q = 2 * N;

img = double(img);
% 添加必要数量的0
img_fp = zeros(P, Q);
img_fp(1:M, 1:N) = img(1:M, 1:N);

figure(2);
imshow(img_fp, []);
title('image B');

% 用(-1)^(x+y)乘以图像的结果
for x = 1:P
    for y = 1:Q
        img_fp(x, y) = img_fp(x, y) .* (-1)^(x+y);
    end
end

% figure(3);
% imshow(img_fp, []);
% title('image C');

% 对图像做快速傅里叶变换
img_Fp = fft2(img_fp);

% figure(4);
% imshow(img_Fp, []);
% title('image D');

% H = ones(P, Q);
% H(P/2, Q/2) = 0;

% H = zeros(P, Q);
% H(P/2, Q/2) = 1;
% H(P/2 - 1, Q/2 - 1) = 1;

r = 30;
H = ones(P, Q);
for x = 1:P
    for y = 1:Q
        d = sqrt((x-M)^2 + (y-N)^2);
        if d > r
            H(x, y) = 0;
        else
            H(x, y) = 1;
        end

%         if x == P/2 && y == Q/2
%             H(x, y) = 1;
%         else
%             H(x, y) = 1;
%         end
    end
end

figure(5);
imshow(H, []);
title('image E');

img_G = img_Fp .* H;

% figure(6);
% imshow(img_G, []);
% title('image F');

img_g = ifft2(img_G);
img_g = real(img_g);

for x = 1:P
    for y = 1:Q
        img_g(x, y) = img_g(x, y) .* (-1)^(x+y);
    end
end

% figure(7);
% imshow(img_g, []);
% title('image G');

img_o = img_g(1:M, 1:N);

figure(8);
imshow(img_o, []);
title('result H');

imwrite(img_o, 'result.jpg');

实验结果：

上图是原图像；

上图是步骤1和2，扩展后的图像fp(x,y)；

上图是滤波器函数H(u,v)，中心一个圆内为1，圆外都是0，表示一个低通滤波器；

上图是逆DFT计算得到的gp(x,y)；

上图为最终的结果，从gp(x,y)提取出的g(x,y)图像；
明显得，可以发现，低通滤波器使原图像变得模糊了。