Gabor函式
　　Gabor變換屬於加窗傅立葉變換，Gabor函式可以在頻域不同尺度、不同方向上提取相關的特徵。另外Gabor函式與人眼的生物作用相仿，所以經常用作紋理識別上，並取得了較好的效果。二維Gabor函式可以表示為：
　　其中：
　　v的取值決定了Gabor濾波的波長，u的取值表示Gabor核函式的方向，K表示總的方向數。引數決定了高斯視窗的大小，這裡取。程式中取4個頻率（v=0, 1, ..., 3），8個方向（即K=8，u＝0, 1, ... ,7），共32個Gabor核函式。不同頻率不同方向的Gabor函式可通過下圖表示：
　　圖片來源：GaborFilter.html
　　圖片來源：http://www.bmva.ac.uk/bmvc/1997/papers/033/node2.html
　　三、程式碼實現
　　Gabor函式是復值函式，因此在運算過程中要分別計算其實部和虛部。程式碼如下：

　　private void CalculateKernel(int Orientation, int Frequency)
　　{
　　double real, img;
　　for(int x = -(GaborWidth-1)/2; x<(GaborWidth-1)/2+1; x++)
　　for(int y = -(GaborHeight-1)/2; y<(GaborHeight-1)/2+1; y++)
　　{
　　real = KernelRealPart(x, y, Orientation, Frequency);
　　img = KernelImgPart(x, y, Orientation, Frequency);
　　KernelFFT2[(x+(GaborWidth-1)/2) + 256 * (y+(GaborHeight-1)/2)].Re = real;
　　KernelFFT2[(x+(GaborWidth-1)/2) + 256 * (y+(GaborHeight-1)/2)].Im = img;
　　}
　　}
　　private double KernelRealPart(int x, int y, int Orientation, int Frequency)
　　{
　　double U, V;
　　double Sigma, Kv, Qu;
　　double tmp1, tmp2;
　　U = Orientation;
　　V = Frequency;
　　Sigma = 2 * Math.PI * Math.PI;
　　Kv = Math.PI * Math.Exp((-(V+2)/2)*Math.Log(2, Math.E));
　　Qu = U * Math.PI / 8;
　　tmp1 = Math.Exp(-(Kv * Kv * ( x*x + y*y)/(2 * Sigma)));
　　tmp2 = Math.Cos(Kv * Math.Cos(Qu) * x + Kv * Math.Sin(Qu) * y) - Math.Exp(-(Sigma/2));
　　return tmp1 * tmp2 * Kv * Kv / Sigma;
　　}
　　private double KernelImgPart(int x, int y, int Orientation, int Frequency)
　　{
　　double U, V;
　　double Sigma, Kv, Qu;
　　double tmp1, tmp2;
　　U = Orientation;
　　V = Frequency;
　　Sigma = 2 * Math.PI * Math.PI;
　　Kv = Math.PI * Math.Exp((-(V+2)/2)*Math.Log(2, Math.E));
　　Qu = U * Math.PI / 8;
　　tmp1 = Math.Exp(-(Kv * Kv * ( x*x + y*y)/(2 * Sigma)));
　　tmp2 = Math.Sin(Kv * Math.Cos(Qu) * x + Kv * Math.Sin(Qu) * y) - Math.Exp(-(Sigma/2));
　　return tmp1 * tmp2 * Kv * Kv / Sigma;
　　}
 

有了Gabor核函式後就可以採用前文中提到的“離散二維疊加和卷積”或“快速傅立葉變換卷積”的方法求解Gabor變換，並對變換結果求均值和方差作為提取的特徵。32個Gabor核函式對應32次變換可以提取64個特徵（包括均值和方差）。由於整個變換過程程式碼比較複雜，這裡僅提供測試程式碼供下載。該程式碼僅計算了一個101×101尺寸的Gabor函式變換，得到均值和方差。程式碼採用兩種卷積計算方式，從結果中可以看出，快速傅立葉變換卷積的效率是離散二維疊加和卷積的近50倍。
-----------------------------------------------------------------------

最近忙著論文，需要Gabor濾波程式碼，可是網上總找不到合適的程式碼，於是就自己編了一個，不當之處請指點。參考論文為 L. Wiskott，J. M. Fellous，N. Kruger，C. v. d.Malsburg. Face Recognition by Elastic Bunch Graph Matching，IEEE Trans. On PAMI，Vol.19，No.7，pp775-779，1997

