数字图像处理与Python实现-图像变换-傅立叶变换

数字图像处理中，傅立叶变换是一种将图像从空间域转换到频域的重要工具。它通过将图像分解为不同频率的正弦和余弦波，揭示图像的频率特性，便于进行图像增强、去噪、压缩等操作。快速傅立叶变换（FFT）通过减少计算量，显著提高了傅立叶变换的效率。二维傅立叶变换将图像扩展到两个维度，广泛应用于图像频域滤波，能够增强或抑制特定频率的成分，从而实现图像增强或去噪。Python代码展示了如何利用FFT对图像进行变换，并通过频谱移位和显示技术，直观展示了傅立叶变换的应用效果。

数字图像处理与Python实现-图像变换-傅立叶变换

• 数字图像处理与Python实现-图像变换-傅立叶变换

  • 1. 前言

  • 2. 傅立叶变换定义及基本概念

  • 3. 傅立叶变换的特性

  • 3. 离散傅立叶变换

    • 3.1 离散傅立变换的定义
    • 3.2 离散傅立叶变换的性质

  • 4. 快速傅立叶变换（FFT）

  • 5.二维离散傅立叶变换

  • 6. 图像处理FFT实现

1. 前言

傅立叶变换被公认为强大的数学工具，它归类于正交变换。傅立叶变换在一维信号处理中得到了广泛应用。傅立叶变换在图像处理中同样也是非常有效的工具。

2. 傅立叶变换定义及基本概念

在数学领域，傅立叶变换具有明确的定义。定义为f(x)的函数，其中f(x)满足狄里赫莱条件：

(1) 包含有有限个不连续点；
(2) 包含有有限个极限点；
(3) 如果函数绝对可积，则可得出以下两式成立：

F(u) = \int_{-\infty}^{+\infty}{f(x)e^{-j2\pi ux}}dx \tag{2-1}

f(x) = \int_{+\infty}^{-\infty}{F(u)e^{j2\pi ux}}du \tag{2-2}

在时域和频域之间进行转换时，设\omega=2\pi u，从而将上述两个公式转化为：F(u) = \int_{-\infty}^{+\infty}{f(x)e^{-j\omega x}}dx \tag{2-3}

f(x) = \int_{+\infty}^{-\infty}{F(u)e^{j\omega x}}du \tag{2-4}

这些公式(2-3)和(2-4)通常被称为傅立叶变换对。函数f(x)的傅立叶变换通常是一个复量，其表达式如下：F(\omega) = R(\omega) + jI(\omega) \tag{2-5}

或写成指数形式：

F(\omega) = |F(\omega)|e^{\phi(\omega)} \tag{2-6}

|F(\omega)| = \sqrt{R^2(\omega) + I^2(\omega)} \tag{2-7}

\phi(\omega) = \arctan{\frac{I(\omega)}{R(\omega)}} \tag{2-8}

其中，|F(\omega)|为f(x)的傅立叶幅度谱，\phi(\omega)为相位谱。

该变换可轻易扩展至二维函数。满足狄利克雷条件的f(x,y)的二维傅立叶变换可表示为：

其中，F(u,v)表示变换后的函数。

F(u,v) = \int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}{f(x,y)e^{-j2\pi(ux+vy)}}dxdy \tag{2-9}

f(x,y) = \int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}{F(u,v)e^{j2\pi(ux+vy)}}dudv \tag{2-10}

二维傅立叶变换的幅度和相位如下：

|F(u,v)| = \sqrt{R^2(u,v) + I^2(u,v)} \tag{2-11}

\phi(u,v) = \arctan{\frac{I(u,v)}{R(u,v)}} \tag{2-12}

E(u,v) = R^2(u,v) + I^2(u,v) \tag{2-13}

其中，F(u,v)为幅度谱，\phi(u,v)为相伴谱，E(u,v)为能量谱。

3. 傅立叶变换的特性

