傅里叶变换通俗解释及快速傅里叶变换的python实现

　　通俗理解傅里叶变换，先看这篇文章傅里叶变换的通俗理解！

　　接下来即是使用python进行傅里叶FFT-频谱分析：

1、一些关键概念的引入

一、离散傅里叶变换（DFT）

       离散傅里叶变换(discrete Fourier transform) 傅里叶分析方法是信号分析的最基本方法，傅里叶变换是傅里叶分析的核心，经过它把信号从时间域变换到频率域，进而研究信号的频谱结构和变化规律。可是它的致命缺点是：计算量太大，时间复杂度过高，当采样点数过高的时候，计算缓慢，由此出现了DFT的快速实现，即下面的快速傅里叶变换FFT。

二、快速傅里叶变换（FFT）

       计算量更小的离散傅里叶的一种实现方法。详细细节这里不作描述。

三、采样频率以及采样定理

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

采样定理：所谓采样定理 ，又称香农采样定理，奈奎斯特采样定理，是信息论，特别是通信与信号处理学科中的一个重要基本结论。采样定理指出，若是信号是带限的，而且采样频率高于信号带宽的两倍，那么，原来的连续信号能够从采样样本中彻底重建出来。

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

　　fs>2*fmax

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

四、如何理解采样定理？

      在对连续信号进行离散化的过程当中，不免会损失不少信息，就拿一个简单地正弦波而言，若是我1秒内就选择一个点，很显然，损失的信号太多了，光着一个点我根本不知道这个正弦信号究竟是什么样子的，天然也没有办法根据这一个采样点进行正弦波的还原，很明显，我采样的点越密集，那越接近原来的正弦波原始的样子，天然损失的信息越少，越方便还原正弦波。故而

       采样定理说明采样频率与信号频率之间的关系，是连续信号离散化的基本依据。 它为采样率创建了一个足够的条件，该采样率容许离散采样序列从有限带宽的连续时间信号中捕获全部信息。

 

2、使用scipy包实现快速傅里叶变换

      本节不会说明FFT的底层实现，只介绍scipy中fft的函数接口以及使用的一些细节。

一、产生原始信号——原始信号是三个正弦波的叠加

import numpy as np
from scipy.fftpack import fft,ifft
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个采样点，同样意思的）
x=np.linspace(0,1,1400)      
 
#设置须要采样的信号，频率份量有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)

这里原始信号的三个正弦波的频率分别为，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组数据，以下：

 

 

二、快速傅里叶变换

其实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获得的复数的模（即绝对值）就是对应的“振幅谱”，复数所对应的角度，就是所对应的“相位谱”，如今能够画图了。

 

三、FFT的原始频谱

N=1400
x = 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的一半，前半部分便可。

 

四、将振幅谱进行归一化和取半处理

先进行归一化

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()

这就是咱们最终的结果，须要的“振幅谱”。

 

3、完整代码

import numpy as np
from scipy.fftpack import fft,ifft
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个采样点，同样意思的）
x=np.linspace(0,1,1400)      
 
#设置须要采样的信号，频率份量有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)                          #快速傅里叶变换
 
N=1400
x = 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()

 

疑问解答，关于归一化和取一半处理需看快速傅里叶变换在信号处理中的应用，具体为：

　　傅里叶变换FT（Fourier Transform）是一种将信号从时域变换到频域的变换形式。它在声学、信号处理等领域有普遍的应用。计算机处理信号的要求是：在时域和频域都应该是离散的，并且都应该是有限长的。而傅里叶变换仅能处理连续信号，离散傅里叶变换DFT（Discrete Fourier Transform）就是应这种须要而诞生的。它是傅里叶变换在离散域的表示形式。可是通常来讲，DFT的运算量是很是大的。在1965年首次提出快速傅里叶变换算法FFT（Fast Fourier Transform）以前，其应用领域一直难以拓展，是FFT的提出使DFT的实现变得接近实时。DFT的应用领域也得以迅速拓展。除了一些速度要求很是高的场合以外，FFT算法基本上能够知足工业应用的要求。因为数字信号处理的其它运算均可以由DFT来实现，所以FFT算法是数字信号处理的重要基石。

　　傅立叶原理代表：任何连续测量的时序或信号，均可以表示为不一样频率的正弦波信号的无限叠加。而根据该原理创立的傅立叶变换算法利用直接测量到的原始信号，以累加方式来计算该信号中不一样正弦波信号的频率、振幅和相位。如图1所示，即为时域信号与不一样频率的正弦波信号的关系。图中最右侧展现的是时域中的一个信号，这是一个近似于矩形的波，而图的正中间则是组成该信号的各个频率的正弦波。从图中咱们能够看出，即便角度几乎为直角的正弦波，其实也是由众多的弧度圆滑的正弦波来组成的。在时域图像中，咱们看到的只有一个矩形波，咱们无从得知他是由这些正弦波组成。但当咱们经过傅里叶变换将该矩形波转换到频域以后，咱们可以很清楚的看到许多脉冲，其中频域图中的横轴为频率，纵轴为振幅。所以能够经过这个频域图像得知，时域中的矩形波是由这么多频率的正弦波叠加而成的。

