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

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

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

一、一些关键概念的引入

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 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组数据,如下:

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=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的一半,前半部分即可。

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

先进行归一化

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 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,就能在一定程度上提高频率分辨率。

  

原文地址:https://www.cnblogs.com/tianqizhi/p/10850377.html

时间: 2024-11-02 15:57:52

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

研究傅里叶变换的一本好书<<快速傅里叶变换及其C程序>>

快速傅里叶变换及其C程序 <快速傅里叶变换及其C程序>是中国科学技术大学出版社出版的.本书系统地介绍了傅里叶变换的理论和技术,内容包括傅里叶变换(FT)的定义.存在条件及其性质,离散傅里叶变换(DFT)的定义.性质及由离散引起的频谱混叠和渗漏,快速傅里叶变换(FFT)算法的基本原理和复序列基2算法及其实用程序,并以此为基础,给出了实序列DFT.正弦变换.余弦变换.傅里叶级数.谱函数近似.功率谱估计.卷积和相关等的快速算法和实用程序,给出了 2D—DFT的行列算法.二维实序列2D—DFT的行列算

FFT算法实现——基于GPU的基2快速傅里叶变换

最近做一个东西,要用到快速傅里叶变换,抱着蛋疼的心态,自己尝试写了一下,遇到一些问题. 首先看一下什么叫做快速傅里叶变换(FFT)(来自Wiki): 快速傅里叶变换(英语:Fast Fourier Transform, FFT),是离散傅里叶变换的快速算法,也可用于计算离散傅里叶变换的逆变换.快速傅里叶变换有广泛的应用,如数字信号处理.计算大整数乘法.求解偏微分方程等等. 对于复数串行,离散傅里叶变换公式为: 直接变换的计算复杂度是O(n^2).快速傅里叶变换可以计算出与直接计算相同的结果,但只

快速傅里叶变换(FFT)算法【详解】

快速傅里叶变换(Fast Fourier Transform)是信号处理与数据分析领域里最重要的算法之一.我打开一本老旧的算法书,欣赏了JW Cooley 和 John Tukey 在1965年的文章中,以看似简单的计算技巧来讲解这个东西. 本文的目标是,深入Cooley-Tukey  FFT 算法,解释作为其根源的“对称性”,并以一些直观的python代码将其理论转变为实际.我希望这次研究能对这个算法的背景原理有更全面的认识. FFT(快速傅里叶变换)本身就是离散傅里叶变换(Discrete

快速傅里叶变换(FFT)

快速傅里叶变换(FFT)算法[详解] 快速傅里叶变换(Fast Fourier Transform)是信号处理与数据分析领域里最重要的算法之一.我打开一本老旧的算法书,欣赏了JW Cooley 和 John Tukey 在1965年的文章中,以看似简单的计算技巧来讲解这个东西. 本文的目标是,深入Cooley-Tukey  FFT 算法,解释作为其根源的"对称性",并以一些直观的python代码将其理论转变为实际.我希望这次研究能对这个算法的背景原理有更全面的认识. FFT(快速傅里叶

【笔记篇】(理论向)快速傅里叶变换(FFT)学习笔记w

现在真是一碰电脑就很颓废啊... 于是早晨把电脑锁上然后在旁边啃了一节课多的算导, 把FFT的基本原理整明白了.. 但是我并不觉得自己能讲明白... Fast Fourier Transformation, 快速傅里叶变换, 是DFT(Discrete Fourier Transform, 离散傅里叶变换)的快速实现版本. 据说在信号处理领域广泛的应用, 而且在OI中也有广泛的应用(比如SDOI2017 R2至少考了两道), 所以有必要学习一波.. 划重点: 其实学习FFT最好的教材是<算法导论

准零基础搞懂FFT快速傅里叶变换及其实现程序(二)

上一篇文章我们了解了DFT的原理,FFT是基于DFT的更适合计算机运算的算法,本文我们就正式开始学习FFT的原理. 首先我么先来宏观的看一下FFT.如果我们把整个FFT的算法看成一个黑盒子的话,那么它的输入就是时间波形信号,比如声音波形(横轴为时间,纵轴为振幅).外什么FFT要比DFT速度更快呢?下面(图1)解释了FFT和DFT的(对于计算机的)算法复杂度 图1 从上面的数学表达式可以看出,一个1024采样点的FFT比DFT块了102.4倍.如果傅里叶变换的数量级更大,FFT的速度优势会更明显.

快速傅里叶变换中的位码倒置算法

最近一直在看傅里叶变换,看到FFT算法,其实算法的关键之一,蝶形运算,只要看懂了,编码实现并不难.反倒是其中位码倒序的环节,看很容易看懂,但是编码实现不是那么容易的.在网上参考了很多资料后,决定把下面这个算法分享给大家,在这里要感谢百度文库用户letsgotoyy123提供的<快速傅里叶变换FFT及其应用>一文(http://wenku.baidu.com/view/7e9f6e6ea1c7aa00b42acb84)还要感谢星星同学,在一些基础知识上的指点,谢谢. #include<io

【转】快速傅里叶变换(FFT)详解

目录 前言 多项式 系数表示法 点值表示法 复数 向量 圆的弧度制 平行四边形定则 复数 运算法则 单位根 单位根的性质 快速傅里叶变换 快速傅里叶逆变换 理论总结 递归实现 迭代实现 本文只讨论FFT在信息学奥赛中的应用 文中内容均为个人理解,如有错误请指出,不胜感激 回到顶部 前言 先解释几个比较容易混淆的缩写吧 DFT:离散傅里叶变换->O(n2)O(n2)计算多项式乘法 FFT:快速傅里叶变换->O(n?log(n)O(n?log?(n)计算多项式乘法 FNTT/NTT:快速傅里叶变换

短时傅里叶变换(Short Time Fourier Transform)原理及 Python 实现

原理 短时傅里叶变换(Short Time Fourier Transform, STFT) 是一个用于语音信号处理的通用工具.它定义了一个非常有用的时间和频率分布类, 其指定了任意信号随时间和频率变化的复数幅度. 实际上,计算短时傅里叶变换的过程是把一个较长的时间信号分成相同长度的更短的段, 在每个更短的段上计算傅里叶变换, 即傅里叶频谱. 短时傅里叶变换通常的数学定义如下: 其中, DTFT (Decrete Time Fourier Transform) 为离散时间傅里叶变换.  其数学公