Main Content

傅里叶变换

傅里叶变换的定义

傅里叶变换是将图像表示为不同幅值、频率和相位的复指数之和。傅里叶变换在广泛的图像处理应用中起着至关重要的作用,包括增强、分析、还原和压缩。

如果 f(m,n) 是两个离散空间变量 mn 的函数,则 f(m,n)二维傅里叶变换由如下关系定义:

F(ω1,ω2)=m=n=f(m,n)ejω1mejω2n.

变量 ω1ω2 是频率变量;其单位是弧度/采样。F(ω12) 通常称为 f(m,n)频域表示。F(ω12) 是复数值函数,在 ω1ω2 中均呈现周期性,期间为 2π。由于具有周期性,通常只显示范围 πω1,ω2π。请注意,F(0,0)f(m,n) 的所有值的总和。因此,F(0,0) 通常称为傅里叶变换的常量分量DC 分量。(DC 表示直流电;这是一个电气工程术语,指恒压电源,不同于电压呈正弦变化的电源。)

变换的逆运算是可应用于变换后的图像以生成原始图像的运算。二维傅里叶逆变换由下式给出:

f(m,n)=14π2ω1=ππω2=ππF(ω1,ω2)ejω1mejω2ndω1dω2.

粗略地说,此方程意味着 f(m,n) 可以表示为无限个不同频率的复指数(正弦曲线)的总和。F(ω12) 给出了 12) 频率的贡献的幅值和相位。

傅里叶变换的可视化

为说明傅里叶变换,我们以下面的函数 f(m,n) 为例,它在矩形区域内等于 1,在其他位置等于 0。为了简化图,f(m,n) 显示为连续函数,即使变量 mn 是离散的。

矩形函数

Plot of rectangular function f(m,n) in the spatial domain

下图以网格图形式显示上图所示的矩形函数的傅里叶变换的幅值,即

|F(ω1,ω2)|,

。幅值的网格图是可视化傅里叶变换的常用方法。

矩形函数的幅值图像

Mesh plot of the magnitude of the Fourier transform of the rectangular function f(m,n), plotted as a function of the horizontal frequency, ω1, and vertical frequency, ω2

图中心的峰值是 F(0,0),这是 f(m,n) 中所有值的总和。该图还显示,F(ω12) 在高水平频率下比在高垂直频率下具有更多能量。这反映出 f(m,n) 的水平横截面是窄脉冲,而垂直横截面是宽脉冲。窄脉冲比宽脉冲具有更多高频成分。

另一种可视化傅里叶变换的常见方法是将

log|F(ω1,ω2)|

显示为图像,如图所示。

矩形函数的傅里叶变换的对数

2-D pot of the log of the Fourier transform of the rectangular function f(m,n), plotted as a function of the horizontal frequency, ω1, and vertical frequency, ω2

F(ω12) 非常接近 0 的区域,使用对数有助于显示傅里叶变换的细节。

其他简单形状的傅里叶变换示例如下所示。

一些简单形状的傅里叶变换

2-D plots of the Fourier transforms of a rectangle tilted forty-five degrees, a circle centered at (0, 0), and an X shape centered at (0, 0)

离散傅里叶变换

在计算机上使用傅里叶变换通常涉及一种称为离散傅里叶变换 (DFT) 的变换形式。离散变换是一种其输入和输出值均为离散样本的变换,便于计算机操作。使用这种形式的变换有两个主要原因:

  • DFT 的输入和输出均为离散值,便于计算机操作。

  • 有一种计算 DFT 的快速算法,称为快速傅里叶变换 (FFT)。

DFT 通常针对仅在有限区域 0mM10nN1 上为非零值的离散函数 f(m,n) 定义。二维 M×N DFT 和逆 M×N DFT 的关系由下式给出:

F(p,q)=m=0M1n=0N1f(m,n)ej2πpm/Mej2πqn/N   p=0, 1, ..., M1q=0, 1, ..., N1

f(m,n)=1MNp=0M1q=0N1F(p,q)ej2πpm/Mej2πqn/N   m=0, 1, ..., M1 n=0, 1, ..., N1

F(p,q)f(m,n) 的 DFT 系数。零频率系数 F(0,0) 通常称为“DC 分量”。DC 是一个电气工程术语,表示直流电。(请注意,MATLAB® 中的矩阵索引始终从 1 开始,而不是从 0 开始;因此,矩阵元素 f(1,1)F(1,1) 分别对应于数学量 f(0,0)F(0,0)。)

MATLAB 函数 fftfft2fftn 分别实现用于计算一维 DFT、二维 DFT 和 N 维 DFT 的快速傅里叶变换算法。函数 ifftifft2ifftn 计算逆 DFT。

与傅里叶变换的关系

DFT 系数 F(p,q) 是傅里叶变换 F(ω12) 的采样。

F(p,q)=F(ω1,ω2)|ω1=2πp/Mω2=2πq/Np=0,1,...,M1q=0,1,...,N1

