python绘制频谱图_使⽤python进⾏傅⾥叶FFT-频谱分析详细
教程
⽬录
⼀、⼀些关键概念的引⼊
1、离散傅⾥叶变换(DFT)
2、快速傅⾥叶变换(FFT)
3、采样频率以及采样定理
4、如何理解采样定理？
⼆、使⽤scipy包实现快速傅⾥叶变换
1、产⽣原始信号——原始信号是三个正弦波的叠加
2、快速傅⾥叶变换
3、FFT的原始频谱
4、将振幅谱进⾏归⼀化和取半处理
三、完整代码
⼀、⼀些关键概念的引⼊
1、离散傅⾥叶变换(DFT)
离散傅⾥叶变换(discrete Fourier transform) 傅⾥叶分析⽅法是信号分析的最基本⽅法，傅⾥叶变换是傅⾥叶分析的核⼼，通过它把信号从时间域变换到频率域，进⽽研究信号的频谱结构和变化规律。但是它的致命缺点是：计算量太⼤，时间复杂度太⾼，当采样点数太⾼的时候，计算缓慢，由此出现了DFT的快速实现，即下⾯的快速傅⾥叶变换FFT。
2、快速傅⾥叶变换(FFT)
计算量更⼩的离散傅⾥叶的⼀种实现⽅法。详细细节这⾥不做描述。
3、采样频率以及采样定理
采样频率：采样频率，也称为采样速度或者采样率，定义了每秒从连续信号中提取并组成离散信号的采样个数，它⽤赫兹(Hz)来表⽰。采样频率的倒数是采样周期或者叫作采样时间，它是采样之间的时间间隔。通俗的讲采样频率是指计算机每秒钟采集多少个信号样本。
采样定理：所谓采样定理 ，⼜称⾹农采样定理，奈奎斯特采样定理，是信息论，特别是通讯与信号处理学科中的⼀个重要基本结论。采样定理指出，如果信号是带限的，并且采样频率⾼于信号带宽的两倍，那么，原来的连续信号可以从采样样本中完全重建出来。
定理的具体表述为：在进⾏模拟/数字信号的转换过程中，当采样频率fs⼤于信号中最⾼频率fmax的2倍时，即
fs>2*fmax
采样之后的数字信号完整地保留了原始信号中的信息，⼀般实际应⽤中保证采样频率为信号最⾼频率的2.56～4倍；采样定理⼜称奈奎斯特定理。
4、如何理解采样定理？
在对连续信号进⾏离散化的过程中，难免会损失很多信息，就拿⼀个简单地正弦波⽽⾔，如果我1秒内就选择⼀个点，很显然，损失的信号太多了，光着⼀个点我根本不知道这个正弦信号到底是什么样⼦
的，⾃然也没有办法根据这⼀个采样点进⾏正弦波的还原，很明显，我采样的点越密集，那越接近原来的正弦波原始的样⼦，⾃然损失的信息越少，越⽅便还原正弦波。故⽽
采样定理说明采样频率与信号频率之间的关系，是连续信号离散化的基本依据。 它为采样率建⽴了⼀个⾜够的条件，该采样率允许离散采样序列从有限带宽的连续时间信号中捕获所有信息。
⼆、使⽤scipy包实现快速傅⾥叶变换
本节不会说明FFT的底层实现，只介绍scipy中fft的函数接⼝以及使⽤的⼀些细节。
1、产⽣原始信号——原始信号是三个正弦波的叠加
import numpy as npfrom scipy.fftpack import fft,ifftimport matplotlib.pyplot as pltfrom matplotlib.pylab import mpl
600y=7*np.sin(2*np.pi*200*x) + 5*np.sin(2*np.pi*400*x)+3*np.sin(2*np.pi*600*x)
这⾥原始信号的三个正弦波的频率分别为，200Hz、400Hz、600Hz，最⼤频率为600赫兹。根据采样定理，fs⾄少是600赫兹的2倍，这⾥选择1400赫兹，即在⼀秒内选择1400个点。
原始的函数图像如下：plt.figure()plt.plot(x,y) plt.title('原始波形') plt.figure()plt.plot(x[0:50],y[0:50]) plt.title('原始部分波形(前50组样本)')plt.show()
由图可见，由于采样点太过密集，看起来不好看，我们只显⽰前⾯的50组数据，如下：
2、快速傅⾥叶变换
其实scipy和numpy⼀样，实现FFT⾮常简单，仅仅是⼀句话⽽已，函数接⼝如下：
from scipy.fftpack import fft,ifft
from numpy import fft,ifft
其中fft表⽰快速傅⾥叶变换，ifft表⽰其逆变换。具体实现如下：
fft_y=fft(y) #快速傅⾥叶变换print(len(fft_y))print(fft_y[0:5])'''运⾏结果如下：1400[-4.18864943e-12+0.j 9.66210986e-05-0.04305756j 3.86508070e-04-0.08611996j 8.69732036e-04-0.12919206j 1.54641157e-03-0.17227871j]'''
我们发现以下⼏个特点：
(1)变换之后的结果数据长度和原始采样信号是⼀样的
(2)每⼀个变换之后的值是⼀个复数，为a+bj的形式，那这个复数是什么意思呢？
我们知道，复数a+bj在坐标系中表⽰为(a,b)，故⽽复数具有模和⾓度，我们都知道快速傅⾥叶变换具有
“振幅谱”“相位谱”，它其实就是通过对快速傅⾥叶变换得到的复数结果进⼀步求出来的，
那这个直接变换后的结果是不是就是我需要的，当然是需要的，在FFT中，得到的结果是复数，
(3)FFT得到的复数的模(即绝对值)就是对应的“振幅谱”，复数所对应的⾓度，就是所对应的“相位谱”，现在可以画图了。
3、FFT的原始频谱N=1400x = np.arange(N) # 频率个数 abs_y=np.abs(fft_y) # 取复数的绝对值，即复数的模(双边频
谱)angle_y=np.angle(fft_y) #取复数的⾓度 plt.figure()plt.plot(x,abs_y) plt.title('双边振幅谱(未归⼀化)')
plt.figure()plt.plot(x,angle_y) plt.title('双边相位谱(未归⼀化)')plt.show()
显⽰结果如下：
注意：我们在此处仅仅考虑“振幅谱”，不再考虑相位谱。
我们发现，振幅谱的纵坐标很⼤，⽽且具有对称性，这是怎么⼀回事呢？
关键：关于振幅值很⼤的解释以及解决办法——归⼀化和取⼀半处理
⽐如有⼀个信号如下：
Y=A1+A2*cos(2πω2+φ2)+A3*cos(2πω3+φ3)+A4*cos(2πω4+φ4)
经过FFT之后，得到的“振幅图”中，
第⼀个峰值(频率位置)的模是A1的N倍，N为采样点，本例中为N=1400，此例中没有，因为信号没有常数项A1
第⼆个峰值(频率位置)的模是A2的N/2倍，N为采样点，
第三个峰值(频率位置)的模是A3的N/2倍，N为采样点，
第四个峰值(频率位置)的模是A4的N/2倍，N为采样点，
依次下去......
考虑到数量级较⼤，⼀般进⾏归⼀化处理，既然第⼀个峰值是A1的N倍，那么将每⼀个振幅值都除以N即可
FFT具有对称性，⼀般只需要⽤N的⼀半，前半部分即可。
4、将振幅谱进⾏归⼀化和取半处理
先进⾏归⼀化matplotlib中subplot
normalization_y=abs_y/N #归⼀化处理(双边频谱)plt.figure()plt.plot(x,normalization_y,'g')plt.title('双边频谱(归⼀
化)',fontsize=9,color='green')plt.show()
结果为：
现在我们发现，振幅谱的数量级不⼤了，变得合理了，接下来进⾏取半处理：half_x = x[range(int(N/2)
)] #取⼀半区间
normalization_half_y = normalization_y[range(int(N/2))] #由于对称性，只取⼀半区间(单边频
谱)plt.figure()plt.plot(half_x,normalization_half_y,'b')plt.title('单边频谱(归⼀化)',fontsize=9,color='blue')plt.show()
这就是我们最终的结果，需要的“振幅谱”。
三、完整代码
import numpy as npfrom scipy.fftpack import fft,ifftimport matplotlib.pyplot as pltfrom matplotlib.pylab import mpl
600y=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) #快速傅
⾥叶变换 N=1400x = np.arange(N) # 频率个数half_x = x[range(int(N/2))] #取⼀半区间 abs_y=np.abs(fft_y) # 取复数的绝对值，即复数的模(双边频
谱)angle_y=np.angle(fft_y) #取复数的⾓度normalization_y=abs_y/N #归⼀化处理(双边频谱) normalization_half_y = normalization_y[range(int(N/2))] #由于对称性，只取⼀半区间(单边频谱) plt.subplot(231)plt.plot(x,y) plt.title('原始波形')
plt.subplot(232)plt.plot(x,fft_y,'black')plt.title('双边振幅谱(未求振幅绝对值)',fontsize=9,color='black')
plt.subplot(233)plt.plot(x,abs_y,'r')plt.title('双边振幅谱(未归⼀化)',fontsize=9,color='red')
plt.subplot(234)plt.plot(x,angle_y,'violet')plt.title('双边相位谱(未归⼀化)',fontsize=9,color='violet')
plt.subplot(235)plt.plot(x,normalization_y,'g')plt.title('双边振幅谱(归⼀化)',fontsize=9,color='green')
plt.subplot(236)plt.plot(half_x,normalization_half_y,'blue')plt.title('单边振幅谱(归⼀化)',fontsize=9,color='blue') plt.show()

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。