傅里叶变换 - 3

最新推荐文章于 2021-04-30 19:53:16 发布

yyx_daytoy

最新推荐文章于 2021-04-30 19:53:16 发布

阅读量431

点赞数

分类专栏：傅里叶变换文章标签： python

本文链接：https://blog.csdn.net/yyx_daytoy/article/details/113843161

版权

傅里叶变换专栏收录该内容

4 篇文章 0 订阅

订阅专栏

一多项式的表达方式

系数表示

$p\left ( x \right ) = a_{0} * x ^{0} + a_{1} * x ^{1} + a_{2} *x ^{2} + a_{3} *x ^{3} + ... + a_{N-1} *x ^{N-1}$

数值表式

对于维度是N-1的多项式，对x取N个不同的值，计算对应的N个p值，得到N组对应的值

$\left ( x_{0}, p\left ( x_{0} \right ) \right ), \left ( x_{1}, p\left ( x_{1} \right ) \right ), \left ( x_{2}, p\left ( x_{2} \right ) \right ),...,\left ( x_{N-1}, p\left ( x_{N-1} \right ) \right )$

说明：

文章中 $N = 2^{n}$

二离散傅里叶变换

对于多项式的数值表式，x取值设置为N次复方根，可以得到离散傅里叶变换，也就是

$\left \{ \left ( w_{N}^{0}, p\left ( w_{N}^{0} \right ) \right ), \left ( w_{N}^{1}, p\left ( w_{N}^{1} \right ) \right ), \left ( w_{N}^{2}, p\left ( w_{N}^{2} \right ) \right ), ...,\left ( w_{N}^{N-1}, p\left ( w_{N}^{N-1} \right ) \right )\right \}$

这个简单举例说明一下：以 $p\left ( w_{n}^{2} \right )$ 为例

离散傅里叶公式

$\dpi{120} X(k) = \sum_{n=0}^{N-1} x(n) \left \{ cos(\frac{2\pi kn}{N}) - i*sin(\frac{2\pi kn}{N})\right \}$

那么

$X\left ( 2 \right ) = \sum_{n=0}^{N-1} x\left ( n\right )\left \{ cos\left ( \frac{2\pi n*2}{N} \right ) - i*sin\left ( \frac{2\pi n*2}{N}\right )\right \}$

$p\left ( w_{n}^{2} \right ) = a_{0} + a_{1} * w_{N}^{2} + a_{2} * \left ( w_{N}^{2} \right )^{2} +a_{3} * \left ( w_{N}^{2} \right )^{3}+ ... + a_{N} * \left ( w_{N}^{2} \right )^{N-1}$

$= a_{0} + a_{1} * w_{N}^{2} + a_{2} * \left ( w_{N}^{4} \right )+a_{3} * \left ( w_{N}^{6} \right ) + ... + a_{N} * \left ( w_{N}^{2\left ( N-1 \right )} \right )$

$=\sum_{n=0}^{N-1} a_{n} * \left ( w_{N}^{2n} \right )$

$=\sum_{n=0}^{N-1} a_{n} * \left \{ cos\left ( \frac{2\pi *2n }{N} \right ) - i * sin\left ( \frac{2\pi *2 n }{N}\right )\right \}$

这里用到复方根性质5，说明了 $X\left ( 2 \right )$ 与 $p\left ( w_{n}^{2} \right )$ 相同，可以用上面的数值对表示离散傅里叶变换。具体的证明过程大家可以尝试下，等有时间我把这部分的证明改正一下。

三快速傅里叶变换推导

将多项式构造成两个多项式，一个只包含x的偶数指数部分对应的系数，另一个只包含x的奇数指数部分对应的系数，也就是

$p^{even}\left ( x \right ) = a_{0} + a_{2}*x^{1} + a_{4}*x^{2} +...+ a_{N-2}*x^{\frac{N-2}{2}}$

