python傅里叶变换FFT绘制频谱图

本文实例为大家分享了python傅里叶变换FFT绘制频谱图的具体代码,供大家参考,具体内容如下

频谱图的横轴表示的是 频率, 纵轴表示的是振幅

#coding=gbk

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt 

#依据快速傅里叶算法得到信号的频域
def test_fft():
 sampling_rate = 8000 #采样率
 fft_size = 8000  #FFT长度
 t = np.arange(0, 1.0, 1.0/sampling_rate)
 x = np.sin(2*np.pi*156.25*t) + 2*np.sin(2*np.pi*234.375*t)+ 3*np.sin(2*np.pi*200*t)
 xs = x[:fft_size]

 xf = np.fft.rfft(xs) / fft_size #返回fft_size/2+1 个频率

 freqs = np.linspace(0, sampling_rate/2, fft_size/2+1) #表示频率
 xfp = np.abs(xf) * 2 #代表信号的幅值,即振幅

 plt.figure(num='original', figsize=(15, 6))
 plt.plot(x[:100])

 plt.figure(figsize=(8,4))
 plt.subplot(211)
 plt.plot(t[:fft_size], xs)
 plt.xlabel(u"时间(秒)", fontproperties='FangSong')
 plt.title(u"156.25Hz和234.375Hz的波形和频谱", fontproperties='FangSong')

 plt.subplot(212)
 plt.plot(freqs, xfp)
 plt.xlabel(u"频率(Hz)", fontproperties='FangSong')
 plt.ylabel(u'幅值', fontproperties='FangSong')
 plt.subplots_adjust(hspace=0.4)
 plt.show()

test_fft()
# np.clip(a, a_min, a_max, out) 输出与a 的shape一样,大于等于a_min,小于等于a_max的数,即在 [a_min, a_max]之间的数
a = np.arange(10)
print(a)
print(a.shape)
# [0 1 2 3 4 5 6 7 8 9]
b = np.empty((10,))
np.clip(a, 3, 8, out=b)
print(b)
# [3. 3. 3. 3. 4. 5. 6. 7. 8. 8.]
c = np.clip(a, 4, 10)
print(c)
# [4 4 4 4 4 5 6 7 8 9]
#a_min, a_max也可以输入与a 相同shape的数组
d = np.arange(4)
d1 = np.clip(d, [-1, 1, -3, 2], 2)
print(d)
print(d1)
# [0 1 2 3] #原数组
# [0 1 2 2] 

print(np.log10(1000))

def test_fft():
# FFT变换是针对一组数值进行运算的,这组数的长度N必须是2的整数次幂,例如64, 128, 256等等; 数值可以是实数也可以是复数,
# 通常我们的时域信号都是实数,因此下面都以实数为例。我们可以把这一组实数想像成对某个连续信号按照一定取样周期进行取样而得来,
# 如果对这组N个实数值进行FFT变换,将得到一个有N个复数的数组,我们称此复数数组为频域信号,此复数数组符合如下规律:
#
# 下标为0和N/2的两个复数的虚数部分为0,
# 下标为i和N-i的两个复数共轭,也就是其虚数部分数值相同、符号相反。
 np.random.seed(66)
 X = np.random.rand(8)
 print(X)
#  [0.15428758 0.13369956 0.36268547 0.67910888 0.19445006 0.25121038
# 0.75841639 0.55761859]
 xf = np.fft.fft(X)
 print(xf)
#  [ 3.0914769 +0.j   -0.20916178+0.39291702j -0.77236422+0.85181752j
#  0.12883683-0.39854483j -0.15179792+0.j   0.12883683+0.39854483j
#  -0.77236422-0.85181752j -0.20916178-0.39291702j]
 #通过快速傅里叶变换的逆变换 ifft 还原成原来的值
 X1 = np.fft.ifft(xf)
 print(X1)
# [0.15428758+0.00000000e+00j 0.13369956-2.00387919e-16j
# 0.36268547+1.66533454e-16j 0.67910888+1.51815661e-16j
# 0.19445006+0.00000000e+00j 0.25121038-1.51815661e-16j
# 0.75841639-1.66533454e-16j 0.55761859+2.00387919e-16j] 

# 下面让我们来看看FFT变换之后的那些复数都代表什么意思。
#
# 首先下标为0的实数表示了时域信号中的直流成分的多少
# 下标为i的复数a+b*j表示时域信号中周期为N/i个取样值的正弦波和余弦波的成分的多少, 其中a表示cos波形的成分,b表示sin波形的成分
 X = np.ones(8)
 x2 = np.fft.fft(X) / len(X) # 为了计算各个成分的能量多少,需要将FFT的结果除以FFT的长度
 print(x2)
