快速傅里叶变换

作者: 寽虎非虫003 | 来源:发表于2021-07-16 18:29 被阅读0次

零、序言

0.1主要参考资料

[1]timothy sauer，裴玉茹.《数值分析(第二版)》[M]. 北京: 机械工业出版社, 2014.10 : P467~P415.
[2]冈萨雷斯.《数字图像处理(第四版)》[M]. 北京: 电子工业出版社, 2020.5: P138~P212.
[3]丘维声. 《高等代数》[M]. 北京: ，： .

吐槽一句，简书的公示显示很不稳定，每次刷新都可能是不同的公式不能正常显示。

0.2目录

零、序言
一、复数
二、一维离散傅里叶变换
三、二维离散傅里叶变换
四、快速傅里叶变换
五、现有代码实验
六、cpu实现
七、gpu实现
八、总结

一、复数

复数表示： $z=a+bi$
模： $|z| = \sqrt{a^2+b^2}$
共轭复数： $\overline{z}=a-bi$

欧拉公式： $e^{i\theta}=\cos{\theta}+i\sin{\theta}$
$z = e^{i\theta}$ 在复平面上的单位圆上。
任一复数 $z = a+bi$ 都可以用极坐标形式表示成 $z=a+bi=re^{i\theta}$ ，其中 $r=|z|= \sqrt{a^2+b^2}, \theta = \arctan{\frac{b}{a}}$ 。

先复习下三角函数的公式：
$\cos{\alpha + \beta} = \cos{\alpha}\cos{\beta} - \sin{\alpha}\sin{\beta}\\ \sin{\alpha + \beta} = \sin{\alpha}\cos{\beta} + \cos{\alpha}\sin{\beta}$
复数乘法：
$\begin{equation} \begin{array}{ll} e^{i\theta}e^{i\gamma} & = & (\cos{\theta}+i\sin{\theta})(\cos{\gamma}+i\sin{\gamma})\\ & = & \cos{\theta}\cos{gamma}+i(\sin{\theta}\cos{\gamma}+\cos{\gamma}\sin{\theta})\\ & = & \cos{(\theta + \gamma )} + i\sin{(\theta + \gamma )}\\ & = &e^{i(\theta +\gamma )} \end{array} \end{equation}$

