数字图像的频域运算

实验六 数字图像的频域运算

一、实验目的

1. 掌握图像傅里叶变换的概念及其计算方法；
2. 熟练编程实现图像的傅里叶变换；
3. 编程实现图像的频域滤波；

二、实验环境

1. PC计算机
2. MatLab软件/语言包括图像处理工具箱(Image Processing Toolbox)
3. 实验所需要的图片

• 实验原理

自己写图像傅里叶变换的原理和流程！以及频域滤波的步骤！

傅里叶变换的原理：

傅里叶变换认为一个周期函数(信号)包含多个频率分量，任意函数（信号）f(t)可通过多个周期函数（基函数）相加而合成。从物理角度理解傅里叶变换是以一组特殊的函数（三角函数）为正交基，对原函数进行线性变换，物理意义便是原函数在各组基函数的投影。傅里叶变换是一种函数在空间域和频率的变换，从空间域到频率域的变换是傅里叶变换，而从频率域到空间域的转换叫做傅里叶的反变换。

傅里叶变换是一种时域与频域的转换关系。

傅里叶变换的流程：首先需要构建一个正弦函数f=sin（2*pi*f*t）不难看出本例中的正弦函数的频率为1，其中fs是我们的采样频率，当采样频率为100的时候 我们的采样周期1/fs 为0.01s 也就是每隔0.01s取一个点。然后恢复幅度轴，创建频率轴

频域滤波的步骤：

1、给定一幅大小为M*N的输入图像f(x,y), 由P >= 2M-1,Q >= 2N-1;得到填充参数P,Q，典型的我们选择，P=2M，Q=2N；

如果滤波的目的仅是粗糙的视觉分析，可以跳过此步骤；

2、对f(x,y) 添加必要数量的0，形成大小为P*Q的图像 fp(x,y);

3、用 乘以fp(x,y)移到其变换中心；

4、计算来自步骤3图像的DFT，得到F(u,v)；

5、生成一个实的、对称的滤波函数H(u,v)，其大小为P*Q，中心在（P/2，Q/2）处。用阵列相乘形成乘积G(u,v)=H(u,v)F(u,v)；

1. 得到处理后的图像  其中为忽略由于计算不准确导致的寄生复分量，选择了实部，下标p表示我们处理的是填充后的阵列；

还有一种处理方法是取模值；magnitude()函数

7、通过从的左上限提取M*N，得到最终的处理结果。

• 实验步骤和结果

1. 对图像进行傅里叶变换，并显示其频谱图

1. 对比中心化和不中心化的频谱图，频谱图中心化采用两种思路：

1. 用原始图像乘以（-1）x+y，就是课堂上讲解的步骤，进行平移

a代码如下：

i=imread('E:\ins风\yun.tif');

subplot(1,2,1);imshow(i);title('20101110227原始图');

size_i=size(i);

for a=1:size_i(1,1)

for b=1:size_i(1,2)

i(a,b)=i(a,b)*(-1)^(a+b);

end

end

ii=fft2(i);

for a=1:size_i(1,1)

for b=1:size_i(1,2)

im(a,b)=sqrt((real(ii(a,b)))^2+(imag(ii(a,b)))^2)/256^2;

end

end

subplot(1,2,2);imshow(im);title('用原始图像乘以 （-1)x+y');

1. 用fftshift函数

b代码如下：

img=imread('E:\ins风\yun.tif');

subplot(2,2,1);imshow(img);title('20101110227原图');

f=rgb2gray(img);    %对于RGB图像必须做的一步，也可以用im2double函数

F=fft2(f);          %傅里叶变换

F1=log(abs(F)+1);   %取模并进行缩放

subplot(2,2,2);imshow(F1,[]);title('傅里叶变换频谱图');

Fs=fftshift(F);      %将频谱图中零频率成分移动至频谱图中心

S=log(abs(Fs)+1);    %取模并进行缩放

subplot(2,2,3);imshow(S,[]);title('频移后的频谱图');

对比使用和不使用log变换的频谱图

代码如下：

clear

clc

img=imread('E:\ins风\yun.tif');

subplot(2,2,1);

imshow(img);

title('20101110227原图');

f=rgb2gray(img);    %对于RGB图像必须做的一步，也可以用im2double函数

F=fft2(f);          %傅里叶变换

subplot(2,2,2);imshow(F,[]);

title('傅里叶变换频谱图');

F2 = fftshift(log(1+abs(F)));% log变化

 subplot(2, 2, 3);

 imshow(F2, []);

 title('log变化频谱图');

Fs=fftshift(F);      %将频谱图中零频率成分移动至频谱图中心

subplot(2,2,4);

imshow(Fs,[]);

title('unlog频谱图');

生成图像如下：

1. 对图像进行用高斯低通滤波器进行滤波

1.  对比填充和不填充的滤波效果，填充的思路有三种：（a必须完成，b、c任选其一即可）

1. 按照课堂讲解，取填充大小为2M,2N，对原始的时域图像进行填充，使得填充后的图像左上角为原始图像，其余部分全部为0.
2. 取填充大小为2M,2N，用FFT2函数的变量，在进行傅里叶变换的时候，完成填充。
3. 填充的尺寸用paddedsize函数得到（这是课本上最基本的理论），再用FFT2函数的变量，在进行傅里叶变换的时候，完成填充。

1.  对比移中和不移中的滤波效果，进行分析。注意：移中时，要保证图像和滤波器一致。

代码如下：

    f = imread('E:\ins风\yun.tif'); %读入原图像