(1). 具有可分离性：
F(u,v)可以表示为：
\begin{aligned} F(u,v) &= \int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}{f(x,y)e^{-j2\pi(ux+vy)}}dxdy \\ &=\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}{f(x,y)e^{-j2\pi(ux)} \cdot e^{-j2\pi(vx)} }dxdy \\ & = \int_{-\infty}^{+\infty} \lbrack \int_{-\infty}^{+\infty}{f(x,y)e^{-j2\pi ux}}\rbrack \cdot e^{-j2\pi(vx)} \\ & = \int_{-\infty}^{+\infty} \left\{ \mathcal{F}_x[f(x,y)] \right\} \cdot e^{-j2\pi(vx)} dy \\ & = \mathcal{F}_y \{ \mathcal{F}_x[f(x,y)]\} \end{aligned}
该二维傅里叶变换F(u,v)可以表示为两个一维傅里叶变换的组合。

(2). 线性 ：
\mathcal{F}[a_1f_1(x,y) + a_2f_2(x,y)] = a_1\mathcal{F}[f_1(x,y)] + a_2\mathcal{F}[f_2(x,y)]

共轭对称性质：其傅立叶变换为F(u,v)，其共轭函数为F^{*}(-u,-v)，则有：F(u,v) = F^*(-u,-v)。

旋转特性：当空间域函数绕原点旋转θ₀角时，其傅里叶变换在变换域中也绕原点旋转了θ₀角，即：f(r, \theta + \theta_0) \iff F(k, \phi + \theta_0)。

(5). 比例变换性质：
对于傅里叶变换F(u,v)及其缩放因子a和b，有以下性质：
其傅里叶变换为aF(u,v)；
缩放后的函数f(ax,by)的傅里叶变换为1/(|ab|)F(u/a,v/b)。

(6).帕塞瓦尔(Parseval)定理：
这一性质也被称为能量守恒定理。当F(u,v)为f(x,y)的傅里叶变换时，
\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}{|f(x,y)|^2}dxdy = \int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}{|F(u,v)|^2}dudv
这表明变换前后系统的能量保持不变。

(7).相关定理：
若f(x)和g(x)分别表示一维时域中的函数，而f(x,y)和g(x,y)则为二维空域中的函数，那么，我们定义以下表达式为相关运算：
f(x) \circ g(x) = \int_{-\infty}^{+\infty}{f(\alpha)g(x+\alpha)}d\alpha
f(x,y) \circ g(x,y) = \int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}{f(\alpha,\beta)g(x+\alpha,y+\beta)}d\alpha d\beta
其中，运算\circ代表相关运算，由此可导出傅里叶变换的一个关键性质，即相关定理：
f(x,y) \circ g(x,y) \iff F(u,v) \cdot G^*(u,v)
f(x,y) \cdot g^*(x,y) \iff F(u,v) \circ G(u,v)
此处，F(u,v)表示f(x,y)的傅里叶变换，G(u,v)表示g(x,y)的傅里叶变换，G^*(u,v)为G(u,v)的共轭复数，g^*则表示g(x,y)的共轭复数。

(8).卷积定理：
对于一维时域函数f(x)和g(x)，以及二维空域函数f(x,y)和g(x,y)，我们定义卷积运算的数学表达式为：
f(x) * g(x) = \int_{-\infty}^{+\infty}{f(\alpha)g(x-\alpha)}d\alpha，其中，\alpha为积分变量。

在公式中，卷积操作由*符号表示。基于此，我们可得，傅里叶域中的卷积定理如下：f(x,y) * g(x,y) \iff F(u,v) \cdot G(u,v)。

f(x,y) \cdot g(x,u) \iff F(u,v) * G(u,v)
式中，F(u,v)和G(u,v)分别是f(x,y)和g(x,y)的傅立叶变换。

3. 离散傅立叶变换

3.1 离散傅立变换的定义

假设x(n)为一离散时间序列，则其DFT正变换定义式如下：
X(m) = \sum_{n=0}^{N-1}{x(n)e^{-j \frac{2\pi mn}{N}}}
其反变换定义式如下：
x(n) = \frac{1}{N}\sum_{m=0}^{N-1}{e^{j\frac{2\pi mn}{N}}}
由此可见，DFT是一种直接处理离散时间信号的傅立叶变换工具。由此可知，

