傅立叶变换与频率域滤波

实验四傅立叶变换与频率域滤波

实验目的

通过本次实验，实现以下几个目标：

1.理解傅立叶变换；

2.熟悉MATLAB中各种傅立叶变换相关的函数；

3.掌握频域滤波的步骤以及MATLAB的实现方法；

4.理解频域滤波器与空域滤波器的关系。

实验内容

一、傅立叶变换及傅立叶反变换

1.傅立叶变换相关函数

MATLAB提供了几个和傅里叶变换相关的函数。其说明如下：

F=fft2(f); 二维傅立叶变换

abs(F); 获得傅立叶频谱

fftshift(F); 将变换的原点移至频率矩形的中心

ifft2(F); 二维傅立叶反变换

real(ifft2(F)); 提取变换后的实部

imag(ifft2(F)); 提取变换后的虚部

2.傅里叶频谱

傅里叶频谱反映了图像的频率成分。下面的例子对课本中123页和

125页的图Fig4.03(a) 和图Fig4.04(a)进行傅立叶变换，得到傅立叶

频谱。显示傅立叶频谱时，使用了对数变换以获得更好效果。

f=imread('Fig4.03(a).jpg');

F=fft2(double(f));

F=fftshift(F);

figure(1), imshow(f);

figure(2), imshow(log(abs(F)+1),[ ]);

f=imread('Fig4.04(a).jpg');

F=fft2(double(f));

F=fftshift(F);

figure(1), imshow(f);

figure(2), imshow(log(abs(F)+1),[ ]);

为了更好地理解频谱，显示下面三个图像（x6.jpg，x60.jpg，y6.jpg）的傅里叶频谱，观察并比较、分析结果。

显示频谱时使用下面的语句来做灰度变换可找出其主要的频率成分。

figure(2), imshow(log(abs(F)+1).^4,[ ]); %先对数、再幂次变换

这些都是用Matlab 制作的空间域存在明显周期的图像，生成图像x60.jpg的代码如下。

for i=1:300

for j=1:300

if(mod(j-1,60)<30) c(i,j)=255;

else c(i,j)=0;

end

end

end

imwrite(c,'x60.jpg')

3．傅立叶变换对

将一幅图像进行傅立叶变换，再进行傅立叶反变换，可以得到原始图像。下面的例子对课本中134页的图Fig4.11(a)进行傅立叶变换，然后再进行傅立叶反变换，观察并了解实现过程。

f=imread('Fig4.11(a).jpg');

F=fft2(f);

ff=ifft2(F);

fa= abs (real (ff));

figure(1),imshow(f);

figure(2),imshow(F);

figure(3),imshow(ff,[0 255]);

figure(4),imshow(fa,[0 255]);

理论上，经过傅里叶变换和傅里叶反变换，应该得到原始图像。

比较fa和f，基本差不多。但用下面的语句可以看到他们之间有细微的差别，为什么？

imshow(fa-double(f),[]);

下面的例子将一幅图像进行傅立叶变换和傅立叶反变换，并将0频率分量移到矩阵的中心，并将四幅图放在一起，便于观察。

f=imread('Fig4.11(a).jpg');

F=fft2(f);

F=fftshift(F);

ff=ifft2(F);

fa= abs (real (ff));

subplot(2,2,1);

imshow(f);

subplot(2,2,2);

imshow(log(abs(F)+1),[]);

subplot(2,2,3);

imshow(ff,[0 255]);

subplot(2,2,4);

imshow(fa,[0 255]);

二、频率域滤波

频域滤波的步骤在课本中124页有详细描述。按照该步骤，在MATLAB中很容易编程实现频域滤波。由于滤波器就是频率域中的函数，关键是如何构造滤波器函数。

1. 低通滤波

低通滤波是使低频部分通过，而使高频部分受到抑制，从而使图像变得平滑。常用的低通滤波有理想低通滤波、巴特沃斯低通滤波和高斯低通滤波。

下面的例子实现了对课本135页图Fig4.11(a)进行理想低通滤波处理，改变截止频率D0重复实验，对结果进行分析和比较。

f=imread('Fig4.11(a).jpg');

% 傅里叶变换，移位

F=fft2(double(f));

F=fftshift(F);

% 构造理想低通滤波器，并用它滤波

[height width]=size(F);

