Python 实现简单的快速傅里叶变换

---

前言

 相比于MATLAB自带的FFT函数以及详尽的官方文档来说，python在傅里叶变换这个方面相比就不是那么简单了，处处需要使用Help查看相关函数的定义。但是本质来说，都是傅里叶变换，只是编程语言不同而已。

---

一、使用到的python库

import numpy as np
from scipy.fftpack import fft, ifft,fftshift
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl

 有关下面几个库，怎么下载，我提供两种方式
①打开Pycharm，依次打开File->Settings->搜索Python Interpreter->通过添加库进行安装。
②Pip模式安装，有时候国内网络直连往往存在安装不上的问题，我们便需要通过cmd命令通过pip-install的方式来进行。下面，我首先给出一些国内的镜像站供大家参考

[1] 阿里云 http://mirrors.aliyun.com/pypi/simple/
[2] 豆瓣http://pypi.douban.com/simple/
[3] 清华大学 https://pypi.tuna.tsinghua.edu.cn/simple/
[4] 中国科学技术大学 http://pypi.mirrors.ustc.edu.cn/simple/
[5] 华中科技大学http://pypi.hustunique.com/

 怎么使用呢？下面我给出一个例子（在cmd命令下进行）:

pip install -i https://pypi.tuna.tsinghua.edu.cn/simple selenium

二、全部示例代码及解释

1.代码

import numpy as np
from scipy.fftpack import fft, ifft,fftshift
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl
# 前面是对于库函数的导入
mpl.rcParams['font.sans-serif'] = ['SimHei']  # 显示中文
mpl.rcParams['axes.unicode_minus'] = False  # 显示负号
# 采样点选择1400个，因为设置的信号频率分量最高为600赫兹，根据采样定理知采样频率要大于信号频率2倍，所以这里设置采样频率为1400赫兹（即一秒内有1400个采样点，一样意思的）
fs=1400
N=1400
x = np.linspace(0, 1, fs)#生成时间刻度,生成原则是0到1，产生fs-1个样点，将区间分为fs份，共fs个点
# 设置需要采样的信号，频率分量有200，400和600
y = 7 * np.sin(2 * np.pi * 200 * x) + 5 * np.sin(2 * np.pi * 400 * x) + 3 * np.sin(2 * np.pi * 600 * x)
fft_y = fft(y)  # 快速傅里叶变换
flabel=fs*np.arange(N)/N
flabel1 = flabel-fs/2 #双边谱频率刻度
flabel2 = flabel[range(int(N/2))]#单边谱频率刻度
bothside = fftshift(np.abs(fft_y) ) # 取复数的绝对值，即复数的模(双边频谱)，fft作用于原频域样点，使之满足双边的条件
singleside = np.abs(fft_y)#单边谱
angle_y = np.angle(fft_y)  # 取复数的角度
normalization_y = bothside / N  # 归一化处理（双边频谱）
normalization_singleside=singleside/N;
normalization_half_y = normalization_singleside[range(int(N / 2))]  # 由于对称性，只取一半区间（单边频谱）
plt.subplot(231)
plt.plot(x, y)
plt.title('原始波形')
plt.subplot(232)
plt.plot(flabel1, fft_y, 'black')
plt.title('双边振幅谱(未求振幅绝对值)', fontsize=9, color='black')
plt.subplot(233)
plt.plot(flabel1, bothside, 'r')
plt.title('双边振幅谱(未归一化)', fontsize=9, color='red')
plt.subplot(234)
plt.plot(flabel1, angle_y, 'violet')
plt.title('双边相位谱(未归一化)', fontsize=9, color='violet')
plt.subplot(235)
plt.plot(flabel1, normalization_y, 'g')
plt.title('双边振幅谱(归一化)', fontsize=9, color='green')
plt.subplot(236)
plt.plot(flabel2, normalization_half_y, 'blue')
plt.title('单边振幅谱(归一化)', fontsize=9, color='blue')
plt.show()

2.部分函数的解释

①fft
 FFT就是DFT的快速计算，定义式如下
D F T { x [ n ] } = ∑ n = 0 N − 1 x [ n ] e − j 2 π N n \mathrm{DFT}\left\{ x\left[ n \right] \right\} =\sum_{n=0}^{N-1}{x\left[ n \right] e^{-j\frac{2\pi}{N}n}} DFT{x[n]}=n=0∑N−1​x[n]e−jN2π​n
②fftshit
 这个函数在Matlab也有出现，主要是做什么的呢，你们是否还记得FFT(DFT)计算的结果是从$\left[0，2\pi \right]$ ，以$\pi$为周期。对应的频率范围为0频到最高频，再从负频率到0，fftshift做的就是把频率点，按照负频率到正频率重新排布，那么我们生成对于的频率标度与之对应就可以形成对一个信号的双边谱分析。