我对快速离散傅里叶变换(FFT)的理解

最新推荐文章于 2024-06-25 11:15:35 发布

置顶 cracy3m

最新推荐文章于 2024-06-25 11:15:35 发布

阅读量1.9k

点赞数

分类专栏：光学相干层析成像

本文链接：https://blog.csdn.net/cracy3m/article/details/109098258

版权

光学相干层析成像专栏收录该内容

8 篇文章 9 订阅

订阅专栏

我对快速离散傅里叶变换(FFT)的理解

记录一下自己对快速离散傅里叶变换 (FFT) 的理解

离散傅里叶变换(DFT)是假设个数为N个的任意离散数据 $x_{0},x_{1},x_{2}...,x_{N-1}$

可以变换为 $x_{n}=\sum_{k=0}^{N-1}A_{k}e^{j\frac{2\pi kn}{N}}$ ,

其中 $A_{k}$ 为对应原始数据k倍频率分量的复数常量(幅值和相位)， $e^{j\frac{2\pi kn}{N}}$ 是随着n变化 k倍频率的复数周期函数，

求 $A_{k}$ 的公式为 $A_{k}=\frac{1}{N} \sum_{n=0}^{N-1}x_{n} * e^{-j\frac{2\pi k n}{N}}$ k=0,1,2,3....N-1

如果将 $e^{-j\frac{2\pi k n}{N}}$ 以 $W_{N}^{kn}$ 标记来代替， $N*A_{k}$ 用 X(k) 标记来代替，公式可写成 $X(k)=\sum_{n=0}^{N-1}x_{n} * W_{N}^{kn}$ ，

$W_{N}^{kn}$ 的一些特性

$W_{N}^{k(n+N)}=e^{-j\frac{2\pi k n}{N}}*e^{-j\frac{2\pi k N}{N}}=e^{-j\frac{2\pi k n}{N}}*e^{-j2\pi k }=e^{-j\frac{2\pi k n}{N}}=W_{N}^{kn}$ $=W_{N}^{(k+N)n}$

$W_{N}^{k(n+x)}=W_{N}^{kn}*W_{N}^{kx}$ $W_{N}^{(k+x)n}=W_{N}^{kn}*W_{N}^{nx}$

$W_{N}^{kn*x}=W_{N/x}^{kn}$

$W_{N}^{kn/x}=W_{N*x}^{kn}$

$W_{N}^{k(n+N/2)}=W_{N}^{kn}*W_{2}^{K}$

$W^{n}=1$ $W^{N/2}_{N}=W_{2}=-1$

将公式 $X(k)=\sum_{n=0}^{N-1}x_{n} * W_{N}^{kn}$ 按n值分别为偶数奇数分成两部分如下一步步化简

$X(k)=\sum_{r=0}^{\frac{N}{2}-1}x_{2r} * W_{N}^{k2r} + \sum_{r=0}^{\frac{N}{2}-1}x_{2r+1} * W_{N}^{k(2r+1)}$

$=\sum_{r=0}^{\frac{N}{2}-1}x_{2r} * W_{N/2}^{kr} + \sum_{r=0}^{\frac{N}{2}-1}x_{2r+1} * W_{N}^{k2r}* W_{N}^{k}$

$=\sum_{r=0}^{\frac{N}{2}-1}x_{2r} * W_{N/2}^{kr} + W_{N}^{k}\sum_{r=0}^{\frac{N}{2}-1}x_{2r+1} * W_{N/2}^{kr}$

可以看到 $\sum_{r=0}^{\frac{N}{2}-1}x_{2r} * W_{N/2}^{kr}$ $\sum_{r=0}^{\frac{N}{2}-1}x_{2r+1} * W_{N/2}^{kr}$ 与 $\sum_{n=0}^{N-1}x_{n} * W_{N}^{kn}$ 形式上是相同的

