傅里叶变换公式表_详尽的快速傅里叶变换推导与Unity中的实现

最新推荐文章于 2020-12-11 18:27:08 发布

weixin_39636164

最新推荐文章于 2020-12-11 18:27:08 发布

阅读量913

点赞数

文章标签：傅里叶变换公式表傅里叶描述子欧氏距离快速傅里叶变换fft 快速傅里叶变换matlab代码

图形学中偶尔会遇到傅里叶变换，我对这个一直处于半懂不懂的状态。网上找了不少实现，每种实现都不太一样，但是又不懂为什么不一样也能正确。为了让自己彻底理解FFT，因此花了一些时间自己进行推导。

本文还有下篇，关于更快速的快速傅里叶变换的实现和性能分析：

https://zhuanlan.zhihu.com/p/211268502zhuanlan.zhihu.com

不想看文字和公式，直接想跑代码，看这里：

https://github.com/MioMioReimu/RadixNFFTgithub.com

一、连续傅里叶变换和离散傅里叶变换

假设

是复平面上的勒贝格可积的函数，则定义

的连续傅里叶变换

是对于任意实数

在这里，我们解释

为角频率，

是复数，并且是信号在该频率成分处的幅度和相位。

傅里叶逆变换是

怎么理解傅里叶变换这2个公式呢。傅里叶提出，任何连续周期函数，都可以由不同频率不同相位不同幅度的正弦波组合而成。上述的傅里叶变换公式表达的是，如何求给定的频率，它在信号

中的幅度和相位。

在计算机中，我们无法处理连续的积分行为，因此只能进行离散化的采样了。对于原始信号

的离散采样

, 我们定义，它的离散傅里叶变换（DFT）为

离散傅里叶逆变换（IDFT）为

这里k表达的是离散的频率，n是离散的信号采样点，

。

二、离散傅里叶变换的快速算法

快速傅里叶变换利用的是分治法的思想。

2.1 奇偶分割

我们在推导复杂的公式之前，先给出一张分治的原理示意图。

Fig1. DFT8 分治为2个 DFT4

图1中，DFT8按照奇偶分解为2个DFT4，2个DFT4先算好，然后再利用蝶形变换简单处理就可以得到DFT8的结果。这就是FFT的核心思想。

下面，我们用公式来推导这个最核心的蝶形变换。

假设

, 我们将

按n的奇偶分为2组：令

，则

我们可以进行一些函数代换，假设

和

分别是

的傅里叶变换, 其中

。则上式等于

上述式子中，我们成功的将

转换成了2个子傅里叶变换之和，我们可以了解到上述式中k的范围在

中，因为两个子傅里叶变换的k的区间只有这么大。

我们求出了离散傅里叶变换分解式

，但是k的区间还不够，还差一半。

2.2 求下一半的公式

下面将描述如何求取

。

由

由欧拉方程

，另外

属于正整数，可知

因此

同理，

, 因此

离散傅里叶变换的另一半的分解式

注意到，这里公式的2种方式都是可以滴，有的代码实现就是用第二个子式实现的。

2.3 蝶形运算

为了简化下面的一些符号的表示，我们令

我们定义，蝶形运算公式

它的图示为

Fig2. 单个蝶形变换

上图中，

和

分别是

的奇数部分和偶数部分的子傅里叶变换的结果，他们的大小为

。

2.4 两个DFT4利用蝶形变换生成DFT8

如下图所示，我们利用2个黑盒生成的DFT4生成DFT8. 首先，我们先取DFT8的原数据的奇数和偶数的信号，然后用某种黑盒过程生成DFT4，然后再利用蝶形变换生成DFT8.

Fig3. DFT8 to DFT4

2.5 利用4个DFT2生成2个DFT4，再生成DFT8

上述说的黑盒过程，可以进一步利用蝶形变换分解。下图的绿色部分标记出了新增的蝶形变换部分。注意到，上面的序号发生了变化，每个4个DFT的奇数取出来。

Fig4. DFT8 to DFT4 to DFT2

2.6 将DFT2直接转换为蝶形变换