• (1)、连续时间周期信号:处理时间连续并且具有周期性的信号，其频域上离散，非周期。
• (2)、连续时间非周期信号:处理时间连续但是不具有周期性的信号，其频域上连续，非周期。
• (3)、离散时间非周期信号:处理时间离散，不具有周期性的信号，其频域上连续，有周期性。
• (4)、离散时间周期信号:处理时间离散，具有周期性的信号，其频域上离散，有周期性。
  离散傅立叶变换也可以使用如下公式表示：
  如果x(t)是连续的非周期函数，其频谱X(f)就是连续的非周期谱，对周期函数来说，其傅立叶变换可以写成：
  X(m) = \frac{1}{T}\int_{-\frac{T}{2}}^{+\frac{T}{2}}{x(t)e^{-j2\pi mt}}dt

x(t)的表示式为X(m)乘以指数函数的总和。
当x(t)为离散函数时，其频谱分析结果表现为一个周期性的连续频谱。
可以表示为：X(f) = \sum_{n=-\infty}^{\infty}{x(n)e^{-j2\pi nf}}

x(n) = \frac{1}{f}\int_{-\frac{f}{2}}^{+\frac{f}{2}}{X(f)e^{j2\pi nf}}df

基于以上两种情况，若函数x(n)既为离散型又为周期型，则其傅里叶变换必然表现为离散型周期频谱，具体表达式如下：
对于连续时间函数的傅里叶变换，我们采用求和运算来代替积分运算，从而得到离散傅里叶变换（DFT）公式：
X(k) = \sum_{n=0}^{N-1}{x(n)e^{-j2\pi kn/N}}
通过逆变换可以恢复原信号，其表达式为：
x(n) = \frac{1}{N}\sum_{k=0}^{N-1}{X(k)e^{j2\pi kn/N}}
这一过程表明，函数由连续向离散方向转变，使得傅里叶变换从积分运算转化为求和运算。

3.2 离散傅立叶变换的性质

• (1). 线性 ：
  如果时间序列x(n)与y(n)各有傅立叶变换X(m)和Y(m)，则
  ax(n) + by(n) \iff aX(m) + bY(m)

• (2) 对称性 ：
  如果：x(n) \iff X(m)
  则 \frac{1}{N}X(n) \iff x(-m)

• (3). 时间移位 ：
  如果x(n) \iff X(m)，序列向左或向右移动k位，则
  x(n-k) \iff X(m) \cdot W^{km}

• (4). 频率移位 ：
  如果x(n) \iff X(m)，
  则 x(n) \cdot W^{-km} \iff X(m-k)

• (5). 周期性 ：
  如果x(n) \iff X(m)，
  则 x(n\pm rN) = x(n)

• (6). 偶函数 ：
  如果x_e(n) = x_e(-n)，
  则 X_e(m) = \sum_{n=0}^{N-1}{x_e(n)[\cos(\frac{2\pi}{N}mn)^{-jsin(\frac{2\pi}{N}mn)}]}

• (7). 奇函数：
  当满足x₀(n) = -x₀(-n)时，
  复数形式的傅里叶系数X₀(m)可表示为：
  X₀(m) = -jΣₙ₌₀^{N-1} x₀(n)·sin(2πmn/N)。
  考虑到x₀(n)为奇函数，而cos(2πmn/N)为偶函数，
  由此可知，实部之和为零，因此：
  X₀(m) = -jΣₙ₌₀^{N-1} x₀(n)·sin(2πmn/N)。

(8).卷积定理：
根据卷积的定义可知：
x(n) * y(n) = \sum_{h=0}^{N-1}{x(h)y(n-h)} 则，
x(n)与y(n)的傅里叶变换的乘积等于X(m)与Y(m)的卷积。
x(n)与y(n)的乘积的傅里叶变换等于X(m)与Y(m)的卷积。

