dsp 有fft程序吗_Python 中 FFT 快速傅里叶分析

行文思路：

采样频率和采样定理
生成信号并做FFT 变换
频率分辨率和显示分辨率
FFT 归一化操作
对噪声信号进行FFT
导入自定义模块
总结

一，相关定理介绍

1，采样频率

采样频率，也称为采样速度或者采样率，定义了每秒从连续信号中提取并组成离散信号的采样个数，它用赫兹（Hz）来表示。采样频率的倒数是采样周期或者叫作采样时间，它是采样之间的时间间隔。通俗的讲采样频率是指计算机每秒钟采集多少个信号样本。

2，采样定理

所谓采样定理，又称香农采样定理，奈奎斯特采样定理，是信息论，特别是通讯与信号处理学科中的一个重要基本结论。采样定理指出，如果信号是带限的，并且采样频率高于信号带宽的两倍，那么，原来的连续信号可以从采样样本中完全重建出来。

定理的具体表述为：在进行模拟/数字信号的转换过程中，当采样频率fs大于信号中最高频率fmax的2倍时,即 fs>2*fmax

Note：采样之后的数字信号完整地保留了原始信号中的信息，一般实际应用中保证采样频率为信号最高频率的2.56～4倍；采样定理又称奈奎斯特定理。

二，FFT 变换

numpy中有一个fft的库，scipy中也有一个fftpack的库，各自都有fft函数，两者的用法基本是一致的。下面内容是利用 scipy 库中的FFT 对信号进行快速傅里叶变化。

1，生成信号做FFT

import numpy as np
from scipy.fftpack import fft
import matplotlib.pyplot as plt
Fs =10e3                                   # 采样频率
f1 =390                                    # 信号频率1
f2 = 2e3                                   # 信号频率2
t=np.linspace(0,1,Fs)                      # 生成 1s 的时间序列
y=2*np.sin(2*np.pi*f1*t)+5*np.sin(2*np.pi*f2*t)+noise2

2，插值补零

进行FFT 运算，实际就一个关键的函数fft（），在进行变换之前先做小的处理。求出信号长度的下一个最近二次幂。类似于MATLAB 中的Nextpower2(), 为了更加精确快速的进行傅里叶变换。

L = len (data) # 信号长度
N =np.power(2,np.ceil(np.log2(L))) # 下一个最近二次幂

求 2 的下一个整数次幂，有更简单快速的方法。这里采用最为易懂的方式实现的。

3，fft() 参数说明

下面为fft（）函数的参数说明。比较重要的是关注下输入参数的类型为数组，其中第二个参数n 指的是做FFT 的点数，如果n 小于信号的长度，输入无效。当n大于信号长度时，超出的部分补零。

Note: 补零并没有提高信号实际的频率分辨率(df=Fs/N)，只是提高了显示分辨了（df=Fs/L）。要想提高实际分辨率可行的办法是:

相同时间下提高采样率
相同采样率下增长采样时间。

下面为完整的FFT 函数。

def FFT (Fs,data):L = len (data)                          # 信号长度N =np.power(2,np.ceil(np.log2(L)))      # 下一个最近二次幂FFT_y = (fft(y,N))/L*2                  # N点FFT，但除以实际信号长度 LFre = np.arange(int(N/2))*Fs/N          # 频率坐标FFT_y = FFT_y[range(int(N/2))]          # 取一半return Fre, FFT_y

4，归一化操作

通常有N个点做FFT就有 N 个复数点与之对应，此时幅值是对复数取绝对值运算。除第一个点之外，实际的幅值是信号的实际长度L，再乘以2。所以一般对信号先做FFT 再取绝对值然后除以L乘以2。

Note: 第一个频率点我信号中的直流分量，对于交流信号来讲直流分量为0，FFT变换后第一个频点为0。

下面以一个例子说明归一化操作。

例如 A1 = 3, A2 =2, A3=5

y=3+2*np.sin(2*np.pi*f1*t)+5*np.sin(2*np.pi*f2*t).

例子中信号的长度为10e3个点。

归一化后图形如下，可以看到第一个峰值为原来的 L 倍

第二个峰为原来的 N/2

第三个峰为原来的N/2

第一个峰为直流分量，其余的峰为信号中包含的频率。

对信号归一化并乘以2 之后为：

信号中加入随机噪声

noise1 = np.random.random(10e3)
# 0-1 之间的随机噪声
noise2 = np.random.normal(1,10,10e3)
#产生的是一个10e3的高斯噪声点数组集合（均值为：1，标准差：10）

当加入 0-1 之间的随机噪声时，可以看到频率0点的数值变大，说明噪声部分相当于直流分量。在加入高斯噪声后，其发现其均值可以理解为直流分量，改变均值的时候频率点0 发生相应变化。标准差表示了噪声之间的离散程度。

加入 0-1 之间随机噪声频谱

加入高斯噪声均值为1，标准差为10 频谱

三，导入自定义模块

1，执行文件和包含模块的文件夹同一个目录下

FFT 模块文件夹中包含：

运行主函数 main.py，将定义的函数写在 __init__.py 中。

通过导入自定义函数：

import FFT
Fre, FFT_y1 = FFT.FFT(Fs,y) # 第一个FFT 为文件名 第二个FFT为函数名
from FFT import FFT         # 另一种导入
Fre, FFT_y1 = FFT(Fs,y)

2，导入的模块不在统一目录下

即执行文件在 main文件夹，模块文件在FFT文件夹，FFT 和 main 的父文件夹为 Python。

这时候需要通过导入sys 模块设定导入模块的路径

import sys
sys.path.append('E:Python')

这里导入路径写到模块的父文件夹即可。

总结：

FFT 变化是信号从时域变化到频域的桥梁，是信号处理的基本方法。本文讲述了利用Python SciPy 库中的fft() 函数进行傅里叶变化，其关键是注意信号输入的类型为np.array 数组类型，以及FFT 变化后归一化和取半操作，得到信号真实的幅度值。注意如果FFT点数和实际信号长度不一样，则归一化时除以信号的实际长度而不是 FFT的点数。

完整代码：

import numpy as np
from scipy.fftpack import fft
import matplotlib.pyplot as plt
Fs =10e3      # 采样频率
f1 =390       # 信号频率1
f2 = 2e3      # 信号频率2
t=np.linspace(0,1,Fs)   # 生成 1s 的实践序列
noise1 = np.random.random(10e3)      # 0-1 之间的随机噪声
noise2 = np.random.normal(1,10,10e3)
#产生的是一个10e3的高斯噪声点数组集合（均值为：1，标准差：10）
y=2*np.sin(2*np.pi*f1*t)+5*np.sin(2*np.pi*f2*t)+noise2def FFT (Fs,data):L = len (data)                        # 信号长度N =np.power(2,np.ceil(np.log2(L)))    # 下一个最近二次幂FFT_y1 = np.abs(fft(data,N))/L*2      # N点FFT 变化,但处于信号长度Fre = np.arange(int(N/2))*Fs/N        # 频率坐标FFT_y1 = FFT_y1[range(int(N/2))]      # 取一半return Fre, FFT_y1Fre, FFT_y1 = FFT(Fs,y)
plt.figure
plt.plot(Fre,FFT_y1)
plt.grid()
plt.show()

导入模块参考blog.csdn.net