通过上述三个递归步骤，可以了解到FFT的详细做法了。利用分治思想，一直2分，直到只有N=2时，即可直接变为蝶形变换。

Fig5. DFT8 to DFT4 to DFT2 to DFT1

2.7 序列输入顺序

我们可以注意到，最原始的DFT8在完全分解成蝶形变换时，

，序列发生了变化，因为除了N=2之外，每一次分解为蝶形变换时，因为取奇偶导致顺序发生了错位。这种看起来没啥规律，但是很明显，这种规律交换的蝶形运算，最后这个序号肯定有规律的。规律刚好是二进制的位反 (bitreverse)。如下图所示

三、FFT具体实现

我们注意到，上述过程 FFT的具体蝶形如何做，有很多种方案，下面这个链接的PDF里总结了好多种。

Fast Fourier Transform

下面我着重讲述2个最常见的类型，第一个是最原版的先换顺序然后进行蝶形变换。第二个是GPU中最常见的FFT实现。首先，为了简化细节，我直接将第二章描述的蝶形变换内部细节标记为一个黑盒，如下图。

3.1 经典快速傅里叶变换算法(Cooley-Tukey FFT)

经典FFT算法，就是先进行一次序号的位反计算出新的序号，然后按照鲽形图进行递归或者非递归计算。下图是Cooley-Tukey FFT的DFT8的详细流程图

Fig6. Cooley-Tukey FFT N=8

我们可以注意到，在Cooley-Tukey FFT的实现中，除了一开始需要映射表的序号之外，所有的蝶形变换中，输入的两个元素的序号和输出的2个元素的序号没有发生变化。例如DFT2->DFT4中，4和6元素的下标分别是1, 3, 他们进行蝶形变换后输出的下标也是1, 3。

3.2 适合GPU的FFT(Stockham FFT)

上述的FFT，输入数据要进行一次序号重映射行为，在GPU端，这种映射可通过查表实现。但是，有一种Stockham FFT，它提出了一种巧妙的映射方式，使得id能直接计算映射，这在GPU中相比查表会更加快捷，对于大型FFT来说，构建表传输表都比较慢。下图来自其他地方，这两个图简单明了的描述了区别。

Cooley-Tukey vs Stockham

下图是Stockham FFT的详细流程图。

Stockham FFT N=8

下面，我们解析一下Stockham FFT的流程图我们从第0轮开始分析，首先，数据被分成8组，每个数据单独成组。然后，(0,4)(1,5)(2,6)(3,7)分别是DFT2的奇偶子项，他们分别做蝶形变换，完成了第0轮。第一轮中，数据被分成了4组，

, 其中

与

分别是DFT4的奇偶子项，然后他们分别做蝶形变换，完成第一轮。以此类推。。。

下面，我们解析以下如何计算Stockham FFT的映射。

3.2.1 Stockham蝶形变换的输入下标的求取

我们观察可以发现以下特征：

所有蝶形变换的2个输入对的下标都符合
的形式，

是数据的长度。例如上面的图里N=8.
所有蝶形变换的第一个输入的下标组合起来刚好是
的形式。

根据这2点，我们可以在代码里确定蝶形变换的输入的下标了。

假设某个蝶形变换左边输入的下标为
，则第二个输入的下标是

, 其中

。

注意到每一轮FFT合并过程中，都执行了

次蝶形变换，因此我们可以定义：

假设某个蝶形变换左边输入的下标为
，则称这个蝶形变换为第

个蝶形变换。

然后我们来求，蝶形变换的输出的下标。

3.2.2 Stockham蝶形变换的输出下标的求取

可以观察到：

每个蝶形变换的2个输出的下标间隔会随着每一轮蝶形变换而呈指数扩张。例如 DFT1 -> DFT2中，2个输出的下标间隔是1. DFT2->DFT4中，2个输出下标的间隔是2。 DFT4->DFT8中，2个输出下标间隔是4.
每一轮内，每组DFT内的蝶形变换，都是左边输出排在右边输出前面。然后，每一轮内不同组之间延后摆放。

因此，我们可以求出以下结论：

