生成采样信号

import numpy as np

f0,f1 = 0.5,2 # 最高频率为f1
T = 1/f0 # 采样时间为最低频率对应的周期
fs = m*f1 # 采样频率为最高频率的m倍
dt = 1/fs # 采样间隔

# 生产采样信号
t = np.arange(0,T,dt)
y = 3 + 2*np.sin(2*np.pi*f0*t+np.pi/2) + np.sin(2*np.pi*f1*t+np.pi/3) 

 编写DFT程序

# 进行dft
N = len(y)
X_dft = np.zeros([N],dtype=complex) # 生成N个元素的向量
for k in range(N):
    for n in range(N):
        # 累加
        X_dft[k] = X_dft[k] + y[n]*np.exp(-2j*np.pi*k*n/N)
# 舍去
X_dft = np.round(X_dft,4)
print(X_dft)

运行已有的FFT函数

进行FFT变换

X_fft = np.fft.fft(y)
X_fft = np.round(X_fft,4)
print(X_fft)

幅值

# 取模后，第一个元素除以N，其余元素除以N/2，才是幅值
X_fft_abs = abs(X_fft)/len(y)*2
X_fft_abs[0] = X_fft_abs[0]/2
X_fft_abs = np.round(X_fft_abs,4)

相位

# 结果加pi/2才是相位
X_fft_phase = np.angle(X_fft) + np.pi/2
X_fft_phase[0] = 0

频率序列

freq = np.fft.fftfreq(len(y),dt)

我们所要的信息是:
f=0,A=3
f=0.5,A=2,Φ=1.5708
f=2,A=1,Φ=1.0472
其他值不是我们关心的。

比如下面这行代码没啥意思，虽然在这没啥意思，但如果要编写程序自动找出幅值与之对应的相位以及频率，这么做会方便查找。

X_fft_phase[0] = 0

对比 已编好的FFT与自编的DFT

第一组 m = 3，N=12

第二组 m = 4，N=16=2^4

因为时间太短不明显，将采样时间调长，T=1/f0*100

第一组，fft用时0.0019s，dft用时5s，第二组，fft用时0.0s，dft用时9s
这只是一次统计，但可以说明fft速度就是快，dft速度特别慢。因为速度太快，所以没有显示出变换序列是2的整数次方，fft会更快。

理论上对于fft算法，变换序列是2的整数次方，会更快。

DFT的快速算法有好多种，有名的是库利－图基快速傅里叶变换算法（Cooley-Tukey）即FFT算法。

https://www.wanweibaike.com/wiki-快速傅里叶变换#其他算法

术语

 混叠：采样率低于最高频率的2倍造成的（欠采样造成的）

能量泄露：采样时间过短造成的（用于FFT的信号过短造成的）

栅栏效应：虽然采样率满足最高频率的2倍，但还是太低造成的

总之，采样率尽量高、采样时间尽量长、采样时间尽量为各分量周期的倍数（整周期采样），就不会发生这些，分解得到的幅值、相位越接近实际。

疑问：在此实验中，采样时间为1/f0，我的采样率为2*f1时，就会出现混叠，导致f=2对应的幅值能量泄露，因此改为3*f1，不知道什么问题？

十大著名算法——出自Guest Editors Introduction to the top 10 algorithms

蒙特卡洛的Metropolis算法
线性规划中的单纯形法
Krylov子空间迭代法
矩阵计算的分解方法
Fotran最优化编译器
计算特征值的QR算法
用于排序的quicksort算法
FFT
整数关系探测
快速多级子法