$p^{odd}\left ( x \right ) = a_{1} + a_{3}*x^{1} + a_{5}*x^{2} +...+ a_{N-1}*x^{\frac{N-2}{2}}$

于是有

$p\left ( x \right ) = p^{even}\left ( x^{2}\right ) + x * p^{odd}\left ( x^{2} \right )$

这样就将一个大问题，分解成两个规模较小的问题，解决两个规模较小的问题后，合起来就解决大规模问题。使用复方根的目的是让 $p^{even}\left ( x \right )$ 和 $p^{odd}\left ( x \right )$ 继续分解下去

在此讨论一下这些复方根的特点：

特点一：把复方根平分为两半，即

$\dpi{120} \left ( w_{N}^{0} \right ), \left ( w_{N}^{1} \right ), ..., \left ( w_{N}^{\frac{N}{2}-1} \right ),$

$\left ( w_{N}^{\frac{N}{2}} \right ), \left ( w_{N}^{\frac{N}{2}+1} \right ), ..., \left ( w_{N}^{\frac{N}{2} + \left ( \frac{N}{2}-1 \right )} \right )$

根据性质5： $w_{N}^{k+j} = w_{N}^{k} * w_{N}^{j}$ 上面的值可以转换为

$\dpi{120} \left ( w_{N}^{0} \right ), \left ( w_{N}^{1} \right ), ..., \left ( w_{N}^{\frac{N}{2}-1} \right ),$

$\left ( w_{N}^{\frac{N}{2}} * w_{N}^{0}\right ), \left ( w_{N}^{\frac{N}{2}} * w_{N}^{1}\right ), ..., \left ( w_{N}^{\frac{N}{2}} * w_{N}^{\frac{N}{2}-1} \right )$

根据性质2： $w_{N}^{\frac{N}{2}} = -1$ 上面的值可以转换为

$\dpi{120} \left ( w_{N}^{0} \right ), \left ( w_{N}^{1} \right ), ..., \left ( w_{N}^{\frac{N}{2}-1} \right ),$

$\left ( -w_{N}^{0} \right ), \left ( -w_{N}^{1} \right ), ..., \left ( -w_{N}^{\frac{N}{2}-1} \right )$

也就是说，在复方根中，后一半的数值是前一半数值的相反数。

特点二：将x取值为复方根，计算每个根的平方

$\left ( w_{N}^{0} \right )^{2}, \left ( w_{N}^{1} \right )^{2}, ..., \left ( w_{N}^{\frac{N}{2}-1} \right )^{2},$

$\left ( w_{N}^{\frac{N}{2}} \right )^{2}, \left ( w_{N}^{\frac{N}{2}+1} \right )^{2}, ..., \left ( w_{N}^{\frac{N}{2} + \left ( \frac{N}{2}-1 \right )} \right )^{2}$

将复方根分为上下两部分比较，根据复方根性质3: ${\left ( w_{N}^{k+\frac{N}{2}}\right )^{2}} = w_{N}^{2k+N} = w_{N}^{2k} = \left ( w_{N}^{k} \right )^{2}$ ，上下两部分平方值相等

四快速傅里叶变换python代码实现

from math import cos, pi, sin
import numpy as np 


def create_complex_roots(N):
    '''
    生成N次复方根
    '''
    roots = []
    for k in range(N):
        root_k = complex(cos(2*k*pi/N), -sin(2*k*pi/N))
        roots.append(root_k)
    return roots


def fast_fourier_transform(coefficients):
    if len(coefficients) <= 1:
        return coefficients
    degree = len(coefficients)
    n_complex_roots = create_complex_roots(degree)    #1 获取N次复方根
    even_coefficients = coefficients[0:degree:2]      #2 获取偶数项系数
    odd_coefficients = coefficients[1:degree:2]       #3 获取奇数项系数
    even_fft = fast_fourier_transform(even_coefficients)  #4 计算偶数项fft
    odd_fft = fast_fourier_transform(odd_coefficients)    #5 计算奇数项fft

    fft_first_part = []
    fft_second_part = []
    for k in range(int(degree/2)): #6 将奇数项和偶数项合并起来
        w = n_complex_roots[k]
        fft_first_part.append(even_fft[k] + w * odd_fft[k])  #7
        fft_second_part.append(even_fft[k] - w * odd_fft[k])  #8 -w是因为后一半根是前一半根的相反数
    fft = []
    fft.extend(fft_first_part)
    fft.extend(fft_second_part)
    return fft


