官方公众号:52phm,与工业互联网人一起成长

2021-12-14 00:04:26    互联网    1109    当前专栏:数字信号处理    分类:算法开发    本站官网:www.52phm.cn   

公众号 ...

使用python(scipy和numpy)实现快速傅里叶变换(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-120.j   9.66210986e-05-0.04305756j   3.86508070e-04-0.08611996j 
   8.69732036e-04-0.12919206j    1.54641157e-03-0.17227871j]
'''

我们发现以下几个特点:

(1)变换之后的结果数据长度和原始采样信号是一样的

(2)每一个变换之后的值是一个复数,为abj的形式,那这个复数是什么意思呢?

我们知道,复数abj在坐标系中表示为(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=A1A2cos(2πω2φ2)A3cos(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()

版权声明:遵循 CC 4.0 BY-SA 版权协议
原文链接:https://blog.csdn.net/qq_27825451/article/details/88553441

免责声明:信息来源于互联网聚合整理或者作者原创,非原创文章标示"互联网",原创文章标示"作者昵称",本文仅代表文章作者的个人观点,且文中陈述文字和内容的时效性、真实性与本站无关,请读者仅作参考,如发现有害或侵权内容,请联系我们将在第一时间进行核实处理。


[推荐] 基于python的快速傅里叶变换FFT

2021-12-14 14:06:34    互联网    610    分类:算法开发    专栏:数字信号处理   


[推荐] window Pycharm及python安装详细教程

2021-12-15 20:40:08    互联网    642    分类:开发环境    专栏:下载安装   


[推荐] linux-ubuntu下pycharm下载安装教程(社区版)

2021-12-16 17:49:15    互联网    608    分类:开发环境    专栏:下载安装   



转发此文章到社区


进群

让志同道合读者学习交流



window系统Python安装教程(新手)

       第一次接触Python,可能是爬虫或者是信息AI开发的小朋友,都说Python 语言简单,那么多学一些总是有好处的,下面从一个完全不懂的Python 的小白来安装Python 等一系列工作的记录,并且遇到的问题也会写出,让完全不懂的小白也可上手安装,并且完成第一个Hello world代码。[Python 安装]进入Python的官方下载页面http://www.python.org/download/出现很

2021-12-15 14:47:52    互联网    1184    分类:开发环境    专栏:下载安装   


Linux 安装最新版本python3

新安装了Linux系统(CentOS 6),发现已安装的python版本是2.6. 在网上搜索研究之后总结了一下怎么在保留python2的同时安装最新版的python3。1. 查看 Python 的版本号:命令行输入以下命令就可以查看python版本:#python -V 或# python --version2. 下载3.x新版本可以访问python的官方网站查看最新的python版本以及下载...

2021-12-16 17:35:07    互联网    644    分类:开发环境    专栏:下载安装   


突变点检测:Magnitude of trend之Sen's slope(python)

# Sen's slopeimport numpy as npfrom pandas import Seriesfrom scipy.stats import normdef sens_slope_trend_detection(inputdata,conf_level=0.95): inputdata = Series(inputdata) n = inputda...

2021-12-21 11:05:07    互联网    421    分类:算法开发    专栏:工业异常检测   


突变点检测:Buishand U test突变点检测(python)

import numpy as npimport pandas as pddef Buishand_U_change_point_detection(inputdata): inputdata = np.array(inputdata) inputdata_mean = np.mean(inputdata) n = inputdata.shape[0] k...

2021-12-21 11:07:15    互联网    354    分类:算法开发    专栏:工业异常检测   


python简介和下载安装-小知带你入门python

### 1、python介绍Python由荷兰数学和计算机科学研究学会的吉多·范罗苏姆 于1990 年代初设计,作为一门叫做ABC语言的替代品。Python提供了高效的高级数据结构,还能简单有效地面向对象编程。Python语法和动态类型,以及解释型语言的本质,使它成为多数平台上写脚本和快速开发应用的编程语言, 随着版本的不断更新和语言新功能的添加,逐渐被用于独立的、大型项目的开发。

2022-01-16 11:22:32    博客笔记    799    分类:算法开发    专栏:python   


python编码和声明-小知带你入门python

python编码和声明-小知带你入门python</br>python编码和声明 1、window环境下,在完成hello world运行后,已经初步掌握了python运行的机制,下面将介绍python编码方式。一般在pycharm编辑器中新建.py文件后,可以看到首行出现一行注释,例如

2022-01-17 23:15:26    博客笔记    660    分类:算法开发    专栏:python   


python输入input输出print函数-小知带你入门python

python输入input输出print函数python输入输出语句分别对应着input()和print()函数,下面分别对这两个函数进行介绍和实操。1、输入语句input()新建一个tmp.py文件当我们输入字符串时,返回数据类型是字符串# -*- coding: utf-8 -*-name = input("请输入您的名字:")print(name, type(name))运行tmp.py后,会提示手动输入信息,回车后就可以得到输出结果:请输入您的名字:小知小知 &lt;c

2022-01-17 23:25:01    博客笔记    645    分类:算法开发    专栏:python   


python变量名命名规则-小知带你入门python

python变量名命名规则变量名是指在代码中自定义的变量参数,比如a,b,_h, h1等,都可以称为变量名。下面介绍变量名的命名规范:变量名第一个字符必须是字母表中字母或下划线 _ ;变量名可以由字母、数字和下划线组成;变量名对大小写敏感,即大小写变量名含义不同;变量名不能与python内置关键字(比如def、if、else等)相同。1、变量名命名错误例子def = "晓物智联" # def 是python关键字,不能被当成变量名5ab = "晓物智联" # 变量名第一个字符不能为数

2022-01-17 23:36:20    博客笔记    628    分类:算法开发    专栏:python   


python模块import导入-小知带你入门python

python模块import导入模块导入分为python内置模块(或者第三方模块)导入和自定义模块导入。模块导入方法可以总结为5大导入方法,下面以python内置模块角度来举例介绍5大导入方法。(1)import… 方法表示导入某个模块# 导入numpy模块import numpyarr = numpy.array([1, 2, 4])print(arr)# [1 2 4](2)import…as…方法import xx as yy,意思是把xx作为yy表示,相当于xx可以使用yy来

2022-01-17 23:59:41    博客笔记    1082    分类:算法开发    专栏:python   


pyquickdoc快速生成python函数接口文档

`pyquickdoc` 支持快速生成 `Python` 库的函数接口文档,并且支持导出 `html` 与 `markdown` 文件,因此根据需求就写了一个能够自动生成 python 函数库的接口文档的程序,支持中文和英文函数接口获取和下载,可导出接口文档为 html 或者 markdown 格式,下面将介绍 pyquickdoc 的使用方法。

2022-04-05 12:09:04    博客笔记    1621    专栏:pyquickdoc