H(1: height,1: width)=0;

x0= height /2; y0= width /2;

for x=1:height

for y=1:width

if(sqrt((x- x0)*(x- x0)+(y-y0)*(y- y0))<30) % 注1

H(x,y)=1;

end

FF(x,y)=F(x,y)*H(x,y);

end

end

% 傅里叶反变换

g=ifft2(FF);

% 显示并比较结果

figure(1), imshow(f);

figure(2), imshow(log(abs(F)+1),[ ]);

figure(3), imshow(log(abs(FF)+1),[ ]);

figure(4),imshow(abs(real(g)),[ ]);

2. 高通滤波

高通滤波是使高频部分通过，而使低频部分受到抑制，从而使图像变得平滑。常用的高通滤波有理想高通滤波、巴特沃斯高通滤波和高斯高通滤波。

下面的例子实现了对课本145页图Fig4.11(a)进行理想高通滤波处理。改变截止频率D0重复实验，对结果进行分析和比较。

f=imread('Fig4.11(a).jpg');

% 傅里叶变换，移位

F=fft2(double(f));

F=fftshift(F);

% 理想高通滤波器，滤波

[height width]=size(F);

H(1: height,1: width)=1;

x0= height /2; y0= width /2;

for x=1:height

for y=1:width

if(sqrt((x- x0)*(x- x0)+(y-y0)*(y- y0))<30) % 注1

H(x,y)=0;

end

FF(x,y)=F(x,y)*H(x,y);

end

end

% 傅里叶反变换

g=ifft2(FF);

% 显示并比较结果

figure(1), imshow(f);

figure(2), imshow(log(abs(F)+1),[ ]);

figure(3), imshow(log(abs(FF)+1),[ ]);

figure(4),imshow(abs(real(g)),[]);

注1：将语句改为(x- x0)*(x- x0)+(y-y0)*(y- y0))<900免去了复杂的开平方运算sqrt，可大大减少计算量。

三、频域滤波器的空间形式

由卷积定理我们知道，空间域滤波和频域滤波之间存在对应关系。频域滤波器在空间域存在对应的空间域滤波器。

下面的例子得到了课本136页理想低通滤波器的空间形式。

F=imread('Fig4.13(a).jpg');

f=ifft2(double(F));

f=fftshift(f);

figure(1), imshow(F);

figure(2),imshow(abs(real(f)),[]);

下面的例子得到了课本137页巴特沃斯低通滤波器的空间形式。

F=imread('Fig4.14(b).jpg');

f=ifft2(double(F));

f=fftshift(f);

figure(1), imshow(F);

figure(2),imshow(abs(real(f)),[]);

下面的例子得到了课本143页理想高通滤波器的空间形式。

F=imread('Fig4.22(b).jpg');

f=ifft2(double(F));

f=fftshift(f);

figure(1), imshow(F);

figure(2),imshow(abs(real(f)),[]);

注：为了能够更好地观察空间形式的图像，可以使用一些灰度变换的技巧

附录

习题4.9参考代码

f=imread('FigProb4.09(left).jpg');

% 傅里叶变换，移位

F=fft2(double(f));

F=fftshift(F);

% 取复共轭

[height width]=size(F);

for x=1:height

for y=1:width

FF(x,y)= complex(real(F(x,y)),-imag(F(x,y)));

end

end

% 傅里叶反变换

g=ifft2(FF);

% 显示并比较结果

figure(1),imshow(f);

figure(2),imshow(abs(real(g)),[]);

高斯低通滤波器

f=imread('Fig4.11(a).jpg');

% 傅里叶变换，移位

F=fft2(double(f));

F=fftshift(F);

% 高斯低通滤波器，滤波

[height width]=size(F);

x0= height /2; y0= width /2;

twoD02=2*15*15;

H(1: height,1: width)=0;

for x=1:height

for y=1:width

D2=(x- x0)*(x- x0)+(y-y0)*(y- y0);

H(x,y)=exp(-D2/twoD02);

FF(x,y)=F(x,y)*H(x,y);

end

end

% 傅里叶反变换

g=ifft2(FF);

% 显示并比较结果

figure(1), imshow(f);

figure(2), imshow(log(abs(F)+1),[ ]);

figure(3), imshow(log(abs(FF)+1),[ ]);

figure(4),imshow(abs(real(g)),[]);