$x_{2r}$ 是 $x_{0},x_{2},x_{4}...,x_{N-2}$ $x_{2r+1}$ 是 $x_{1},x_{3},x_{5}...,x_{N-1}$

当 $k=0,1,2,...\frac{N}{2}-1$ 时， $\sum_{r=0}^{\frac{N}{2}-1}x_{2r} * W_{N/2}^{kr}$ $\sum_{r=0}^{\frac{N}{2}-1}x_{2r+1} * W_{N/2}^{kr}$ 实际上是在原始数据 $x_{0},x_{1},x_{2}...,x_{N-1}$ 抽取偶数和奇数序列的离散傅里叶变换分别标记为 $X^{-1}_{e}(k)$ 和 $X^{-1}_{o}(k)$ (-1是指上标，并不是指求幂)

$X^{-1}_{e}(k)=\sum_{r=0}^{\frac{N}{2}-1}x_{2r} * W_{N/2}^{kr}$ $X^{-1}_{o}(k)=\sum_{r=0}^{\frac{N}{2}-1}x_{2r+1} * W_{N/2}^{kr}$

所以 $X(k)=\sum_{r=0}^{\frac{N}{2}-1}x_{2r} * W_{N/2}^{kr} + W_{N}^{k}\sum_{r=0}^{\frac{N}{2}-1}x_{2r+1} * W_{N/2}^{kr}=X^{-1}_{e}(k)+W_{N}^{k}*X^{-1}_{o}(k)$

那么当 $k=\frac{N}{2},\frac{N}{2}+1,...N-1$ 时呢？？？

记 $k'=0,1,2,...\frac{N}{2}-1$ $X(k)=\sum_{r=0}^{\frac{N}{2}-1}x_{2r} * W_{N/2}^{kr} + W_{N}^{k}\sum_{r=0}^{\frac{N}{2}-1}x_{2r+1} * W_{N/2}^{kr}$ 可以转换成