1. 第m轮(m从0开始)，该轮内的子DFT的大小是
, 输出的DFT的大小是

。

2. 第m轮(m从0开始)，第
个蝶形变换，它在某个子DFT内的下标是

。

3. 第m轮(m从0开始)，第
个蝶形变换的左边输出的下标为

4. 右边输出的下标为

3.2.3 Stockham根据输出下标求取输入的下标

在我浏览的一些实现中，GPU是通过输出的下标进行分发任务的，也就是说，每个GPU函数中输出的下标确定了，这时候我们需要反向求出蝶形变换的2个输入下标，以及该输出在蝶形变换中的正负。我们通过计算和分析可得出

假设第m轮(m从0开始)，第x个输出下标所属的
1. 蝶形变换的序号是

2. 根据蝶形变换下标的定义可知左边的输入的下标就是蝶形变换的序号
3. 右边的输入的下标是

4. 输出的正负性是
为正，否则为负。(一般这种实现，不需要利用正负性，直接利用公式

的第二种形式即可)

四、Unity中实现Radix2的Stockham FFT

4.1 Compute Shader简介

Unity中我使用了Compute Shader来实现FFT。Compute Shader跟CUDA编程的概念基本一一对应。

首先，需要暴露给Host端的函数，需要用修饰为kernel函数。kernel函数都需要指明一个WorkGroup的大小。

然后Host端设置参数和数据，最后Dispatch N个并行任务。

如下面的代码所示，如果dit_x_radix2函数Dispatch 4x4x1次，则这个函数会被并行的执行4x4x1次WorkGroup，每个WorkGroup的都是8x8x1的大小。 dit_x_radix2函数会被执行 (4x8)x(4x8)x(1x1)次。

#define NTHREADS 8
[numthreads(NTHREADS, NTHREADS, 1)]
void dit_x_radix2(uint3 id: SV_DispatchThreadID) {
    // code...
}

更多关于Compute Shader的详细入门教程，知乎也有不少，这里先不介绍了。

4.2 Stockham FFT实现

有了上述第三章的公式，我们可以很容易的实现下面的代码。

// FFT.compute
#pragma kernel idft_scale
#pragma kernel dit_x_radix2
#pragma kernel dit_y_radix2
#pragma kernel dit_x_radix2 IDFT
#pragma kernel dit_y_radix2 IDFT

#define NTHREADS 2

#define PI 3.14159265359f

uint n;     // fft 整个大小
uint p;     // 本轮蝶形变换的子FFT的大小，它从 1,2,4 ... To n/2

float scale;  // idft scale

// Texture for FFT
Texture2D<float2> fftin;

// Create a RenderTexture with enableRandomWrite flag and set it
// with cs.SetTexture
RWTexture2D<float2> fftout;

inline float2 mul(float2 a, float2 b) {
    return float2(
        a.x * b.x - a.y * b.y,
        a.x * b.y + a.y * b.x
        );
}

inline float2 exp_theta(float theta) {
    theta *= PI;
    return float2(cos(theta), sin(theta));
}

[numthreads(NTHREADS, NTHREADS, 1)]
void dit_x_radix2(uint3 id: SV_DispatchThreadID) {
    // id.x是本轮蝶形变换的第几个
    // 无论是哪一轮的蝶形变换，蝶形变换的个数都是 n/2
    int2 l_in = id.xy;
    int2 r_in = int2(id.x + int(n / 2), id.y);
    // k是本轮蝶形变换在当前某个Sub FFT中的偏移
    int k = id.x & (p - 1);  // id.x mod p


    int2 l_out = int2((id.x << 1) - k, id.y);
    int2 r_out = int2(l_out.x + p, id.y);
   
#ifdef IDFT
    float theta = float(k) / p;
#else
    float theta = -float(k) / p;
#endif
    float2 f0 = fftin[l_in];
    float2 f1 = mul(fftin[r_in], exp_theta(theta));

    fftout[l_out] = f0 + f1;
    fftout[r_out] = f0 - f1;
}