首先實現濾波器：

function [bank] = do_createfilterbank(imsize,varargin)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 函式實現：建立Gabor 濾波組
%
%%% 必選引數：
%   imsize - 影象大小
%%% 可選引數：
%   freqnum — 頻率數目
%   orientnum — 方向數目
%   f       —   頻率域中的取樣步長
%   kmax    —   最大的取樣頻率
%   sigma   —   高斯窗的寬度與波向量長度的比率
%
%%% 返回結果：
%   bank
%           .freq        —     濾波頻率
%           .orient      —     濾波方向
%           .filter      —     Gabor濾波
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

conf = struct(...,
    'freqnum',3,...
    'orientnum',6,...
    'f',sqrt(2),...
    'kmax',(pi/2),...
    'sigma',(sqrt(2)*pi) ...
    );

conf = do_getargm(conf,varargin);

bank = cell(1,conf.freqnum*conf.orientnum);
for f0=1:conf.freqnum
    fprintf('處理頻率 %d \n', f0);
    for o0=1:conf.orientnum
        [filter_,freq_,orient_] = do_gabor(imsize,(f0-1),(o0-1),conf.kmax,conf.f,conf.sigma,conf.orientnum);
        bank{(f0-1)*conf.orientnum + o0}.freq = freq_; %以orient增序排列
        bank{(f0-1)*conf.orientnum + o0}.filter = filter_;
        bank{(f0-1)*conf.orientnum + o0}.orient = orient_;
    end
end

for ind = 1:length(bank)
    bank{ind}.filter=fftshift(bank{ind}.filter);
end

function [filter,Kv,Phiu] = do_gabor (imsize,nu,mu,Kmax,f,sigma,orientnum)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 函式實現： 建立Gabor濾波
%
%%% 引數：
%   imsize     : 濾波的大小（即影象大小）
%   nu         : 頻率編號 [0 ...freqnum-1];
%   mu         : 方向編號 [0...orientnum-1]
%   Kmax       ： 最大的取樣頻率
%   f          ： 頻率域中的取樣步長
%   sigma      : 高斯窗的寬度與波向量長度的比率
%   orientnum : 方向總數
%
%%% 返回值：
%   filter :    濾波
%   Kv    :    頻率大小
%   Phiu :    方向大小
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rows = imsize(1);
cols = imsize(2);
minrow = fix(-rows/2);
mincol = fix(-cols/2);
row = minrow + (0:rows-1);
col = mincol + (0:cols-1);
[X,Y] = meshgrid(col,row);

Kv = Kmax/f^nu;
Phiu = pi * mu /orientnum;
K = Kv * exp(i * Phiu);

F1 = (Kv ^ 2)/ (sigma^2) * exp(-Kv^2 * abs(X.^2 + Y.^2) / (2*sigma^2)) ;
F2 = exp(i * (real(K) * X + imag(K) * Y)) - exp(-sigma^2/2);
filter = F1.* F2;

Gabor 濾波實現(1)已經建立了Gabor濾波組，現在可以使用該濾波組對影象進行轉換，得到振幅和相位。

function [result] = do_filterwithbank(im,bank)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 函式實現：對影象使用Gabor濾波組進行轉換換
%
%%% 引數：
%   im       — 被轉換的影象
%   bank        — 由函式do_createfilterbank得到的濾波組
%
%%% 返回：
%   result                     — 影象被轉換後的結果
%       .amplitudes            — 不同畫素點的振幅向量
%       .phases                — 不同畫素點的相位向量
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[N1 N2] = size(im);
N3 = length(bank);
phases = zeros(N1,N2,N3);
amplitudes = zeros(N1,N2,N3);
imagefft = fft2(im);

for ind = 1:N3
    fprintf('正在處理濾波 %d \n',ind);
     temp = ifft2(imagefft .* bank{ind}.filter);
     phases(:,:,ind) = angle(temp);
     amplitudes(:,:,ind) = abs(temp);
end
result.phases = phases;
result.amplitudes = amplitudes;

整個程式可以如下使用。 im = imread('image.jpg'); im = rgb2gray(im); bank = do_createfilterbank(size(im)); result = do_filterwithbank(im,bank);