离散傅里叶变换的可视化

  1. 构造矩阵 f,它类似于傅里叶变换的定义的示例中的函数 f(m,n)。前面提到过,f(m,n) 在矩形区域内等于 1,在其他位置等于 0。使用二值图像来表示 f(m,n)

    f = zeros(30,30);
    f(5:24,13:17) = 1;
    imshow(f,"InitialMagnification","fit")

    Binary image representation of f(m,n)

  2. 使用以下命令计算并可视化 f 的 30×30 DFT。

    F = fft2(f);
    F2 = log(abs(F));
    imshow(F2,[-1 5],"InitialMagnification","fit");
    colormap(jet); colorbar

    在无填充情况下计算的离散傅里叶变换

    2-D plot of the 30-by-30 discrete Fourier transform of the binary rectangular function

    该图不同于傅里叶变换的可视化中显示的傅里叶变换。首先,傅里叶变换的采样要粗略得多。其次,零频系数显示在左上角,而不是传统的中心位置。

  3. 为了获得傅里叶变换的更精细采样,在计算 f 的 DFT 时可向在其中填零。使用以下命令,可一步执行填零和 DFT 计算。

    F = fft2(f,256,256);

    该命令在计算 DFT 之前用零将 f 填充为 256×256。

    imshow(log(abs(F)),[-1 5]); colormap(jet); colorbar

    在填充情况下计算的离散傅里叶变换

    2-D plot of the 30-by-30 discrete Fourier transform of the binary rectangular function with zero padding. The transform with zero padding has a finer frequency resolution.

  4. 然而,零频系数仍显示在左上角,而不是中心位置。您可以使用函数 fftshift 解决此问题,该函数交换 F 的象限,使零频系数位于中心位置。

    F = fft2(f,256,256);F2 = fftshift(F);
    imshow(log(abs(F2)),[-1 5]); colormap(jet); colorbar

    生成的图与傅里叶变换的可视化中所示的图相同。

傅里叶变换的应用

本节介绍傅里叶变换的许多图像处理相关应用中的一些应用。

线性滤波器的频率响应

线性滤波器的冲激响应的傅里叶变换可得到滤波器的频率响应。函数 freqz2 计算并显示滤波器的频率响应。高斯卷积核的频率响应表明,此滤波器允许低频通过,但会使高频衰减。

h = fspecial("gaussian");
freqz2(h)

高斯滤波器的频率响应

Mesh plot of the magnitude of the frequency response of a Gaussian filter.

有关线性滤波、滤波器设计和频率响应的详细信息,请参阅Design Linear Filters in the Frequency Domain

使用傅里叶变换执行快速卷积

此示例说明如何使用傅里叶变换对两个矩阵执行快速卷积。傅里叶变换的一个关键特性是两个傅里叶变换相乘对应于相关联的空间函数的卷积。此特性与快速傅里叶变换一起构成了快速卷积算法的基础。

注意:基于 FFT 的卷积方法最常用于大型输入。对于小型输入,使用 imfilter 函数通常更快。

创建两个简单的矩阵 ABA 是 M×N 矩阵,B 是 P×Q 矩阵。

A = magic(3);
B = ones(3);

AB 执行填零,使它们至少为 (M+P-1)×(N+Q-1)。(通常填零后的 AB 的大小要为 2 的幂,这时 fft2 的执行速度最快。)该示例将矩阵填充为 8×8。

A(8,8) = 0;
B(8,8) = 0;

使用 fft2 函数计算 AB 的二维 DFT。将两个 DFT 相乘,并使用 ifft2 函数计算结果的逆二维 DFT。

C = ifft2(fft2(A).*fft2(B));

提取结果的非零部分,并删除舍入误差导致的虚部。

C = C(1:5,1:5);
C = real(C)
C = 5×5

    8.0000    9.0000   15.0000    7.0000    6.0000
   11.0000   17.0000   30.0000   19.0000   13.0000
   15.0000   30.0000   45.0000   30.0000   15.0000
    7.0000   21.0000   30.0000   23.0000    9.0000
    4.0000   13.0000   15.0000   11.0000    2.0000

执行基于 FFT 的相关性分析来定位图像特征

此示例说明如何使用傅里叶变换来执行相关性分析,这与卷积密切相关。相关性可用于定位图像中的特征。在这种情况下,相关性通常称为模板匹配

将示例图像读入工作区。

bw = imread('text.png');

通过从图像中提取字母“a”来创建匹配模板。请注意,您也可以使用 imcrop 函数的交互式语法来创建模板。

a = bw(32:45,88:98);

通过将模板图像旋转 180 度,然后使用基于 FFT 的卷积技术,计算模板图像与原始图像的相关性。(如果将卷积核旋转 180 度,卷积等效于相关性。)要将模板与图像匹配,请使用 fft2ifft2 函数。在生成的图像中,明亮的峰值对应于字母的出现次数。

C = real(ifft2(fft2(bw) .* fft2(rot90(a,2),256,256)));
figure
imshow(C,[]) % Scale image to appropriate display range.

Figure contains an axes object. The axes object contains an object of type image.

要查看模板在图像中的位置,请找到最大像素值,然后定义一个小于该最大值的阈值。阈值化图像将这些峰值的位置显示为阈值化相关性图像中的白色斑点。(为了使位置在此图中更容易看到,此示例扩大了阈值化图像以放大点的大小。)

max(C(:))
ans = 68
thresh = 60; % Use a threshold that's a little less than max.
D = C > thresh;
se = strel('disk',5);
E = imdilate(D,se);
figure
imshow(E) % Display pixels with values over the threshold.

Figure contains an axes object. The axes object contains an object of type image.

另请参阅

| |

相关主题