%f = mat2gray(f,[0 255]);

D0=20;

subplot(3,3,1),imshow(f),title('20101110227原图');

s=fftshift(fft2(f));

subplot(3,3,2),imshow(log(abs(s)),[]),title('原图频谱');

% 1.给定一幅大小为M*N的输入图像f(x,y),得到填充参数P = 2M,Q = 2N

[m,n] = size(f);

P = 2 * m;

Q = 2 * n;

% 2.对f(x,y)添加必要数量的0，形成大小为P * Q的填充后的图像fp(x,y)

fp = zeros(P,Q);

fp(1:m,1:n) = f(1:m,1:n);

ofp = zeros(P,Q);

ofp(1:m,1:n) = f(1:m,1:n);

% 3.用(-1)^(x+y)乘以fp(x,y)移到其交换的中心

for i = 1 : m

    for j = 1 : n

        fp(i,j) = double(fp(i,j)*(-1)^(i+j));

    end

end

% 4.计算来自步骤3的图像的DFT，得到F(u,v)

F = fft2(fp,P,Q);

oF = fft2(ofp,P,Q);

% 5.生成一个实的、对称的滤波函数H(u,v),其大小为P*Q，中心在（P/2,Q/2）处。

% 用阵列相乘形成乘积G(u,v) = H(u,v)F(u,v);即G(i,k)=H(i,k)F(i,k)

H = zeros(P,Q);

a = 2 * D0^2;

for u = 1 :P

    for v = 1:Q

        D = (u-(m+1.0))^2+(v-(n+1.0))^2;

        H(u,v) = exp((-D)/a);

    end

end

G = H.*F;

oH = zeros(P,Q);

oa = 2 * D0^2;

for ou = 1 :P

    for ov = 1:Q

        oD = ou^2+ov^2;

        oH(ou,ov) = exp((-oD)/oa);

    end

end

oG = oH.*oF;

% 6.得到处理后的图像

gp = ifft2(G);

gp = real(gp);

ogp = ifft2(oG);

ogp = real(ogp);

for i = 1 : m

    for j = 1 : n

        gp(i,j) = double(gp(i,j)*(-1)^(i+j));

    end

end

% 7.通过从gp(x,y)的左上象限提取M*N区域，得到最终处理结果g(x,y)

g(1:m,1:n) = gp(1:m,1:n);

og(1:m,1:n) = ogp(1:m,1:n);

subplot(3,3,3),imshow(g,[]),title('填充后的高斯低通滤波');

subplot(3,3,7),imshow(og,[]),title('不移中填充后的高斯低通滤波');

s1=fftshift(fft2(g));

subplot(3,3,4),imshow(log(abs(s1)),[]),title('填充后的高斯低通滤波频谱');

image_fft = fft2(f);%用傅里叶变换将图象从空间域转换为频率域

image_fftshift = fftshift(image_fft);

%将零频率成分(坐标原点)变换到傅里叶频谱图中心

[width,high] = size(f);

D = zeros(width,high);

%创建一个width行，high列数组，用于保存各像素点到傅里叶变换中心的距离

for i=1:width

    for j=1:high

        D(i,j) = sqrt((i-width/2)^2+(j-high/2)^2);

        %像素点(i,j)到傅里叶变换中心的距离

        H = exp(-1/2*(D(i,j).^2)/(D0*D0));

        %高斯低通滤波函数

        image_fftshift(i,j)= H*image_fftshift(i,j);

        %将滤波器处理后的像素点保存到对应矩阵

    end

end

image_result = ifftshift(image_fftshift);%将原点反变换回原始位置

image_result = uint8(real(ifft2(image_result)));

subplot(3,3,5),imshow(image_result),title('无填充的高斯低通滤波');

s2=fftshift(fft2(image_result));

subplot(3,3,6),imshow(log(abs(s2)),[]),title('无填充的高斯低通滤波频谱');

结果如下：

六、实验思考

我体会到了傅里叶变换的物理意义：傅里叶变换的物理意义#

1、图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。（灰度变化得快频率就高，灰度变化得慢频率就低）

傅立叶变换在实际中有非常明显的物理意义，设f是一个能量有限的模拟信号，则其傅立叶变换就表示f的谱。从纯粹的数学意义上看，傅立叶变换是将一个函数转换为一系列周期函数来处理的。从物理效果看，傅立叶变换是将图像从空间域转换到频率域，其逆变换是将图像从频率域转换到空间域。换句话说，傅立叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，傅立叶逆变换是将图像的频率分布函数变换为灰度分布函数

2、傅立叶变换以前，图像（未压缩的位图）是由对在连续空间（现实空间）上的采样得到一系列点的集合，我们习惯用一个二维矩阵表示空间上各点，则图像可由z=f(x,y)来表示。由于空间是三维的，图像是二维的，因此空间中物体在另一个维度上的关系就由梯度来表示，这样我们可以通过观察图像得知物体在三维空间中的对应关系。为什么要提梯度？因为实际上对图像进行二维傅立叶变换得到频谱图，就是图像梯度的分布图，当然频谱图上的各点与图像上各点并不存在一一对应的关系，即使在不移频的情况下也是没有。傅立叶频谱图上我们看到的明暗不一的亮点，实际上图像上某一点与邻域点灰度值差异的强弱，即梯度的大小，也即该点的频率的大小（差异/梯度越大，频率越高，在频谱图上就越亮（亮点多）。差异/梯度越小，频率越低，在频谱图上就越暗（暗点多）。一般来讲，梯度大则该点的亮度强，否则该点亮度弱。