傅里叶快速算法（FFT）的理解与实现

最新推荐文章于 2024-05-11 10:00:52 发布

大瓜子V5

最新推荐文章于 2024-05-11 10:00:52 发布

阅读量8.7k

点赞数 11

分类专栏：数字图像处理算法分析与设计文章标签：算法图像处理快速傅立叶变换 python

本文链接：https://blog.csdn.net/weixin_43959755/article/details/114278150

版权

数字图像处理同时被 2 个专栏收录

3 篇文章 9 订阅

订阅专栏

算法分析与设计

2 篇文章 0 订阅

订阅专栏

上篇文章中我介绍了离散傅里叶变换（参考这里），所以这次决定在DFT的基础上再实现一下FFT算法。

一、快速傅里叶变换的意义

直接实现DFT的公式在频率域的处理中并不实际。蛮力实现需要约 $(MN)^{2}$ 次乘法和加法运算。对于中等大小的图像（如2048x2048），进行一次二维DFT运算就需要约17万亿次乘法和加法，且还不包括计算一次并存储到查找表中的指数运算。所以如果说人们没有发现能够将计算量降到 $MNlog_{2}MN$ 次乘法和加法的快速傅里叶变换（FFT），那么之前介绍的内容也将没有任何实用价值。

可以说，快速傅里叶变换算法的“发现”，引发了信号处理领域的变革。考虑大小为 $M\times M$ 和 $m\times m$ 的方形图像与核，相比于采用不可分离的核，采用FFT对空间域滤波的计算优势可以定义为

$C_{n}(m)=\frac{M^{2}m^{2}}{2M^{2}log_{2}M^{2}}=\frac{m^{2}}{4log_{2}M}$

如果核是可分离的，那么这一优势变为

$C_{s}(m)=\frac{2M^{2}m}{2M^{2}log_{2}M^{2}}=\frac{m}{2log_{2}M}$

两种情况下，当 $C(m)>1$ 时，FFT方法的计算优势更大（根据计算次数）；而当 $C(m)\leqslant 1$ 时，空间滤波的优势更大。可以说，FFT相较于较大的空间核具有压倒性的优势。

二、算法原理

我们用于实现这一目标的算法是逐次加倍法，这个算法是导致整个产业诞生的最初算法。要求假设的取样数是2的整数次幂，但并不是其他方法的通用要求。推导FFT的完整过程如下：

从最后几个式子中，我们可以惊奇地发现，一个M点DFT可通过把原表达式分成两部分来计算。 $F(u)$ 的前一半和后一半都可以由 $F_{even}(u)$ 和 $F_{odd}(u)$ 直接代入得到，不需要额外的变换计算。而且 $F_{even}(u)$ 和 $F_{odd}(u)$ 本身也可以再分成两个子序列，从而形成一种递归结构。最后这几步，我们用伪代码的形式表现为

在这个循环中，把旋转因子 $\omega = \omega_{n}^{k}$ 乘以 $y_{k}^{[1]}$ ，把所得乘积存入 $t$ 中，然后从 $y_{k}^{[0]}$ 中增加及减去 $t$ ，这一系列操作成为一个蝴蝶操作（bufferfly operation），下图展示了蝴蝶操作的执行步骤。

图4 (a)两个输入向量从左边进入，和与差在右边输出。
(b)一个蝴蝶操作的简化草图。

因为DFT的实际应用中需要尽可能快的速度，所以本次实现采用迭代结构而不是递归结构，这就要求我们将输入向量安排在一次FFT调用相关的各次递归调用中，将输入向量安排成树形结构。如下图所示

观察此树注意到，如果将输入向量按照树中叶的次序进行安排，就可以自底向上对FFT过程进行追踪。每次我们都成对取出元素，利用一次蝴蝶操作计算出每对的DFT，然后用其DFT取代这对元素。之后我们重复这一过程向上进行，直至向量包含两个具有n/2个元素的DFT，这时再综合应用n/2次蝴蝶操作，就可以合成最终具有n个元素的DFT。

为了让输入向量a中的元素符合要求，我们还需要对一开始的元素进行一次位逆序置换。方法是将元素一开始的下标二进制表示以后，将其各位逆序再转化成十进制。例如一开始的二进制序列[000,001,010,011,100,101,110,111]逆序表示后为[000,100,010,110,001,101,011,111]，即十进制[0,4,2,6,1,5,3,7]符合我们所期望的排序。

三、代码实现

接下来展示位逆序置换及一维FFT迭代算法的代码。

