频率域滤波（1）

一、频率域基础

频率域滤波实际上是将图像进行傅里叶变换，然后在变换域进行处理，然后进行傅里叶反变换转换回空间域，原理是用傅里叶变换表示的函数特征完全可以通过傅里叶反变换来重建，而且不会丢失任何信息（因为任何周期或非周期函数都可以表示为不同频率的正弦函数和余弦函数之和的形式）。实际上，空间域滤波和频率域滤波经常是对应的：

 　　

空间滤波实际上是图像与各种空间滤波器（模板）的卷积，而空间卷积的傅里叶变换是频率域中相应变换的乘积，因此频率域滤波可以用图像的傅里叶变换乘以相应的频率域滤波器。

 

连续变量函数的傅里叶变换暂且不提，我们直接来看图像的二维离散傅里叶变换，令f（x,y）表示一幅大小为[M,N]像素的数字图像，由F（u,v）表示图像的二维离散傅里叶变换（DFT）：

指数项由欧拉公式获得，可以将其展开为正弦函数和余弦函数，频率域是使用u和v作为（频率）变量，由F(u，v)构成坐标系。由u和v构成的大小为[M,N]的矩形区域称为频率矩形，大小和输入图像的大小相同。

离散傅里叶反变换(IDFT)的形式为:

频率域原点处变换的值（F(0,0)）称为傅里叶变换的直流分量，不难看出，F(0,0)等于f(x,y)平均值的M*N倍。还需要注意的是即使f(x,y)是实函数，它的变换通常也是复数。因此，傅里叶变换常用的是F(u,v)的频谱和相角，并将其显示为一幅图像。matlab有专有的快速傅里叶变换（FFT）算法实现，可以直接调用：F=fft2(f)。一般来说，使用傅里叶变换滤波时，需要对输入数据进行0填充：F=fft2(f,P,Q)，结果函数的大小为P*Q。令R(u,v)和I(u,v)表示其实部和虚部，则傅里叶频谱（幅度）为：

傅里叶谱可以实验函数abs直接获得：S=abs(F)。变换的相角定义为：

相角可以用atan2调用：phi=atan2(imag(F),real(F))；或者直接使用函数angle：phi=angle(F)。有些地方还会提到功率谱，其实就是频谱的平方。实函数的傅里叶变换是关于原点共轭对称的，因此频谱是关于原点偶对称的，而相角是关于原点奇对称的。由此得出的结论是傅里叶变换在u,v方向上是无穷周期的，而通过反变换得到的图像也是无穷周期的，我们感兴趣的是在区间[0,M-1]上的一个完整周期，变换之前，通过将f(x,y)乘以（-1）^(x+y)可以得到期望的周期，调用fftshift可以实现。最重要的是零频率项与f(x,y)的平均值成正比：

傅里叶变换后，频率矩阵的原点不在中心，可以使用函数fftshift将变换的原点移动到频率矩形的中心：Fc=fftshift(F)；与中心处亮度值占支配地位的8位显示相比，谱中的范围很大，显示不明显，通过对数变换可以解决这个问题：S2=log(1+abs(Fc))。具体的例子如下：

matlab代码如下：

%傅里叶变换
f=imread('G:\数字图像处理（冈萨雷斯）\DIP3E_CH04_Original_Images\DIP3E_Original_Images_CH04\Fig0424(a)(rectangle).tif');
subplot(231);imshow(f);title('原图');
F = fft2(f);           %傅里叶变换
subplot(232);imshow(F,[]);title('傅里叶变换');
S=abs(F);             %傅里叶谱
subplot(233);imshow(S,[]);title('傅里叶频谱')
Fc = fftshift(F);      %将变换的原点移动到频率矩形的中心
Sc=abs(Fc);           %居中的谱
subplot(234);imshow(Sc,[]);title('居中的谱')
S2 = log(1+abs(Sc));
subplot(235);imshow(S2,[]);title('对数变换后的谱')
f2=real(ifft2(F));          %傅里叶反变换
subplot(236);imshow(f2,[]);title('傅里叶反变换')

phi = atan2(imag(F),real(F));
%subplot(236);imshow(phi,[]);

　　频率域变换中的低频与图像中缓慢变化的灰度成分有关，高频由灰度的尖锐过渡造成，如边缘和噪声等，衰减高频而通过低频的滤波器（低通滤波器）将模糊一幅图像，相反的高通滤波器将增强尖锐的细节，但会降低图像的对比度。

二、有填充和无填充的滤波效果

空间滤波由滤波模板卷积一幅图像组成，函数相对位移，直到一个函数全部滑过另一个函数为止，在频率域中，F和H是周期的，表明在离散频率域中执行的卷积也是周期的，保证空间和频率域给出相同结果的唯一办法就是适当的0填充。假设函数f(x,y)和h(x,y)的大小均为M*N，通过对f和g补零，构造大小均为P*Q的填充函数：P>=2M-1,Q>=2N-1。填充函数为：

function PQ = paddedsize(AB,CD,PARAM)

if nargin ==1
    PQ = 2*AB;
elseif nargin ==2&&~ischar(CD)
    PQ = AB+CD-1;
    PQ = 2*ceil(PQ/2);
elseif nargin ==2
    m = max(AB);
    P=2^nextpow2(2*m);
    PQ=[P,P];
elseif (nargin==3)&&strcmpi(PARAM,'pwr2')
    m=max([AB CD]);
    P=2^nextpow2(2*m);
    PQ=[P,P];
else
    error('wrong number of imputs');
end

　　具体是否填充的例子如下：

当滤波器位于虚线图像的一侧时，他将会遇到邻近这一侧的周期分量中的一个相同区域，因为一个恒定区域的均值是同一常数，因此结果的这一部分没有被模糊。代码如下：

f=imread('G:\数字图像处理（冈萨雷斯）\DIP3E_CH04_Original_Images\DIP3E_Original_Images_CH04\Fig0432(a)(square_original).tif');
subplot(221);imshow(f);title('原图');
[M,N]=size(f);
%f = im2double(f);  %转换为浮点类型再计算
[f,revertclass] = tofloat(f);
F = fft2(f);             %傅里叶变换
H = lpfilter('gaussian',M, N,10);      %频域高斯低通滤波器
G=H.*F;                   %滤波（乘积）
g=ifft2(G);          %傅里叶反变换
%g=im2uint8(g);
g=revertclass(g);
subplot(222);imshow(g);title('无填充时频率域低通滤波');

%经过填充后的傅里叶变换及其滤波
PQ = paddedsize(size(f));    %以0填充，
Fp = fft2(f,PQ(1),PQ(2));    %用填充的方式计算傅里叶变换
Hp = lpfilter('gaussian',PQ(1),PQ(2),2*10);   %频域高斯低通滤波器
Gp = Hp.*Fp;
gp = ifft2(Gp);
gpc = gp(1:size(f,1),1:size(f,2)); %只取右上角一个周期的图像
gpc = im2uint8(gpc);
subplot(223);imshow(gpc);title('有填充时频域低通滤波');
subplot(224);imshow(gp);title('未裁剪的全填充图像');