# [1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]

 X = np.arange(0, 2*np.pi, 2*np.pi/8)
 y = np.sin(X)
 x3 = np.fft.fft(y) /len(y)
 print(x3)
# [ 1.43029718e-18+0.00000000e+00j -4.44089210e-16-5.00000000e-01j # 只有下标为 1 的复数的虚部为-0.5,
# 1.53080850e-17-1.38777878e-17j 3.87727691e-17-1.11022302e-16j
# 2.91858728e-17+0.00000000e+00j 0.00000000e+00-1.11022302e-16j
# 1.53080850e-17+1.38777878e-17j 3.44084101e-16+5.00000000e-01j]
 output1 = np.fft.fft(np.cos(X) / len(X))
 print(output1)
# [-4.30636606e-17+0.00000000e+00j 5.00000000e-01-2.66538563e-16j #只有下标为1 的实部为 0.5
# 1.53080850e-17+0.00000000e+00j 5.55111512e-17+1.97149624e-16j
# 1.24474906e-17+0.00000000e+00j -1.11022302e-16+2.05306223e-16j
# 1.53080850e-17+0.00000000e+00j 5.00000000e-01-1.35917284e-16j] 

 #综合的例子
 X = np.arange(0, 2*np.pi, 2*np.pi/128)
 y = 0.3*np.cos(X) + 0.5*np.cos(2*X+np.pi/4) + 0.8*np.cos(3*X-np.pi/3)
 yf = np.fft.fft(y) / len(y)
 print(2*np.abs(yf[1]), np.rad2deg(np.angle(yf[1])))
#  0.30000000000000016 3.3130777931911615e-15   #计算出幅值和相位角
 print(2*np.abs(yf[2]), np.rad2deg(np.angle(yf[2])))
#  0.5000000000000002 44.999999999999986
 print(2*np.abs(yf[3]), np.rad2deg(np.angle(yf[3])))
#  0.7999999999999998 -60.00000000000007

# 周期为128/1.0点的余弦波的相位为0, 振幅为0.3
# 周期为64/2.0点的余弦波的相位为45度, 振幅为0.5
# 周期为128/3.0点的余弦波的相位为-60度,振幅为0.8
# test_fft()

#使用多个正玄波合成三角波
import pylab as pl
# 取FFT计算的结果freqs中的前n项进行合成,返回合成结果,计算loops个周期的波形
def fft_combine(freqs, n, loops=1):
 length = len(freqs) * loops
 data = np.zeros(length)
 index = loops * np.arange(0, length, 1.0) / length * (2 * np.pi)
 for k, p in enumerate(freqs[:n]):
  if k != 0: p *= 2 # 除去直流成分之外,其余的系数都*2
  data += np.real(p) * np.cos(k*index) # 余弦成分的系数为实数部
  data -= np.imag(p) * np.sin(k*index) # 正弦成分的系数为负的虚数部
 return index, data 

# 产生size点取样的三角波,其周期为1
def triangle_wave(size):
 x = np.arange(0, 1, 1.0/size)
 y = np.where(x<0.5, x, 0)
 y = np.where(x>=0.5, 1-x, y)
 return x, y

def test_show():
 fft_size = 256

 # 计算三角波和其FFT
 x, y = triangle_wave(fft_size)
 fy = np.fft.fft(y) / fft_size

 # 绘制三角波的FFT的前20项的振幅,由于不含下标为偶数的值均为0, 因此取
 # log之后无穷小,无法绘图,用np.clip函数设置数组值的上下限,保证绘图正确
 pl.figure()
 pl.plot(np.clip(20*np.log10(np.abs(fy[:20])), -120, 120), "o")
 pl.xlabel("frequency bin")
 pl.ylabel("power(dB)")
 pl.title("FFT result of triangle wave")

 # 绘制原始的三角波和用正弦波逐级合成的结果,使用取样点为x轴坐标
 pl.figure()
 pl.plot(y, label="original triangle", linewidth=2)
 for i in [0,1,3,5,7,9]:
  index, data = fft_combine(fy, i+1, 2) # 计算两个周期的合成波形
  pl.plot(data, label = "N=%s" % i)
 pl.legend()
 pl.title("partial Fourier series of triangle wave")
 pl.show()

# test_show()

以上就是本文的全部内容,希望对大家的学习有所帮助,也希望大家多多支持我们。

时间: 2019-07-16

Python实现快速傅里叶变换的方法(FFT)

本文介绍了Python实现快速傅里叶变换的方法(FFT),分享给大家,具体如下: 这里做一下记录,关于FFT就不做介绍了,直接贴上代码,有详细注释的了: import numpy as np from scipy.fftpack import fft,ifft import matplotlib.pyplot as plt import seaborn #采样点选择1400个,因为设置的信号频率分量最高为600赫兹,根据采样定理知采样频率要大于信号频率2倍,所以这里设置采样频率为1400赫兹(即