$X(k)=\sum_{r=0}^{\frac{N}{2}-1}x_{2r} * W_{N/2}^{(k'+N/2)r} + W_{N}^{(k'+N/2)}\sum_{r=0}^{\frac{N}{2}-1}x_{2r+1} * W_{N/2}^{(k'+N/2)r}$

$=\sum_{r=0}^{\frac{N}{2}-1}x_{2r} * W_{N/2}^{k'r}*W^{rN/2}_{N/2} + W_{N}^{k'}*W^{N/2}_{N}*\sum_{r=0}^{\frac{N}{2}-1}x_{2r+1} * W_{N/2}^{k'r}*W^{rN/2}_{N/2}$

$=\sum_{r=0}^{\frac{N}{2}-1}x_{2r} * W_{N/2}^{k'r}*W^{r} - W_{N}^{k'}*\sum_{r=0}^{\frac{N}{2}-1}x_{2r+1} * W_{N/2}^{k'r}*W^{r}$

$=\sum_{r=0}^{\frac{N}{2}-1}x_{2r} * W_{N/2}^{k'r} - W_{N}^{k'}*\sum_{r=0}^{\frac{N}{2}-1}x_{2r+1} * W_{N/2}^{k'r}$

$=X^{-1}_{e}(k')-W_{N}^{k'}*X^{-1}_{o}(k')$

也就是说当 $k=\frac{N}{2},\frac{N}{2}+1,...N-1$ $k'=0,1,2,...\frac{N}{2}-1$ $X(k)=X^{-1}_{e}(k')-W_{N}^{k'}*X^{-1}_{o}(k')$

小总结

$X(k)=\sum_{n=0}^{N-1}x_{n} * W_{N}^{kn}$ k=0,1,2,3....N-1 可以分成 $k=0,1,2,...\frac{N}{2}-1$ 和 $k=\frac{N}{2},\frac{N}{2}+1,...N-1$ 两部分计算

--------------------------------------------------------------------------------------------------------------------------------------------

当 $k=0,1,2,...\frac{N}{2}-1$ 时 $X(k)=X^{-1}_{e}(k)+W_{N}^{k}*X^{-1}_{o}(k)$

当 $k=\frac{N}{2},\frac{N}{2}+1,...N-1$ $k'=0,1,2,...\frac{N}{2}-1$ $X(k)=X^{-1}_{e}(k')-W_{N}^{k'}*X^{-1}_{o}(k')$

--------------------------------------------------------------------------------------------------------------------------------------------

又或者 $k=0,1,2,...\frac{N}{2}-1$ $X(k)=X^{-1}_{e}(k)+W_{N}^{k}*X^{-1}_{o}(k)$

$X(k+\frac{N}{2})=X^{-1}_{e}(k)-W_{N}^{k}*X^{-1}_{o}(k)$

--------------------------------------------------------------------------------------------------------------------------------------------

$X^{-1}_{e}(k)$ 是抽取的偶序列数据 $x_{0},x_{2},x_{4}...,x_{N-2}$ 的傅里叶变换

$X^{-1}_{o}(k)$ 是抽取的奇序列数据 $x_{1},x_{3},x_{5}...,x_{N-1}$ 的傅里叶变换

--------------------------------------------------------------------------------------------------------------------------------------

$k=0,1,2,...\frac{N}{2}-1$ 可简化成下面的图形，（因此叫做蝶形算法吧）

对FFT 原始数据序列重选择的规律总结

假设有原始数据 x0,x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15

那么它的傅里叶变换 X(0...15) 一级简化如下图所示，由于2的4次方为16，所以可以用4位二进制位表示原始数据的索引号

下图中将原始数据分成奇偶两组即索引号分为 4'bxxx0 和 4'bxxx1 xxx为从 000~111二进制数

二级简化如下图所示索引号从上一次为 4'bxxx0 和 4'bxxx1 的两组再依次分别分成奇偶两组

即 4'bxxx0 -> 4'bxx00, 4'bxx10 4'bxxx1 -> 4'bxx01, 4'bxx11

三级简化如下图所示

索引号从上一次为 4'bxx00, 4'bxx10, 4'bxx01, 4'bxx11 的四组再依次分别分成奇偶两组

即 4'bxx00 -> 4'bx000, 4'bx100 4'bxx10-> 4'bx010, 4'bx110 4'bxx01 -> 4'bx001, 4'bx101 4'bxx11-> 4'bx011, 4'bx111

四级简化如下图由于上一级输入端本来就剩余两个数据，也就是一偶一奇不能再分组，所以索引号与上面一级相同

总结出的原始索引号变化规律如下图所示 xxx代表自增序列

观察最后一次化简的原始数据索引号，实际上是 4'b0000 ~ 4'b11111 自增序列的镜像如下表所示

这种规律同样也适用于多于或小于16个数据的情况，但二进制索引号最大位数要作相应的变化

python 实现FFT示例代码

# coding:utf-8

import numpy as np
from six.moves import xrange
import cmath

# 计算 fft 的镜像二进制索引
mirrorindex8bit = [0] * 256
for i in xrange(256):
    bitstr = ""
    for j in xrange(8):
        bit = (i & (0x01 << j)) != 0
        if bit:
            mirrorindex8bit[i] = mirrorindex8bit[i] | 0x01 << (7 - j)
            bitstr = bitstr + "1"
        else:
            bitstr = bitstr + "0"
    # print(mirrorindex8bit[i])
    # print(i)
    # print(strbit)

mirrorindex16bit = [0] * 65536
for i in xrange(65536):
    mirrorindex16bit[i] = mirrorindex8bit[(i & 0x00FF00) >> 8] | (mirrorindex8bit[i & 0x00FF] << 8)
    # print(mirrorindex16bit[i])

origindata = [1, 2, 5, 9, -5, 8, 7, -1, -3]

# 求取填充长度
padlen = 0
padpowerbit = 0
for i in xrange(16):
    padlen = 2 ** (i + 1)
    if (padlen >= len(origindata)):
        padpowerbit = i + 1
        break
    if (i == (16 - 1)):
        print("data len exceed 65536")
        exit()

# 原始数据 第一次蝶形运算时的 索引号
index = [0] * padlen
for i in xrange(padlen):
    index[i] = mirrorindex16bit[i << (16 - padpowerbit)]
    # print(i, index[i])

# 计算 W 旋转系数
Wk = [0 + 0j] * (padlen // 2)
for i in xrange(len(Wk)):
    Wk[i] = cmath.exp(-1j * 2 * cmath.pi * i / padlen)

# 定义两个复数缓冲区 用于放置变换时的源和目标
buf0 = [0 + 0j] * padlen
buf1 = [0 + 0j] * padlen

src_buf = buf0
dst_buf = buf1

# 将原始数据 放置到源缓冲区
for i in xrange(len(src_buf)):
    if index[i] < len(origindata):
        src_buf[i] = complex(origindata[index[i]])

# 蝶形运算
calctime = padpowerbit  # 总共要运算的级数
groupcnt = padlen // 2  # 每次的分组数
itemingrp = padlen // groupcnt  # 每次的每组里的项数

for level in xrange(calctime):
    for grpi in xrange(groupcnt):
        for i in xrange(itemingrp // 2):  # 每组奇偶一起运算
            # print(grpi * itemingrp + i, grpi * itemingrp + itemingrp // 2 + i)
            odditem = src_buf[grpi * itemingrp + itemingrp // 2 + i] * Wk[i * groupcnt]
            dst_buf[grpi * itemingrp + i] = src_buf[grpi * itemingrp + i] + odditem
            dst_buf[grpi * itemingrp + itemingrp // 2 + i] = src_buf[grpi * itemingrp + i] - odditem
    if level != calctime - 1:  # 不是最后一次都要进行更新
        groupcnt = groupcnt // 2  # 更新组数
        itemingrp = padlen // groupcnt  # 更新组里面的项数
        if level % 2:  # 交换源和目的缓冲区
            src_buf = buf0
            dst_buf = buf1
        else:
            src_buf = buf1
            dst_buf = buf0

# FFT 结果
print("FFT result:")
for i in xrange(len(dst_buf)):
    print("X(%d)=%s" % (i, dst_buf[i].__str__()))

# 进行FFT-1 逆变换
# 将数据 放置到源缓冲区
for i in xrange(len(src_buf)):
    src_buf[i] = dst_buf[index[i]] / padlen

# 计算 W 旋转系数共轭值
for i in xrange(len(Wk)):
    Wk[i] = Wk[i].conjugate()

# 蝶形运算
calctime = padpowerbit  # 总共要运算的级数
groupcnt = padlen // 2  # 每次的分组数
itemingrp = padlen // groupcnt  # 每次的每组里的项数
for level in xrange(calctime):
    for grpi in xrange(groupcnt):
        for i in xrange(itemingrp // 2):  # 每组奇偶一起运算
            # print(grpi * itemingrp + i, grpi * itemingrp + itemingrp // 2 + i)
            odditem = src_buf[grpi * itemingrp + itemingrp // 2 + i] * Wk[i * groupcnt]
            dst_buf[grpi * itemingrp + i] = src_buf[grpi * itemingrp + i] + odditem
            dst_buf[grpi * itemingrp + itemingrp // 2 + i] = src_buf[grpi * itemingrp + i] - odditem
    if level != calctime - 1:  # 不是最后一次都要进行更新
        groupcnt = groupcnt // 2  # 更新组数
        itemingrp = padlen // groupcnt  # 更新组里面的项数
        t_buf = src_buf
        src_buf = dst_buf
        dst_buf = t_buf

# FFT 逆变换结果
print("invert FFT result:")
for i in xrange(len(dst_buf)):
    print("x%d=%s" % (i, dst_buf[i].__str__()))