[numthreads(NTHREADS, NTHREADS, 1)]
void dit_y_radix2(uint3 id: SV_DispatchThreadID) {
    int2 l_in = id.xy;
    int2 r_in = int2(id.x, id.y + int(n / 2));

    int k = id.y & (p - 1);  // id.y mod p

    int2 l_out = int2(id.x, (id.y << 1) - k);
    int2 r_out = int2(id.x, l_out.y + p);

#ifdef IDFT
    float theta = float(k) / p;
#else
    float theta = -float(k) / p;
#endif
    float2 f0 = fftin[l_in];
    float2 f1 = mul(fftin[r_in], exp_theta(theta));

    fftout[l_out] = f0 + f1;
    fftout[r_out] = f0 - f1;
}

[numthreads(NTHREADS, NTHREADS, 1)]
void idft_scale(uint3 id: SV_DispatchThreadID) {
    fftout[id.xy] = fftin[id.xy].xy * scale;
}

上面的Compute Shader在C#中还需要进行一些设置和Dispatch。一个1024x1024的DFT，需要Dispatch 9次x方向，9次y方向，然后还有一次scale操作。

// FFT.cs
public RenderTexture DFT(NormalizeType normalizeType=NormalizeType.DFT)
{
    FFTU4(FFTType.DFT);
    FFTV4(FFTType.DFT);

    if (normalizeType == NormalizeType.DFT)
    {
        Normalize(1.0f / n);
    }
    else if (normalizeType == NormalizeType.SYMMETRIC)
    {
        Normalize(1.0f / Mathf.Sqrt(n));
    }
    else
    {
        SwapTex();
    }

    return fftTextureOut_;
}

public void FFTU(FFTType fftType)
{
    int ditXKernel = fftType == FFTType.DFT ? KERNEL_DIT_X : KERNEL_DIT_X_IDFT;

    uint p = 1;
    if (p < n)
    {
        fftShader_.SetTexture(ditXKernel, SHADER_NAME_FFT_IN, fftTextureIn_);
        fftShader_.SetTexture(ditXKernel, SHADER_NAME_FFT_OUT, fftTextureOut_);
        fftShader_.SetInt(SHADER_NAME_P, (int)p);

        fftTextureOut_.DiscardContents();
        fftShader_.Dispatch(ditXKernel, halfGroupN, groupN, 1);
        SwapTex();
        p <<= 1;

        while (p < n)
        {
            fftShader_.SetTexture(ditXKernel, SHADER_NAME_FFT_IN, fftTextureTmp_);
            fftShader_.SetTexture(ditXKernel, SHADER_NAME_FFT_OUT, fftTextureOut_);
            fftShader_.SetInt(SHADER_NAME_P, (int)p);

            fftTextureOut_.DiscardContents();
            fftShader_.Dispatch(ditXKernel, halfGroupN, groupN, 1);
            SwapTex();
            p <<= 1;
            
        }
    }
}

public void FFTV(FFTType fftType)
{
    int ditYKernel = fftType == FFTType.DFT ? KERNEL_DIT_Y : KERNEL_DIT_Y_IDFT;

    uint p = 1;
    while (p < n)
    {
        fftShader_.SetTexture(ditYKernel, SHADER_NAME_FFT_IN, fftTextureTmp_);
        fftShader_.SetTexture(ditYKernel, SHADER_NAME_FFT_OUT, fftTextureOut_);
        fftShader_.SetInt(SHADER_NAME_P, (int)p);

        fftTextureOut_.DiscardContents();
        fftShader_.Dispatch(ditYKernel, groupN, halfGroupN, 1);
        SwapTex();
        p <<= 1;
    }
}

public void Normalize(float scale)
{
    fftShader_.SetTexture(KERNEL_SCALE, SHADER_NAME_FFT_IN, fftTextureTmp_);
    fftShader_.SetTexture(KERNEL_SCALE, SHADER_NAME_FFT_OUT, fftTextureOut_);
    fftShader_.SetFloat(SHADER_NAME_SCALE, scale);

    fftTextureOut_.DiscardContents();
    fftShader_.Dispatch(KERNEL_SCALE, groupN, groupN, 1);
}

结果