1. 连续傅里叶变换
1.1 一维连续傅里叶变换
-
一维连续傅里叶正变换(\(\text{1-Dimensional Continuous Fourier Transform}\))
对于函数 \(f(x)\),一维连续傅里叶变换有如下定义:
\[\Re: \; F(u) = \int_{- \infty}^{\infty} f(x)e^{-j2 \pi u x} \text{d}x \]其中 \(j^2 = 1\).
-
一维连续傅里叶反变换(\(\text{1-Dimensional Continuous Fourier Inversion Transform}\))
上式中 \(F(u)\) 的反变换,有如下定义:
\[\Re ^{-1}: \; f(x) = \int_{- \infty}^{\infty} f(x)e^{j2 \pi u x} \text{d}u \]
函数 \(f(x)\) 与 \(F(u)\) 被称为一个 傅里叶变换对。对于任意函数 \(f(x)\),其傅里叶变换 \(F(u)\) 是 唯一的;反之亦然。
对于函数 \(f(x)\),其对应的 \(F(u)\) 一般是复函数。由 欧拉公式 \(e^{ix} = \text{cos}x + i\;\text{sin}x\) 有如下定义:
-
实部
\[R(u) = \int_{- \infty}^{\infty} f(x)\text{cos}(2\pi u x)\;\text{d}x \] -
虚部
\[I(u) = - \int_{- \infty}^{\infty} f(x)\text{sin}(2\pi u x)\;\text{d}x \] -
振幅
\[|F(u)| = [R^2(u) + I^2(u)]^{\frac{1}{2}} \] -
能量
\[E(u) = |F(u)|^2 = R^2(u) + I^2(u) \] -
相位
\[\varphi (u) = \text{arctan}\frac{I(u)}{R(u)} \]
1.2 二维连续傅里叶变换
设 \(f(x, y)\) 连续可积,且 \(F(u, v)\) 可积,则存在傅里叶变换对:
-
正变换
\[F\{f(x, y)\} = F(u, v) = \int_{- \infty}^{\infty}f(x, y)e^{-j2\pi (ux + vy)}\; \text{d}x \text{d}y \] -
反变换
\[F^{-1}\{F(u, v)\} = f(x, y) = \int_{- \infty}^{\infty}F(u, v)e^{j2\pi (ux + vy)}\; \text{d}u \text{d}v \]
同样有如下定义:
-
傅里叶频谱
\[|F(u, v)| = [R^2(u, v) + I^2(u, v)]^{\frac{1}{2}} \] -
能量谱
\[E(u, v) = |F(u, v)|^2 = R^2(u, v) + I^2(u, v) \] -
相位谱
\[\varphi (u, v) = \text{arctan}\frac{I(u, v)}{R(u, v)} \]
2. 离散傅里叶变换
2.1 一维离散傅里叶变换
\(f(x)\) 是在时域上等距采样得到的 \(N\) 点离散序列, \(x\) 是离散实变量,\(u\) 为离散频率变量。由此,离散傅里叶变换(\(\text{Discrete Fourier Transform, DFT}\))定义如下:
-
正变换
\[\Re: \; F(u) = \frac{1}{\sqrt{N}}\sum_{x=0}^{N-1}f(x)e^{-j2 \pi \frac{ux}{N}}, \; u = 0, 1, ... ,N - 1 \] -
反变换
\[\Re: \; f(x) = \frac{1}{\sqrt{N}}\sum_{x=0}^{N-1}F(u)e^{j2 \pi \frac{ux}{N}}, \; u = 0, 1, ... ,N - 1 \]
由欧拉公式,可以得到如下定义:
-
傅里叶频谱
\[|F(u)| = [R^2(u) + I^2(u)]^{\frac{1}{2}} \] -
能量谱
\[E(u) = |F(u)|^2 = R^2(u) + I^2(u) \] -
相位谱
\[\varphi (u) = \text{arctan}\frac{I(u)}{R(u)} \]
2.2 二维离散傅里叶变换
同连续函数的傅里叶变换一样,离散函数的傅里叶变换也可推广到二维的情形,其二维离散傅里叶变换定义为:
-
正变换
\[F(u, v) = \frac{1}{\sqrt{MN}}\sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x, y) e^{-j 2 \pi (\frac{ux}{M} + \frac{vy}{N})}, \; u = 0, 1, ... , M - 1, \; v = 0, 1, ... , N - 1 \] -
反变换
\[f(x, y) = \frac{1}{\sqrt{MN}}\sum_{x=0}^{M-1} \sum_{y=0}^{N-1} F(u, v) e^{j 2 \pi (\frac{ux}{M} + \frac{vy}{N})}, \; u = 0, 1, ... , M - 1, \; v = 0, 1, ... , N - 1 \]
同样有如下定义:
-
傅里叶频谱
\[|F(u, v)| = [R^2(u, v) + I^2(u, v)]^{\frac{1}{2}} \] -
能量谱
\[E(u, v) = |F(u, v)|^2 = R^2(u, v) + I^2(u, v) \] -
相位谱
\[\varphi (u, v) = \text{arctan}\frac{I(u, v)}{R(u, v)} \]
频谱:图像的 主要成分是低频信息,它形成了图像的 基本灰度等级,对图像结构的决定作用较小;
中频信息 决定了图像的基本结构,形成了图像的 主要边缘结构;
高频信息 形成了图像的 边缘和细节,是在中频信息上对图像内容的进一步强化。
3. 快速傅里叶变换
3.1 FFT的作用
快速傅里叶变换 (\(\text{Fast Fourier Transform, FFT}\))并不是一种新的变换,它是 离散傅里叶变换(\(\text{DFT}\)) 的一种算法。
这种方法是在分析 离散傅里叶变换(\(\text{DFT}\)) 中的多余运算的基础上,进而消除这些重复工作的思想指导下得到的,所以在运算中大大节省了工作量,达到了快速的目的。
对于一个有限长序列 \(\{f(x)\}, 0 \le x \le n - 1\),其傅里叶变换如下:
\[F(u) = \frac{1}{\sqrt{n}} \sum_{x=0}^{n - 1}f(x) w_n^{ux} \]其中 \(W_n = e^{-j\frac{2 \pi}{n}}\)
计算每个频率的分量,需要进行 \(n\) 次运算,计算所有频率的分量,则需要进行 \(O(n^2)\) 的时间复杂度。当 \(n\) 较大时,这显然会耗费很多时间。
而 快速傅里叶变换 (\(\text{FFT}\)) 可以做到在 \(O(n\text{log}n)\) 的时间复杂度,计算傅里叶变换对应的多项式。
3.2 单位根
以单位圆点为起点,单位圆的 \(n\) 等分点为终点,在单位圆上可以得到 \(n\) 个复数,设幅角为正且最小的复数为 \(w_n\),称为 \(n\) 次 单位根,即:
\[w_n = \text{cos}\frac{2 \pi}{n} + i \;\text{sin} \frac{2 \pi}{n} \]由欧拉公式可得:
\[w_n^k = \text{cos}\frac{2 k \pi}{n} + i \;\text{sin} \frac{2 k \pi}{n} \]特别地,\(w^0_n = w^n_n = 1\)
由周期性可以得到如下性质:
\[\begin{split} w_{cn}^{ck} &= \text{cos}\frac{2 ck \pi}{cn} + i \;\text{sin} \frac{2 ck \pi}{cn} \\\\ &= w^k_n \\\\\\ w_{n}^{k + \frac{n}{2}} &= w^k_n \cdot [\text{cos}(\frac{2 \pi}{n} \cdot \frac{n}{2}) + i \;\text{sin} (\frac{2 \pi}{n} \cdot \frac{n}{2})] \\\\ &= w^k_n \cdot (\text{cos} \pi + i \; \text{sin} \pi) \\\\ & = - w^k_n \end{split} \]3.3 快速傅里叶变换
对于傅里叶变换计算表达式
\[F(k) = \frac{1}{\sqrt{n}} \sum_{x=0}^{n - 1}f(x) w_n^{kx} \]\[w_n^{kx} = e^{-j\frac{2 \pi}{n}kx} \]我们令 \(t = w_n^k, a_x = f(x)\)
由此可化简为:
\[F(k) = \frac{1}{\sqrt{n}} \sum_{x=0}^{n - 1}a_x t^k \]展开傅里叶变换表达式:
\[\begin{split} F(t) &= (a_0 + a_2 t^2 + ... + a_{n-2}t^{n-2}) \\ &+ (a_1t + a_3t^{3} + ... + a_{n-1}t^{(n-1)}) \end{split} \]令
\[\begin{split} G(t) = a_0 + a_2t + a_4t^{2} + ... + a_{n-2}t^{\frac{n}{2} - 1} \\\\ H(t) = a_1 + a_3t + a_5t^{2} + ... + a_{n-1}t^{\frac{n}{2} - 1} \end{split} \]由此可得:
\[F(t) = G(t^2) + t\; H(t^2) \]将 \(w^k_n = t\) 回代到等式中,由 单元根 性质可得:
\[\begin{split} F(w^k_n) &= G(w_n^{2k}) + w^k_n \cdot H(w_n^{2k}) \\\\ &= G(w^k_{\frac{n}{2}}) + w^k_n \cdot H(w^k_{\frac{n}{2}}) \\\\\\ F(w^{k + \frac{n}{2}}_n) &= G(w^{2k + n}_n) + w^{k + \frac{n}{2}}_n \cdot H(w_n^{2k + n}) \\\\ &= G(w^{2k}_n \cdot w^n_n) - w_n^k \cdot H(w^{2k}_n \cdot w^n_n) \\\\ &= G(w^{2k}_n) - w^k_n \cdot H(w^{2k}_n) \\\\ &= G(w_{\frac{n}{2}}^k) - w^k_n \cdot H(w_{\frac{n}{2}}^k) \end{split} \]可以发现,只要求出 \(G(w_{\frac{n}{2}}^k)\) 和 \(H(w_{\frac{n}{2}}^k)\) 就可以进而求出 \(F(w^k_n)\) 和 \(F(w_n^{k + \frac{n}{2}})\)。
对于 \(G(w_{\frac{n}{2}}^k)\) 和 \(H(w_{\frac{n}{2}}^k)\) 的求解,只需要转换求解为 \(G(w^k_{\frac{n}{4}})\) 和 \(H(w^k_{\frac{n}{4}})\) 的问题即可。
直到最后,转换为求解 \(G(w^k_1) = 1\) 和 \(H(w_1^k) = 1\) 的问题。
所以,只需要不断对 \(G, H\) 进行 递归 求解即可。 求解每个分量对应的 \(F(w_n^k)\) 只需要 \(O(\text{log} n)\) 的时间复杂度。求解所有分量则需要 \(O(n \; \text{log} n)\) 的时间复杂度。
涉猎过算法竞赛的都知道,\(O(n \; \text{log} n)\) 比 \(O(n ^ 2)\) 在 \(n\) 较大时会快很多。
4. 傅里叶变换的性质
4.1 可分离性
二维离散傅里叶变换对,可以写成如下形式:
\[\begin{split} & F(u, v) = \frac{1}{\sqrt{MN}}\sum_{x=0}^{M-1}e^{-j 2\pi \frac{ux}{M}}f(x, y)\sum_{y=0}^{N-1}e^{-j 2\pi \frac{vy}{N}} \\\\ & f(x, y) = \frac{1}{\sqrt{MN}}\sum_{u=0}^{M-1}e^{j2 \pi \frac{ux}{M}} \sum_{v=0}^{N-1} F(u, v) e^{j 2\pi \frac{vy}{N}} \\\\ & u = 0, 1, ... , M-1; \; v = 0, 1, ... , N - 1 \end{split} \]一个 \(\text{2-D}\) 傅里叶变换可由连续两次运用 \(1-D\) 傅里叶变换来实现:
\[\begin{split} F(u, y) &= \frac{1}{\sqrt{M}}\sum_{x=0}^{M-1}f(x, y) e^{-j 2 \pi \frac{ux}{M}} \\\\ F(u, v) &= \frac{1}{\sqrt{N}}\sum_{y=0}^{N-1}F(u, y) e^{-j 2 \pi \frac{vy}{N}} \end{split} \]当 \(M = N\) 时,变换过程如图所示:
4.2 平移性质
如果 \(F(u,v)\) 的频率变量 \(u\),\(v\) 各移动了 \(u_0\),\(v_0\) 距离,\(f(x, y)\) 的实变量 \(x\),\(y\) 各移动了 \(x_0\),\(y_0\)距离,则:
\[\begin{split} & f(x, y)e^{j 2 \pi (\frac{u_0x}{M} + \frac{v_0y}{N})} \Leftrightarrow F(u - u_0, v - v_0) \\\\ & f(x - x_0, y - y_0) \Leftrightarrow F(u, v)e^{-j 2 \pi (\frac{ux_0}{M} + \frac{vy_0}{N})} \end{split} \]特别地,当 \(u_0 = \frac{M}{2}, v_0 = \frac{N}{2}\) 时,\(e^{j 2 \pi (\frac{u_0x}{M} + \frac{v_0y}{N})} = e^{j \pi (x + y)} = (-1)^{x + y}\)
则:
\[\begin{split} & f(x, y)(-1)^{x+y} \Leftrightarrow F(u - \frac{M}{2}, v - \frac{N}{2}) \\\\ & f(x - \frac{M}{2}, y - \frac{N}{2}) \Leftrightarrow F(u, v)(-1)^{u + v} \end{split} \]这条性质被称为 移中性 。
如图所示,图像平移 不会改变图像的频谱形状(fftshift后),但是会 改变相位:
相关代码如下:
I = imread('3.png');
I = rgb2gray(I);
I = im2double(I);
F = fftshift(fft2(I));
subplot(2,3,1), imshow(I,[]), title('Original Image');
subplot(2,3,2), imshow(log(1+abs(F)),[]), title('Fourier Transform - Magnitude');
subplot(2,3,3), imshow(angle(F),[-pi pi]), title('Fourier Transform - Phase');
translation = [400, 400];
I_translate = circshift(I, translation);
F_translate = fftshift(fft2(I_translate));
subplot(2,3,4), imshow(I_translate,[]), title('Translated Image');
subplot(2,3,5), imshow(log(1+abs(F_translate)),[]), title('Fourier Transform of Translated Image - Magnitude');
subplot(2,3,6), imshow(angle(F_translate),[-pi pi]), title('Fourier Transform of Translated Image - Phase');
4.3 周期性和共轭对称性
傅里叶变换和反变换的周期性:
\[\begin{split} F(u, v) &= F(u + M, v) = F(u, v + N) = F(u + M, v + M) \\\\ f(x, y) &= f(x + M, y) = f(x, y + N) = f(x + M, y + N) \end{split} \]如果 \(f(x, y)\) 是实函数,则它的傅里叶变换具有共轭对称性:
\[\begin{split} F(u, v) &= F^*(-u, -v) \\\\ |F(u, v)| &= |F(-u, -v)| \end{split} \]其中 \(F^*(u, v)\) 为 \(F(u, v)\) 的 复共轭 。
4.4 旋转性质
借助极坐标 \(x = r \text{cos} \theta, y = r\text{sin}\theta, u = w\text{cos}\phi, v = w\text{sin}\phi\):
\[\begin{split} & f(x, y) \Leftrightarrow F(u, v) \\\\ & f(r\text{cos}\theta, r \text{sin}\theta) \Leftrightarrow F(w \text{cos}\phi, w \text{sin}\phi) \\\\ & f(r, \theta + \theta_0) \Leftrightarrow F(w, \phi + \theta_0) \end{split} \]如图所示,旋转之后,傅里叶频谱也随之旋转,相位也发生了改变。会 改变图像的频谱,也会 改变相位:
相关代码如下:
I = imread('3.png');
I = rgb2gray(I);
I = im2double(I);
F = fftshift(fft2(I));
subplot(2,3,1), imshow(I,[]), title('Original Image');
subplot(2,3,2), imshow(log(1+abs(F)),[]), title('Fourier Transform - Magnitude');
subplot(2,3,3), imshow(angle(F),[-pi pi]), title('Fourier Transform - Phase');
I_translate = imrotate(I, -45, 'bilinear', 'crop');
F_translate = fftshift(fft2(I_translate));
subplot(2,3,4), imshow(I_translate,[]), title('Translated Image');
subplot(2,3,5), imshow(log(1+abs(F_translate)),[]), title('Fourier Transform of Translated Image - Magnitude');
subplot(2,3,6), imshow(angle(F_translate),[-pi pi]), title('Fourier Transform of Translated Image - Phase');
4.5 分配律
傅里叶变换和反变换对 加法满足分配律,根据傅里叶变换对的定义可得到:
\[F\{f_1(x, y) + f_2(x, y)\} = F\{f_1(x, y)\} + F\{f_2(x, y)\} \]乘法则不满足 :
\[F\{f_1(x, y) \cdot f_2(x, y)\} \not = F\{f_1(x, y)\} \cdot F\{f_2(x, y)\} \]4.6 尺度变换
尺度变换描述了函数自变量的尺度变化对其傅里叶变换的作用,有如下表达式:
\[af(x, y) \Leftrightarrow aF(u, v) \]可证明得到:
\[f(ax, by) \Leftrightarrow \frac{1}{|ab|}F(\frac{u}{a}, \frac{v}{b}) \]如图所示:
相关代码:
img = imread('img.jpg');
img_norm = im2double(img);
figure;
subplot(2, 2, 1);
imshow(img_norm);
title('initial image');
img_scaled = img_norm * 0.1;
subplot(2, 2, 2);
imshow(img_scaled);
title('scaled image');
fft_img = fft2(img_norm);
fft_img = fftshift(fft_img);
subplot(2, 2, 3);
imshow(log(1 + abs(fft_img)), []);
title('initial image fft');
fft_scaled_img = fft2(img_scaled);
fft_scaled_img = fftshift(fft_scaled_img);
subplot(2, 2, 4);
imshow(log(1 + abs(fft_scaled_img)), []);
title('scaled image fft');
4.7 平均值
对于一个二维离散函数,其平均值为:
\[\bar f(x, y) = \frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x, y) \]傅里叶变换在原点 \((0, 0)\) 的频谱分量:
\[\begin{split} F(0, 0) &= \frac{1}{\sqrt{MN}}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x, y)e^{-j2 \pi 0} \\\\ &= \sqrt{MN}[\frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x, y)] \\\\ &= \sqrt{MN} \bar f(x, y) \end{split} \]4.8 卷积定理
-
一维傅里叶变换
\[f(x) * g(x) = \int_{- \infty}^{\infty}f(z)g(x - z)\text{d}z \Leftrightarrow F(u)G(u) \] -
二维傅里叶变换
\[f(x, y) * g(x, y) \Leftrightarrow F(u, v)G(u, v) \]