def py_fft(inpt, N):
    fft_result = np.fft.rfft(inpt, N)
    return fft_result


if __name__ == "__main__":
    inpt = [1,2,3,4,5,6,7,8]
    fft = fast_fourier_transform(inpt)
    py_res = py_fft(inpt, len(inpt))
    for item in fft:
        print(item)
    print("\n================================\n")
    for item in py_res:
        print(item)

运行结果是：

(36+0j)
(-4+9.65685424949238j)
(-4+4j)
(-4+1.6568542494923797j)
(-4+0j)
(-4-1.6568542494923806j)
(-3.9999999999999996-4j)
(-3.9999999999999987-9.65685424949238j)

================================

(36+0j)
(-4+9.65685424949j)
(-4+4j)
(-4+1.65685424949j)
(-4+0j)

五 python代码说明

根据代码来理解快速傅里叶变换的算法逻辑。

函数fast_fourier_transform()接收多项式系数组成的数组，如果传入的参数数组只有一个元素，意味着多项式只有一个系数 $a_{0}$ ，代码直接将该元素返回。

如果参数数组元素个数大于1，根据#1创建对应个数的复方根。接着#2和#3抽取偶数项和奇数项对应的系数，然后#4和#5将问题分解成两个规模更小的问题进行递归求解

需要注意的是，能够递归依赖于N次方根的特性。当执行#4语句时，进行递归再次进入函数fast_fourier_transform()，在函数起始处又构造N/2个复方根。在递归前，代码希望对x取值 $w_{N}^{0}, w_{N}^{2}, ..., w_{N}^{N-2}$ 进行计算，递归后代码是对 $w_{\frac{N}{2}}^{0}, w_{\frac{N}{2}}^{1}, ..., w_{\frac{N}{2}}^{\frac{N-2}{2}}$ 进行计算。根据复方根的性质1： $w_{N*d}^{k*d} = w_{N}^{k}$ , 就有 $w_{\frac{N}{2}}^{0}, w_{\frac{N}{2}}^{1}, ..., w_{\frac{N}{2}}^{\frac{N-2}{2}}$ = $w_{\frac{N}{2}*2}^{0*2}, w_{\frac{N}{2}*2}^{1*2}, ..., w_{\frac{N}{2}*2}^{\frac{N-2}{2}*2}$ = $w_{N}^{0}, w_{N}^{2}, ..., w_{N}^{N-2}$ 。所以递归前x的取值和递归时x的取值是一样的，这是由N次复方根的性质决定的。如果x取值不是复方根，无法保证递归前后x的取值一致，这就是要使用复方根进行快速傅里叶变换的原因。

注意#6只循环系数个数的一半，这个因为复方根在取平方时，前一半值与后一半值完全相同。

注意#7和#8，因为只循环系数个数的一半，复方根只有前一半值参与计算。又因为复方根后一半值是前一半值的相反数，w对应复方根前一半值，-w对应复方根后一半值。这就是#7中是+，而#8中是 - 的原因。

此时的FFT的复杂度是O(n*lg(n))。用T(n)表示计算维度是n的多项式所需的时间，那么计算 $p^{evev}$ , $p^{odd}$ 的时间就是 2* T(n/2)，#6对应的循环是O(n)，整个复杂度就是 $T(n) = 2 * T(\frac{n}{2}) + O(n)$ ,使用猜测检验法可以确定 T(n) = O(n*lg(n))。