基于Python的计算物理学习笔记第二章进阶计算方法-CSDN博客

本文链接：https://blog.csdn.net/ihdbdnw/article/details/134988703

一、傅立叶变换

import numpy as np
import matplotlib.pyplot as plt
from scipy import fftpack

1. 离散傅立叶变换（DFT）

将 $\vec a$ 傅立叶变换为 $\vec c$ ，可通过一个傅立叶矩阵进行变换：

$\vec c=\bold F\vec a$

$\bold F=\begin{pmatrix} 1 &1 &... &1 \\ 1& \omega &... &\omega ^{N-1} \\ ... & ...&... &... \\ 1 &\omega ^{N-1} &... &\omega ^{(N-1)^2} \end{pmatrix}$

矩阵元为 $F_{ij}=\omega^{jk},(j,k=0,1,...,N-1)$ ，其中 $\omega\equiv e^{2\pi i\over N}$

也就是说， $F_{ij}=e^{-{2\pi i\over N}jk}$

$c_k=\sum\limits_{j=1}^{N-1}a_je^{-{2\pi i\over N}jk}$

这一方法的时间复杂度为 $O(N^2)$

2. 快速傅立叶变换（FFT）

在DFT的基础上，注意到当j*k分别为N/2、N、3N/2、2N时，有 $\omega^{jk}$ 值的将会依次为 $-i,-1,-i,1$ 。因此我们可以尝试利用 $e^{-{2\pi i\over N}jk}$ 的周期性：

当将 k 为偶数时， $k=2n$ ，根据周期性有 $e^{-{2\pi i\over N}2nj}=e^{-i{2\pi n\over N/2}j}=e^{-i{2\pi n\over N/2}(j+N/2)}$

因此当 k 为偶数时，有 $\omega^{jk}=\omega^{(j+N/2)k}$ ，也就是说第j行与j+N/2行的矩阵元是相同的。

同理，当 k 为奇数时，我们只需要分离出一个 $e^{-{2\pi i\over N}j}$ 的公因数，剩下的部分也满足偶数的条件。

因此，我们只需要将奇偶列分开并各自利用 $\omega^{jk}=\omega^{(j+N/2)k}$ 的性质求解，就能极大的加速算法，将时间复杂度由 $O(N^2)$ 降为 $O(N\log N)$

3. 利用scipy进行傅立叶变换

（1）正变换

x = np.linspace(0,2*np.pi,20000,)

a = 30*np.sin(x) + 20*np.cos(10*x) + 10*np.sin(20*x) + 5*np.cos(30*x)

c = fftpack.fft(a)/len(a) #傅立叶变换，除以len(a)是为了归一化

k = fftpack.fftfreq(len(a)) #傅立叶曲线的频率轴

plt.plot(x,a)

plt.plot(k,abs(c))

fft函数计算的模数太多，以至于分辨不出图像。示例给出的周期函数用不了那么多模。

plt.plot(k[:50],abs(c)[:50]) #手动选取前50个mode