• (9). 相关定理 ：
  如果x(n) \iff X(m)，y(n) \iff Y(m) 则，
  x(n) \circ y(n) \iff X^*(m) \cdot Y^*(m)

• (10). 帕斯维尔定理 ：
  如果x(n) \iff X(m)，则，
  \sum_{n=0}^{N-1}{x^2(n)} = \frac{1}{N}\sum_{m=0}^{N-1}|X(m)|^2

4. 快速傅立叶变换（FFT）

离散傅立叶变换已经成为数字信号处理的重要工具，然而，它的计算量比较大，运算时间长。快速傅立叶变换（FFT）并不是一种新的变换，它是离散傅立叶变换的一种算法。
对于一个有限长序列{x(n)}(0 \leq n \leq N - 1)，它的傅立叶变换表示如下：
X(m) = \sum_{n=0}^{N-1}{x(n)e^{\frac{-j2\pi nm}{N}}} (其中，m=0,1,...,N-1)
因此，傅立叶变换对可以写成：
X(m) = \sum_{n=0}^{N-1}{x(n)W^{nm}}
x(n) = \frac{1}{N}\sum_{m=0}^{N-1}{X(m)W^{-nm}}
以上式子就不再做展开，展开后可以知道，离散傅立叶变换中的乘法运算有许多重复内容，因此离散傅立叶变换可以写成：
X(m) = \sum_{n=0}^{N-1}{x(n)W_{N}^{nm}}
由周期性可以得到：
X(m) = X_1(m) + W_{N}^{m}X_2(m)，m=0,1,...,N-1

通过分析可知，一个N点的离散傅立叶变换可通过两个N/2点的傅立叶变换实现。离散傅立叶变换的计算时间主要取决于乘法操作，其分解后的乘法运算次数将大幅减少。快速傅立叶变换（FFT）主要分为两类：一类是基于时间域的分解方法，另一类是基于频率域的分解方法。FFT的具体实现细节在此不做深入讨论，若对此感兴趣，建议查阅相关专业文献。

5.二维离散傅立叶变换

一幅静止的图像可以看成是二维数据阵列。数字图像处理主要是二维数据处理。二维离散傅立叶变换的定义如下：
正变换为：
F(u,v) = \sum_{x=0}^{M-1}\sum_{y=0}^{N-1}{f(x,y)\exp[-j2\pi(\frac{ux}{M}+\frac{uy}{N})]} (u=0,1,2,...,M-1,v=0,1,2,...,N-1)
反变换为：
F(x,y) = \frac{1}{MN}\sum_{u=0}^{M-1}\sum_{v=0}^{N-1}{F(u,v)\exp[j2\pi(\frac{ux}{M} + \frac{vy}{N})]}(x=0,1,2,...,M-1,y=0,1,2,...,N-1)
式中，符号F(u,v)可称为空间频率。

二维离散傅立叶变换具有以下性质：

• 平移性
• 旋转性

6. 图像处理FFT实现

    import numpy as np
    import scipy
    import cv2
    
    src = np.zeros((256,256),np.uint8);
    
    src[124:132,120:136] = 255
    
    # src = cv2.imread('resources/images/circle3.jpg',0)
    
    srcf = src.astype(np.float32) / 255.0
    
    # 二维FFT变换
    srcfft = np.fft.fft2(srcf)
    
    # 转换成可以显示的数据
    srcabs = np.log( np.abs(srcfft))
    srcabs = np.clip(srcabs * 255,0,255).astype(np.uint8)
     
    # 频谱移到中心
    srcfftshift = np.fft.fftshift(srcfft)
    
    srcfftshiftabs = np.log(np.abs(srcfftshift))
    srcfftshiftabs = np.clip(srcfftshiftabs * 255,0,255).astype(np.uint8)
    
    cv2.imshow('src',src);
    cv2.imshow('srcfft',srcabs)
    cv2.imshow('srcfftshift',srcfftshiftabs)
    
    cv2.waitKey();
    cv2.destroyAllWindows()

程序运行结果：