本原单位根： $w=e^{-i\frac{2\pi}{n}}$ 是 $n$ 次本原单位根。
$\Rightarrow w^n = 1\\$
$\begin{array}{ll} \Rightarrow & 1+w + w^1 + w^2 + \cdots +w^{n-1} & = 0\\ & 1+w^2 + w^4 + w^6 + \cdots +w^{2(n-1)} &= 0\\ & 1+w^3 + w^6 + w^9 + \cdots +w^{3(n-1)} &= 0\\ & & \vdots &\\ & 1+w^{n-1} + w^{(n-1) 2} + w^{(n-1) 3} + \cdots +w^{(n-1)^2} &= 0\\ \end{array}$
$\Rightarrow \sum_{j=0}^{n-1}{w^{jk}} = \left\{ \begin{array}{ll} n, & k,\frac{k}{n}\in Q\\ 0, & others \end{array}\right.$

二、一维离散傅里叶变换

离散傅里叶变换定义： $x=[x_0,x_1,\cdots,x_{n-1}]^T$ 的离散傅里叶变换(DFT)为 $n$ 维向量 $y=[y_0,y_1,\cdots ,y_{}n-1]^T$ ，其中 $w=e^{-i\frac{2\pi}{n}}$
$y_k=\frac{1}{\sqrt{n}}\sum_{j=0}^{n-1}{x_j w^{jk}}$
用矩阵表示如下
$\left[ \begin{array}{ll} & y_0\\ & y_1\\ & y_2\\ & y_3\\ & \vdots \\ & y_{n-1} \end{array}\right] = \left[ \begin{array}{ll} a_0+ib_0\\ a_1+ib_1\\ a_2+ib_2\\ a_3+ib_3\\ \vdots \\ a_{n-1} + ib_{n-1} \end{array}\right] =\frac{1}{\sqrt{n}} \left[ \begin{array}{ll} w^0 & w^0 & w^0 & \cdots & w^0\\ w^0 & w^1 & w^2 & \cdots & w^{n-1}\\ w^0 & w^2 & w^4 & \cdots & w^{2(n-1)}\\ w^0 & w^3 & w^6 & \cdots & w^{3(n-1)}\\ \vdots & \vdots & \vdots & & \vdots\\ w^0 & w^{n-1} & w^{2(n-1)} & \cdots & w^{(n-1)^2}\\ \end{array}\right] \left[ \begin{array}{ll} x_0\\ x_1\\ x_2\\ x_3\\ \vdots \\ x_n-1 \end{array}\right]$
上式中的 $n \times n$ 矩阵称为傅里叶矩阵，傅里叶矩阵是一个范德蒙矩阵。

傅里叶变换的矩阵为：
$F_n = \frac{1}{\sqrt{n}} \left[ \begin{array}{ll} w^0 & w^0 & w^0 & \cdots & w^0\\ w^0 & w^1 & w^2 & \cdots & w^{n-1}\\ w^0 & w^2 & w^4 & \cdots & w^{2(n-1)}\\ w^0 & w^3 & w^6 & \cdots & w^{3(n-1)}\\ \vdots & \vdots & \vdots & & \vdots\\ w^0 & w^{n-1} & w^{2(n-1)} & \cdots & w^{(n-1)^2}\\ \end{array}\right]$
傅里叶逆变换的矩阵为
$F_n^{-1} = \overline{F_n} =\frac{1}{\sqrt{n}} \left[ \begin{array}{ll} w^0 & w^0 & w^0 & \cdots & w^0\\ w^0 & w^{-1} & w^{-2} & \cdots & w^{-(n-1)}\\ w^0 & w^{-2} & w^{-4} & \cdots & w^{-2(n-1)}\\ w^0 & w^{-3} & w^{-6} & \cdots & w^{-3(n-1)}\\ \vdots & \vdots & \vdots & & \vdots\\ w^0 & w^{-(n-1)} & w^{-2(n-1)} & \cdots & w^{- (n-1)^2}\\ \end{array}\right]$

若方阵 $F$ 满足 $\overline{F}^TF = I$ ，则称 $F$ 是酉阵。傅里叶矩阵是酉阵。

例题就不搬了。

复数共轭和乘积运算时可交换的，即先共轭再相乘等于先相乘再共轭。

三、二维离散傅里叶变换

正式开始之前先吐槽一句，冈萨雷斯是真的难，太难了，看一个小时推导就头疼得要死。
顺便说一句，看推导的话，还请移步冈萨雷斯。

二维离散傅里叶变换：
$F_{(u,v)} = \frac{1}{\sqrt{mn}} \sum_{x=0}^{m-1}\sum_{y=0}^{n-1}f_{(x,y)}e^{-i2\pi(\frac{ux}{m} + \frac{vy}{n})}$
二维离散傅里叶逆变换：
$f_{(x,y)} = \frac{1}{\sqrt{mn}} \sum_{x=0}^{m-1}\sum_{y=0}^{n-1}F_{(u,v)}e^{i2\pi(\frac{ux}{m} + \frac{vy}{n})}$
以上写法和冈萨雷斯有所出入，是为了和《数值分析》的表达形式统一。
而且这儿暂时不太好表示成 $w$ 的次方的形式，也不是，不好表示，而是若是要表示的话需要特地说明两个 $w_1$ 和 $w_2$ 分别是几次本原单位根，这步还是放到后面说。

运算可分性：
二维傅里叶变换：
$\begin{array}{ll} F_{(u,v)} & = & \frac{1}{ \sqrt{mn}} \sum_{x=0}^{m-1} \sum_{y=0}^{n-1} f_{(x,y)}e^{-i2 \pi (\frac{ux}{m} + \frac{vy}{n})} \\ & = & \frac{1}{\sqrt{mn}} \sum_{x=0}^{m-1} e^{-i \frac{2\pi ux}{m}} \sum_{y=0}^{n-1} f_{(x,y)}e^{-i \frac{2\pi vy }{n}} \\ & = & \frac{1}{\sqrt{mn}} \sum_{x=0}^{m-1} ( \sum_{y=0}^{n-1}f_{(x,y)}w_n^{vy}) w_m^{ux} \\ \end{array}$
其中 $w_m = e^{-i\frac{2 \pi }{m}}$ ， $w_n = e^{-i \frac{w \pi }{n}}$ ，下同。
同理，二维傅里叶逆变换：
$f_{(x,y)} = \frac{1}{\sqrt{mn}} \sum_{x=0}^{m-1} (\sum_{y=0}^{n-1}F_{(u,v)} w_n^{-vy})w_m^{-ux}$
这一部分理论来源于

矩阵运算形式：千万要注意转置：
二维离散傅里叶变换矩阵形式：
$F_{m \times n} = \frac{1}{\sqrt{mn}} \left[ \begin{array}{ll} w_m^0 & w_m^0 & w_m^0 & \cdots & w_m^0 \\ w_m^0 & w_m^1 & w_m^2 & \cdots & w_m^{m-1} \\ w_m^0 & w_m^2 & w_m^4 & \cdots & w_m^{2(m-1)} \\ w_m^0 & w_m^3 & w_m^6 & \cdots & w_m^{3(m-1)} \\ \vdots & \vdots & \vdots & & \vdots \\ w_m^0 & w_m^{m-1} & w_m^{2(m-1)} & \cdots & w_m^{(m-1)^2} \\ \end{array}\right]_{m \times m } \left( \left[ \begin{array}{ll} w_n^0 & w_n^0 & w_n^0 & \cdots & w_n^0 \\ w_n^0 & w_n^1 & w_n^2 & \cdots & w_n^{n-1} \\ w_n^0 & w_n^2 & w_n^4 & \cdots & w_n^{2(n-1)} \\ w_n^0 & w_n^3 & w_n^6 & \cdots & w_n^{3(n-1)} \\ \vdots & \vdots & \vdots & & \vdots \\ w_n^0 & w_n^{n-1} & w_n^{2(n-1)} & \cdots & w_n^{(n-1)^2}\\ \end{array}\right]_{n \times n } \left[ \begin{array}{ll} f_{(0,0)} & f_{(0,1)} & f_{(0,2)} & \cdots & f_{(0,n-1)} \\ f_{(1,0)} & f_{(1,1)} & f_{(1,2)} & \cdots & f_{(1,n-1)} \\ f_{(2,0)} & f_{(2,1)} & f_{(2,2)} & \cdots & f_{(2,n-1)} \\ f_{(3,0)} & f_{(3,1)} & f_{(3,2)} & \cdots & f_{(3,n-1)} \\ \vdots & \vdots & \vdots & & \vdots \\ f_{(m-1,0)} & f_{(m-1,1)} & f_{(m-1,2)} & \cdots & f_{(m-1,n-1)} \end{array}\right]_{m \times n }^{T} \right)^{T}$
二维离散傅里叶逆变换矩阵形式：
$f_{m \times n} = \frac{1}{\sqrt{mn}} \left[ \begin{array}{ll} w_m^0 & w_m^0 & w_m^0 & \cdots & w_m^0 \\ w_m^0 & w_m^{-1} & w_m^{-2} & \cdots & w_m^{-(m-1)} \\ w_m^0 & w_m^{-2} & w_m^{-4} & \cdots & w_m^{-2(m-1)} \\ w_m^0 & w_m^{-3} & w_m^{-6} & \cdots & w_m^{-3(m-1)} \\ \vdots & \vdots & \vdots & & \vdots \\ w_m^0 & w_m^{-(m-1)} & w_m^{-2(m-1)} & \cdots & w_m^{-(m-1)^2} \\ \end{array}\right]_{m \times m } \left( \left[ \begin{array}{ll} w_n^0 & w_n^0 & w_n^0 & \cdots & w_n^0 \\ w_n^0 & w_n^{-1} & w_n^{-2} & \cdots & w_n^{-(n-1)} \\ w_n^0 & w_n^{-2} & w_n^{-4} & \cdots & w_n^{-2(n-1)} \\ w_n^0 & w_n^{-3} & w_n^{-6} & \cdots & w_n^{-3(n-1)} \\ \vdots & \vdots & \vdots & & \vdots \\ w_n^0 & w_n^{-(n-1)} & w_n^{-2(n-1)} & \cdots & w_n^{-(n-1)^2} \\ \end{array}\right]_{n \times n } \left[ \begin{array}{ll} F_{(0,0)} & F_{(0,1)} & F_{(0,2)} & \cdots & F_{(0,n-1)} \\ F_{(1,0)} & F_{(1,1)} & F_{(1,2)} & \cdots & F_{(1,n-1)} \\ F_{(2,0)} & F_{(2,1)} & F_{(2,2)} & \cdots & F_{(2,n-1)} \\ F_{(3,0)} & F_{(3,1)} & F_{(3,2)} & \cdots & F_{(3,n-1)} \\ \vdots & \vdots & \vdots & & \vdots \\ F_{(m-1,0)} & F_{(m-1,1)} & F_{(m-1,2)} & \cdots & F_{(m-1,n-1)} \end{array}\right]_{m \times n }^{T} \right)^{T}$
也就是说，，二维傅里叶变换可以拆分成两步一维傅里叶变换进行运算。

写对也不容易啊，东西基本都忘得差不多了。

四、快速傅里叶变换

先看一个一维傅里叶变换的例子，计算 $z=M_nx=\sqrt{n}F_nx$ .
这儿本原单位根 $w=w_4=e^{-i \frac{2 \pi }{4}} = -i$ 。
$\left[ \begin{array}{ll} z_0\\ z_1\\ z_2\\ z_3\\ \end{array}\right] = \left[ \begin{array}{ll} w^0 & w^0 & w^0 & w^0\\ w^0 & w^1 & w^2 & w^3\\ w^0 & w^2 & w^4 & w^6\\ w^0 & w^3 & w^6 & w^9\\ \end{array}\right] \left[ \begin{array}{ll} x_0\\ x_1\\ x_2\\ x_3\\ \end{array}\right]$
把矩阵乘积算出来，但是偶数项写在前面，并且提取公因式，还要切记 $w^0=1$ ，有时候不要迷糊。
$z_0 = w^0x_0+w^0x_2+w^0(w^0x_1+w^0x_3)\\ z_1 = w^0x_0+w^2x_2+w^1(w^0x_1+w^2x_3)\\ z_2 = w^0x_0+w^4x_2+w^2(w^0x_1+w^4x_3)\\ z_3 = w^0x_0+w^6x_2+w^3(w^0x_1+w^6x_3)\\$
由于 $w^4=1=w^0$ ，可以得到
$z_0 = w^0x_0+w^0x_2+w^0(w^0x_1+w^0x_3)\\ z_1 = w^0x_0+w^2x_2+w^1(w^0x_1+w^2x_3)\\ z_2 = w^0x_0+w^0x_2+w^2(w^0x_1+w^0x_3)\\ z_3 = w^0x_0+w^2x_2+w^3(w^0x_1+w^2x_3)\\$
后面部分，冈萨雷斯和《数值分析》又不一样了，两个方法由差异，以冈萨雷斯为准好像又更好懂一些。
即，令 $2n =4$ 以上部分又可以写成
$z_k = \sum_{j=0}^{2n-1}x_j w^{jk} = \sum_{j=0}^{n-1}x_{2j} w^{k(2j)} + \sum_{j=0}^{n-1}x_{2j+1} w^{k(2j+1)}$
如此迭代分治下去，就是cpu版本的快速傅里叶变换了。

五、现有代码实验

5.1 使用单位阵测试

输入矩阵为 $I_{4 \times 4}$

5.1.1 MATLAB

Matlab代码如下，其中，fft是一维傅里叶变换，fft2是二维傅里叶变换。

A = eye(4,4)
A_fft = fft(A)
A_fft2 = fft2(A)

得到结果如下：

A =

     1     0     0     0
     0     1     0     0
     0     0     1     0
     0     0     0     1

A_fft =

   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i
   1.0000 + 0.0000i   0.0000 - 1.0000i  -1.0000 + 0.0000i   0.0000 + 1.0000i
   1.0000 + 0.0000i  -1.0000 + 0.0000i   1.0000 + 0.0000i  -1.0000 + 0.0000i
   1.0000 + 0.0000i   0.0000 + 1.0000i  -1.0000 + 0.0000i   0.0000 - 1.0000i

A_fft2 =

     4     0     0     0
     0     0     0     4
     0     0     4     0
     0     4     0     0

5.1.2 opencv的cpu版本

代码如下，放到函数里面就行，

{
   Mat src = Mat::eye(Size(4,4),CV_32FC1);
   Mat FFT;
   dft(src , FFT, 0/*DFT_REAL_OUTPUT*/);
}

结果

捕获fft.PNG

结果是这样的原因是opencv输入实数矩阵的时候的输出时经过特殊编码的，编码方式如下，文档说明中说是按ipp的格式来编码的。可以参考opencv dft文档。

捕获cpufftcode.PNG

5.1.3 opencv的gpu版本

经过查看cv::cuda::dft的源码了解到， opencv的cuda加速的版本的fft就是直接调用了Nvdia的cuFFT接口实现的。可以参考opencv cv:cuda::dft文档和cuFFT输入输出格式。
主要实现语句诸如以下几句

cufftPlan1d(&plan, dft_size_opt.width, dft_type, dft_size_opt.height);
cufftPlan2d(&plan, dft_size_opt.height, dft_size_opt.width, dft_type);

cufftExecC2C(plan, src_cont.ptr<cufftComplex>(), dst.ptr<cufftComplex>(),
                            is_inverse ? CUFFT_INVERSE : CUFFT_FORWARD);
cufftExecC2R(plan, src_cont.ptr<cufftComplex>(), dst.ptr<cufftReal>());
cufftExecR2C(plan, src_cont.ptr<cufftReal>(), dst.ptr<cufftComplex>());

调用也很简单

{
  Mat src = Mat::eye(Size(4,4),CV_32FC1);
  Mat FFT;
  cuda::GpuMat gP,  gFFT;
  gP.upload(src);
  cuda::dft(gP, gFFT, padded2.size(), 0);
  gFFT.download(FFT);
}

得到的结果如下，列数的解释在cuFFT输入输出格式。

捕获cufft.PNG

5.2 使用一般矩阵测试

输入矩阵

$A= \left[ \begin{array}{ll} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 1 & 0 & 0 & 1\\ \end{array}\right]$

5.2.1 Matalb结果

A_fft =

   2.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i
   1.0000 + 1.0000i   0.0000 - 1.0000i  -1.0000 + 0.0000i   0.0000 + 1.0000i
   0.0000 + 0.0000i  -1.0000 + 0.0000i   1.0000 + 0.0000i  -1.0000 + 0.0000i
   1.0000 - 1.0000i   0.0000 + 1.0000i  -1.0000 + 0.0000i   0.0000 - 1.0000i
A_fft2 =

   5.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i
   0.0000 + 1.0000i   0.0000 + 1.0000i   0.0000 + 1.0000i   4.0000 + 1.0000i
  -1.0000 + 0.0000i  -1.0000 + 0.0000i   3.0000 + 0.0000i  -1.0000 + 0.0000i
   0.0000 - 1.0000i   4.0000 - 1.0000i   0.0000 - 1.0000i   0.0000 - 1.0000i

5.2.2 opencv cpu版本

捕获fft1.PNG

5.2.3 opencv gpu版本

捕获cufft1.PNG

六、cpu实现

应该会有一个别人的和一个自己改的的版本，至于opencv的版本写得过于复杂，引用依赖过多，就不拷贝过来做参考了。

6.1 蝴蝶操作

搁置

6.2

搁置

七、gpu实现

7.1 cuFFT头文件

首先说一下，据说cuda提供了cuFFT头文件，可以测试下。也可以不用测试，毕竟opencv的gpu版本已经就是了。另外也有很好的参考资料 CUDA学习笔记3：CUFFT（CUDA提供了封装好的CUFFT库）的使用例子。

7.2 自行实现

搁置
可参考CUDA实现FFT快速傅里叶变换;

八、结论

各家提供的的结果编码方式差异好大，好坑。