zoukankan      html  css  js  c++  java
  • FFT快速傅里叶变换的python实现

      FFT是DFT的高效算法,能够将时域信号转化到频域上,下面记录下一段用python实现的FFT代码。

    # encoding=utf-8
    
    import numpy as np
    import pylab as pl  # 导入和matplotlib同时安装的作图库pylab
    
    
    sampling_rate = 8000  # 采样频率8000Hz
    fft_size = 512   # 采样点512,就是说以8000Hz的速度采512个点,我们获得的数据只有这512个点的对应时刻和此时的信号值。
    t = np.linspace(0, 1, sampling_rate)   # 截取一段时间,截取是任意的,这里取了0~1秒的一段时间。
    
    x = np.sin(2*np.pi*156.25*t) + 2*np.sin(2*np.pi*234.375*t)  # 输入信号序列,人工生成了一段信号序列,范围在0~1秒
    xs = x[:fft_size]   # 由上所述,我们只采样了512个点,所以我们只获得了前512个点的数据
    xf = np.fft.rfft(xs)/fft_size  # 调用np.fft的函数rfft(用于实值信号fft),产生长度为fft_size/2+1的一个复数向量,分别表示从0Hz~4000Hz的部分,这里之所以是4000Hz是因为Nyquist定理,采样频率8000Hz,则能恢复带宽为4000Hz的信号。最后/fft_size是为了正确显示波形能量
    
    freqs = np.linspace(0, sampling_rate//2, fft_size//2 + 1)  # 由上可知,我们得到了数据,现在产生0~4000Hz的频率向量,方便作图
    xfp = 20*np.log10(np.clip(np.abs(xf), 1e-20, 1e1000))  # 防止幅值为0,先利用clip剪裁幅度,再化成分贝
    
    pl.figure(figsize=(8, 4))  # 生成画布
    pl.subplot(211)  # 生成子图,211的意思是将画布分成两行一列,自己居上面。
    pl.plot(t[:fft_size], xs)  # 对真实波形绘图
    pl.xlabel(u"time(s)")
    pl.title(u"The Wave and Spectrum of 156.25Hz and 234.375Hz")
    pl.subplot(212)  # 同理
    pl.plot(freqs, xfp)  # 对频率和幅值作图,xlabel是频率Hz,ylabel是dB
    pl.xlabel(u"Hz")
    pl.subplots_adjust(hspace=0.4)  # 调节绘图参数
    pl.show()

      代码进行了详细标注。有一个小细节是FFT对于取样时间有要求。N点FFT进行精确频谱分析的要求是N个取样点包含整数个取样对象的波形。因此N点FFT能够完美计算频谱,对取样对象的要求是n*Fs/N(n*采样频率/FFT长度)在本例中Fs = 8000Hz,N=512  base_freq=15.625Hz 所以本例中给出了频率为156.25Hz(n=10)和234.375Hz(n=15)做例子。

      效果如下:

  • 相关阅读:
    zabbix监控docker
    Ubuntu下Zabbix结合percona监控mysql数据
    centos7安装ftp
    Ubuntu 16.04 搭建 ELK
    ubuntu网卡配置及安装ssh服务
    CentOS7.5二进制安装MySQL-5.6.40
    生产环境MySQL数据库集群MHA上线实施方案
    Mysql主从复制
    GIt+jenkins代码自动上线
    虚拟机网卡丢失解决方法
  • 原文地址:https://www.cnblogs.com/chester-cs/p/11630779.html
Copyright © 2011-2022 走看看