数值计算利用Python实现快速傅里叶变换(FFT)(不调用函数库)

Python实现快速傅里叶变换(FFT)详解

最新推荐文章于 2024-03-16 11:24:43 发布

原创最新推荐文章于 2024-03-16 11:24:43 发布 · 3.5k 阅读

30 ·

CC 4.0 BY-SA版权

文章标签：

#python #numpy #傅立叶分析

数值计算专栏收录该内容

5 篇文章

订阅专栏

该博客介绍了如何在Python中不依赖函数库实现快速傅里叶变换(FFT)。首先，解释了离散傅里叶变换(DFT)的基本概念和公式，接着详细阐述了FFT算法，降低了计算复杂度。博客还提供了具体的代码实现，并给出了使用实例，例如在[−π,π]上对f(x)=x^2cosx进行16次三角插值多项式的FFT计算。" 112066581,10539879,决策树：优化选择与生活决策,"['机器学习', '决策树算法', '信息论', '数据挖掘', '分类问题']

数值计算-利用Python实现快速傅里叶变换(FFT)(不调用函数库)

快速傅里叶变换介绍

首先介绍离散傅里叶变换（DFT），不管是正变换还是逆变换，计算公式可以归结为
$X(k)=\sum_{n=0}^{N-1}x(n)W_N^{nk}$
其中, $N$ 为抽样点的个数(一般为 $2^n$ ), $x (n)$ 为已知数据（拟合函数时为节点函数值)， $WN=e−i2πnW_N=e^{-i\frac{2\pi}{n}}$ 为旋转因子。

而快速傅里叶变换（FFT），只是对DFT在算法上进行一定的幼化，降低复杂度，同样可以总结为一个公式:
$,NX(k)=\left\{ \begin{aligned} &\sum_{n=0}^{N/2-1}x_1(r)W_{N/2}^{rk}+W_N^{k}\sum_{n=0}^{N/2-1}x_2(r)W_{N/2}^{rk} \quad k=0,1,\cdots,\frac{N}{2}-1\\ &\sum_{n=0}^{N/2-1}x_1(r)W_{N/2}^{rk'}+W_N^{k'}\sum_{n=0}^{N/2-1}x_2(r)W_{N/2}^{rk'} \quad k=k'+\frac{N}{2}=\frac{N}{2},\frac{N}{2}+1,\cdots,N \end{aligned} \right.$

代码实现

import numpy as np
from numpy import pi
from numpy import exp
from numpy import cos, sin
import matplotlib.pyplot as plt
from pylab import mpl


def func(p: float):
    return (p - pi) ** 2 * np.cos((p - pi))  # x = y - pi


def get_function(p: float, c: list, N: int):
    a = [_ for _ in range(0, N)]
    b = [_ for _ in range(0, N)]
    res = 0
    for k in range(0, N):
        a[k] = (2 * c[k]).real
        b[k] = -(2 * c[k]).imag
    # print(a)
    for k in range(1, N // 2):
        res += a[k] * cos(k * p)
        res += b[k] * sin(k * p)

    res += a[0] / 2
    res += b[N // 2] * sin(N // 2 * p)
    return res


def solve_c(y: list, f: list, N: int):
    c = [0 for _ in range(0, N)]

    for k in range(0, N):
        if k < N / 2:
            c[k] = sum(f[0::2][i] / N * exp(-k * i * pi * 2j / (N // 2)) for i in range(0, N // 2)) \
                   + sum(f[1::2][i] / N * exp(-k * i * pi * 2j / (N // 2)) for i in range(0, N // 2)) \
                   * exp(-k * pi * 2j / N)
            # c[k] *= (-1) ** k
        else:
            kk = k - N / 2
            c[k] = sum(f[0::2][i] / N * exp(-kk * i * pi * 2j / (N // 2)) for i in range(0, N // 2)) \
                   - sum(f[1::2][i] / N * exp(-kk * i * pi * 2j / (N // 2)) for i in range(0, N // 2)) \
                   * exp(-kk * pi * 2j / N)

    # print(c)
    return c


def draw(data: np.ndarray, x: list, f: list, N: int):
    c = solve_c(y, f, N)
    x = [_ - pi for _ in y]
    plt.plot(data - pi, get_function(data, c, N), label="%s次三角插值拟合曲线" % str(N // 2), color="blue")
    plt.scatter(x, f, label="节点数据", color="red")
    plt.title("%s次三角插值拟合曲线" % str(N // 2))
    mpl.rcParams['font.sans-serif'] = ['SimHei']
    mpl.rcParams['axes.unicode_minus'] = False
    plt.legend(loc="upper right")

使用实例

考虑 $f(x)=x^2\cos{x}$ 在 $[−π,π][-\pi,\pi]$ 上的16次三角插值多项式

if __name__ == "__main__":
    N = 2 * 16
    y = [2 * pi / N * _ for _ in range(0, N)]
    f = [func(_) for _ in y]
    draw(np.arange(0, 2 * pi + 0.001, 0.001), y, f, N)
    plt.plot(np.arange(-pi, pi + 0.001, 0.001), func(np.arange(0, 2 * pi + 0.001, 0.001)), color="black")
    plt.show()