opencv python 傅里叶变换的使用

理论 傅立叶变换用于分析各种滤波器的频率特性,对于图像,2D离散傅里叶变换(DFT)用于找到频域.快速傅里叶变换(FFT)的快速算法用于计算DFT. 于一个正弦信号,x(t)=Asin(2πft),我们可以说 f 是信号的频率,如果它的频率域被接受,我们可以看到 f 的峰值.如果信号被采样来形成一个离散信号,我们得到相同的频率域,但是在[−π,π] or [0,2π]范围内是周期性的 (or [0,N] for N-point DFT). 可以将图像视为在两个方向上采样的信号.因此,在X和Y方向

Opencv+Python实现图像运动模糊和高斯模糊的示例

运动模糊:由于相机和物体之间的相对运动造成的模糊,又称为动态模糊 Opencv+Python实现运动模糊,主要用到的函数是cv2.filter2D(): # coding: utf-8 import numpy as np import cv2 def motion_blur(image, degree=12, angle=45): image = np.array(image) # 这里生成任意角度的运动模糊kernel的矩阵, degree越大,模糊程度越高 M = cv2.getRotat

OpenCV+python手势识别框架和实例讲解

基于OpenCV2.4.8和 python 2.7实现简单的手势识别. 以下为基本步骤 1.去除背景,提取手的轮廓 2. RGB->YUV,同时计算直方图 3.进行形态学滤波,提取感兴趣的区域 4.找到二值化的图像轮廓 5.找到最大的手型轮廓 6.找到手型轮廓的凸包 7.标记手指和手掌 8.把提取的特征点和手势字典中的进行比对,然后判断手势和形状 提取手的轮廓 cv2.findContours() 找到最大凸包cv2.convexHull(),然后找到手掌和手指的相对位置,定位手型的轮廓和关键点

Opencv+Python 色彩通道拆分及合并的示例

一.图像色彩通道拆分 import cv2 img1 = cv2.imread(r"D:\OpencvTest\example.jpg", cv2.IMREAD_COLOR) # 传入一张彩色图片 b, g, r = cv2.split(img1) cv2.imshow("exampleB", b) # 展示B通道图 cv2.imshow("exampleG", g) cv2.imshow("exampleR", r) B通道

opencv python如何实现图像二值化

这篇文章主要介绍了opencv python如何实现图像二值化,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友可以参考下 代码如下 import cv2 as cv import numpy as np import matplotlib.pyplot as plt # 二值图像就是将灰度图转化成黑白图,没有灰,在一个值之前为黑,之后为白 # 有全局和局部两种 # 在使用全局阈值时,我们就是随便给了一个数来做阈值,那我们怎么知道我们选取的这个数的好坏呢?答

OpenCV python sklearn随机超参数搜索的实现

本文介绍了OpenCV python sklearn随机超参数搜索的实现,分享给大家,具体如下: """ 房价预测数据集 使用sklearn执行超参数搜索 """ import matplotlib as mpl import matplotlib.pyplot as plt import numpy as np import sklearn import pandas as pd import os import sys import tens

opencv python Canny边缘提取实现过程解析

这篇文章主要介绍了opencv python Canny边缘提取实现过程解析,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友可以参考下 Canny是边缘提取算法,在1986年提出的是一个很好的边缘检测器Canny算法介绍 非最大信号抑制: 高低阈值连接: example import cv2 as cv import numpy as np # canny运算步骤:5步 # 1. 高斯模糊 - GaussianBlur # 2. 灰度转换 - cvtCol

opencv python图像梯度实例详解

这篇文章主要介绍了opencv python图像梯度实例详解,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友可以参考下 一阶导数与Soble算子 二阶导数与拉普拉斯算子 图像边缘: Soble算子: 二阶导数: 拉普拉斯算子: import cv2 as cv import numpy as np # 图像梯度(由x,y方向上的偏导数和偏移构成),有一阶导数(sobel算子)和二阶导数(Laplace算子) # 用于求解图像边缘,一阶的极大值,二阶的零点

opencv+python实现均值滤波

本文实例为大家分享了opencv+python实现均值滤波的具体代码,供大家参考,具体内容如下 原理 均值滤波其实就是对目标像素及周边像素取平均值后再填回目标像素来实现滤波目的的方法,当滤波核的大小是3×3 3\times 33×3时,则取其自身和周围8个像素值的均值来代替当前像素值. 均值滤波也可以看成滤波核的值均为 1 的滤波. 优点:算法简单,计算速度快: 缺点:降低噪声的同时使图像产生模糊,特别是景物的边缘和细节部分. 代码 import cv2 as cv import numpy a