　　这就是傅里叶变换的最基本最简单的应用，固然这是从数学的角度去看傅立叶变换。在信号分析过程当中，傅里叶变换的做用就是将组成这个回波信号的全部输入源在频域中按照频率的大小来表示出来。傅里叶变换以后，信号的幅度谱可表示对应频率的能量，而相位谱可表示对应频率的相位特征。通过傅立叶变换能够在频率中很容易的找出杂乱信号中各频率份量的幅度谱和相位谱，而后根据需求，进行高通或者低通滤波处理，最终获得所须要频率域的回波。

　　傅里叶变换在图像处理过程当中也有很是重要的做用，设信号f是一个能量有限的模拟信号，则其傅里叶变换就表示信号f的频谱。从纯粹的数学意义上看，傅里叶变换是将一个函数转换为一系列周期函数来处理的。从物理效果看，傅里叶变换是将图像从空间域转换到频率域，其逆变换是将图像从频率域转换到空间域。换句话说，傅里叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数。傅里叶逆变换是将图像的频率分布函数变换为灰度分布函数。傅里叶频谱图上咱们看到的明暗不一的亮点，其意义是指图像上某一点与邻域点差别的强弱，即梯度的大小，也即该点的频率的大小。通常来说，梯度大则该点的亮度强，不然该点亮度弱。这样经过观察傅里叶变换后的频谱图，也叫功率图，咱们就能够直观地看出图像的能量分布：若是频谱图中暗的点数更多，那么实际图像是比较柔和的，这是由于各点与邻域差别都不大，梯度相对较小；反之，若是频谱图中亮的点数多，那么实际图像必定是尖锐的、边界分明且边界两边像素差别较大的。

　　以信号处理过程当中的一个例子来详细说明FFT的效果：假设采样频率为Fs，信号频率为F，采样点数为N。那么FFT处理以后的结果就是一个点数为N点的复数。每个点就对应着一个频率点，而每一个点的模值，就是该频率值下的幅度特性。假设原始信号的峰值为A，那么在处理后除第一个点以外的其余点的模值就是A的N/2倍。而第一个点就是直流份量，它的模值就是直流份量的N倍。而每一个点的相位呢，就是在该频率下的信号的相位。第一个点表示直流份量(即频率为0Hz)，而最后一个点N的再下一个点(实际上这个点是不存在的，这里是假设的第N+1个点，也能够看作是将第一个点分作两半分，另外一半移到最后)则表示采样频率Fs，这中间被N-1个点平均分红N等份，每一个点的频率依次增长。例如某点n所表示的频率为：Fn=(n-1)*Fs/N。由上面的公式能够看出，Fn 所能分辨到的最小频率为Fs/N，若是采样频率Fs为1024Hz，采样点数N为1024点，则最小分辨率能够精确到1Hz。1024Hz的采样率采样1024点，恰好是1秒，也就是说，采样1秒时间的信号并作FFT处理，则结果能够分析到1Hz；若是采样2秒时间的信号并作FFT处理，则结果能够精确到0.5Hz。

　　假设如今咱们有一个输入信号，该信号总共包含3种成分信号，其一是5V的直流份量；其二是频率为50Hz、相位为-60度、幅度为10V的交流信号；第三个成分信号是频率为100Hz、相位为90度、幅度为5V的交流信号。该输入信号用数学表达式表示以下：

 　　　　S=5+10*cos(2*pi*50*t-pi*60/180)+5*cos(2*pi*100*t+pi*90/180)

　　图2即为S信号的图像表示。如今，咱们以256Hz的采样率Fs对这个信号进行采样，采样点数N一样为256点。根据公式咱们能够算出其频谱图中的频率精度为1Hz。所以对于输入信号频率包含0Hz、50Hz和100hz的复合信号，在其通过FFT处理以后，应该会在频谱图中出现3个峰值，并且频率分别为0Hz、50Hz和100Hz，处理结果如图3所示：

　　结果正如咱们所预料的，对输入信号’S’作FFT处理以后，图3中出现了5个峰值，这是由于对输入信号作256点的FFT处理以后并无第257个频点信息，这也是前文中所提到的第一个点的模值是N倍的缘由。所以，信号的 FFT结果具备必定的对称性。通常状况下，咱们只使用前半部分的结果，即小于采样频率一半的结果。对于图像进行简单处理后，咱们的前半部分的FFT结果如图4所示：

　　从图4中能够看出，三个输入信号频点的幅值依次为1280、1280、640；其余频率所对应的幅值均为0。按照公式，能够计算直流份量（频率为0Hz）的幅值为：1280/N= 1280/256=5;频率为50Hz的交流信号的幅值为：1280/(N/2)= 1280/(256/2)=10;而75Hz的交流信号的幅值为640/(N/2)=640/(256/2)=5。这也正是咱们输入信号中的三个份量的直流份量值，因而可知，从频谱分析出来的幅值是正确的。

　　经过上面的例子能够看出，对于一个输入信号，假如咱们不能肯定该输入信号的频率组成，咱们对其进行FFT处理以后，便能够很轻松的看出其频率份量，而且能够经过简单的计算来获知该信号的幅值信息等。另外，若是想要提升频率分辨率，咱们根据计算公式首先想到的就是须要增长采样点数，但增长采样点数也就意味着计算量增长，这在工程应用中增长了工程难度。解决这个问题的方法有频率细分法，比较简单的方法是采样较短期的信号，而后在后面补充必定数量的0，使其长度达到须要的点数（通常为2的幂次方的点数），而后再作FFT，就能在必定程度上提升频率分辨率。