#位逆序置换
def bit_reverse(l):
    #初始大小
    M = l.shape[0]
    bits = int(np.log2(M))
    
    #创建新的列表用于保存结果
    L = np.zeros(l.shape, dtype=np.complex128)
    
    #遍历数组进行位逆序
    for i in range(M):
        b = int(bin(i)[2:].zfill(bits)[::-1], 2)
        L[b] = l[i]
        
    return L

#快速傅里叶变换
def fft(f):
    #位逆序处理
    F = bit_reverse(f)
    
    #得到长度
    n = f.shape[0]
    
    #log2(n)次奇偶序列合并
    for s in range(1, int(np.log2(n))+1):
        #m为子序列长度
        m = 2**s
        base = np.exp(-1j*2*np.pi/m)
        #遍历每个奇偶序列
        for k in range(0, n-1, m):
            #w初始化
            w = 1
            #蝴蝶操作
            for j in range(int(m/2)):
                t = w*F[k+j+int(m/2)]
                u = F[k+j]
                F[k+j] = u+t
                F[k+j+int(m/2)] = u-t
                w = w*base
    
    return F

然后二维和反变换算法的实现方法与普通DFT类似，代码如下。中心化和展示频谱的方法是通用的，上篇文章中已经介绍过，便不再赘述，都可以点击这里了解。

#二维快速傅里叶变换
def fft_2d(f):
    #得到输入的行数和列数
    M, N = f.shape[0], f.shape[1]
    
    #初始化两个复数类型数组，用于保存结果
    F, F_x = np.zeros(shape=(M,N), dtype=np.complex128), np.zeros(shape=(M,N), dtype=np.complex128)
    
    #逐行逐列进行两次一维傅里叶变换
    for i in range(M):F_x[i,:] = fft(f[i,:])
    for i in range(N):F[:,i] = fft(F_x[:,i])
    
    return F

#快速离散傅里叶反变换
def ifft_2d(F):
    #得到行数和列数
    M, N = F.shape[0], F.shape[1]
    
    #先对F的共轭进行正向离散傅里叶变换，除以常数后再取共轭
    f = np.conjugate(fft_2d(np.conjugate(F)))/M/N
    
    return f

四、效果展示

和上次一样，先展示频谱图。

图6 (a)原图；(b)使用自己完成的FFT函数变换得到频谱图；
(c)使用自己完成的FFT函数和中心化函数变换得到的频谱图；
(d)使用numpy包中FFT函数和中心化函数变换得到的频谱图。

然后是与numpy包函数混用的情况。

图7 (a)原图；
(b)使用自己实现的函数进行正变换后使用numpy包中的函数进行反变换；
(c)使用numpy包中的函数进行正变换后使用自己实现的函数进行反变换。

为了展示FFT速度上的优势，我还使用了这几个方法对512x512像素大小的图片进行了处理，并打印它们的运行的时间进行对比。可以看出，使用自己实现的FFT算法确实也比普通DFT算法快不少，但是较numpy包的实现仍有不小差距（考虑过可能是python不太适合迭代运算的原因，下次可以考虑使用c++实现一下），但总的来说实验还算成功。

图8 各个方法实现DFT的运行时间，从上到下依次是：
自己实现的DFT普通算法、自己实现的FFT算法、numpy包自带的FFT算法

五、参考资料

[1] Rafael C.Gonzalez, Richard E.Woods. 数字图像处理(第四版)[M]. 阮秋琦, 阮宇智, 译. 北京:电子工业出版社, 2020:210-212

[2] (美)Thomas H.Cormen等. 算法导论[M]. 殷建平, 徐云等, 译. 北京:机械工业出版社, 2013:527-539

大瓜子V5

关注

11
点赞
踩
42

收藏

觉得还不错? 一键收藏
0
评论
傅里叶快速算法（FFT）的理解与实现

上篇文章中我介绍了离散傅里叶变换（参考这里），所以这次决定在DFT的基础上再实现一下FFT算法。一、快速傅里叶变换的意义直接实现DFT的公式在频率域的处理中并不实际。蛮力实现需要约次乘法和加法运算。对于中等大小的图像（如2048x2048），进行一次二维DFT运算就需要约17万亿次乘法和加法，且还不包括计算一次并存储到查找表中的指数运算。所以如果说人们没有发现能够将计算量降到次乘法和加法的快速傅里叶变换（FFT），那么之前介绍的内容也将没有任何实用价值。图1 冈萨雷斯书中FFT关于一维DFT的
复制链接

扫一扫