jjzjj

Python中利用FFT(快速傅里叶变换)进行频谱分析

Mobius8086 2023-12-22 原文

本文将从实例的角度出发讲解fft函数的基本使用,不包含复杂的理论推导。

一、基本条件

要对一个信号进行频谱分析,首先需要知道几个基本条件。

  1. 采样频率fs
  2. 信号长度N(信号的点数)

采样频率fs:根据采样定理可知,采样频率应当大于等于被测信号里最高频率的2倍,才能保证不失真,但是实际情况下,我们可能并不知道最高频率是多少,所以这个就是根据一定的经验或者搜索得到的,比如本次所使用到的ECG(心电)信号,最高频率一般不高于100Hz,于是我们设置采样频率为500(原本200Hz就够了,但是实际工程一般会将标准放大3~5倍,以便留有一定的裕量)。

信号长度N:这个一般很容易获得,因为我们经过采样得到的信号都是离散信号,如果是一维的,只需要使用len函数就可以直接获得信号的点数。

二、设计代码实例

在Python的第三方库文件中,numpy和scipy都有fft的函数,本文使用scipy中的fft函数。

from scipy.fftpack import fft

打开fft函数的源文件,可以看到如下内容:

def fft(x, n=None, axis=-1, overwrite_x=False):
    """
    Return discrete Fourier transform of real or complex sequence.

    The returned complex array contains ``y(0), y(1),..., y(n-1)``, where

    ``y(j) = (x * exp(-2*pi*sqrt(-1)*j*np.arange(n)/n)).sum()``.

    Parameters
    ----------
    x : array_like
        Array to Fourier transform.
    n : int, optional
        Length of the Fourier transform. If ``n < x.shape[axis]``, `x` is
        truncated. If ``n > x.shape[axis]``, `x` is zero-padded. The
        default results in ``n = x.shape[axis]``.

阅读一下可以得到输入参数的信息:

x:待进行FFT的信号序列。

n:可选,傅里叶变换的点数,如果n的长度小于序列x长度,则x将会被截断(truncated),如果n大于序列长度,则序列将会补零,默认是两者相等。

于是,我们开始对fft函数根据自己的需要进行函数封装编写:

def FFT(Fs, data):
    """
    对输入信号进行FFT
    :param Fs:  采样频率
    :param data:待FFT的序列
    :return:
    """
    L = len(data)  # 信号长度
    N = np.power(2, np.ceil(np.log2(L)))  # 下一个最近二次幂,也即N个点的FFT
    result = np.abs(fft(x=data, n=int(N))) / L * 2  # N点FFT
    axisFreq = np.arange(int(N / 2)) * Fs / N  # 频率坐标
    result = result[range(int(N / 2))]  # 因为图形对称,所以取一半
    return axisFreq, result

代码解读:

L是输入序列的长度,很容易理解并获得。

N为大于序列长度的第一个2的幂次数,比如3之后第一个2的幂次数为,5之后的第一个2的幂次数为,因为FFT变换一般选择2的幂次进行,但是计算机计算的话,我们也可以不必如此费事,反正不是我们自己手算。

傅里叶变换的核心代码在第三行,abs是取模运算,因为FFT返回的数据是复数,为了便于绘图,需要求模进行,至于为什么后面又除以L又乘以2,是为了保证变换前后的幅值和变换前一致,比如的信号,如果不进行这一步的话,得到的FFT结果在频率为1的那个“柱子”幅度就不是2,而是一个很大的数,会发生变化,至于为什么,这里理论不去深究。

第四行代码时获取绘图所用的频率x轴,Fs/N是频率分辨率,乘以N个点就可以获得一个序列,这个序列就是每个频率分布点的x轴坐标点,N/2是为了截取一半,我们知道FFT绘制出来的图形是左右对称的,因此这里只取前一半,下面一行代码也是一样,取前一半。

举例说明:通过手动设计被测信号和采样频率等条件,掌握对库函数fft的使用和理解,以便迁移到实际工程中。下面这个例子中,我们设计了一个信号频率为390和2000Hz的正弦叠加信号,该信号的样子如图1所示,对其进行FFT之后如图2所示。

import numpy as np
from scipy.fftpack import fft
import matplotlib.pyplot as plt


def FFT(Fs, data):
    """
    对输入信号进行FFT
    :param Fs:  采样频率
    :param data:待FFT的序列
    :return:
    """
    L = len(data)  # 信号长度
    N = np.power(2, np.ceil(np.log2(L)))  # 下一个最近二次幂,也即N个点的FFT
    result = np.abs(fft(x=data, n=int(N))) / L * 2  # N点FFT
    axisFreq = np.arange(int(N / 2)) * Fs / N  # 频率坐标
    result = result[range(int(N / 2))]  # 因为图形对称,所以取一半
    return axisFreq, result


if __name__ == '__main__':
    Fs = 10000  # 采样频率
    f1 = 390  # 信号频率1
    f2 = 2000  # 信号频率2
    t = np.linspace(0, 1, Fs)  # 生成 1s 的时间序列
    # 给定信号
    y = 2 * np.sin(2 * np.pi * f1 * t) + 5 * np.sin(2 * np.pi * f2 * t)
    # 第一步,对没有添加噪声的信号进行FFT,验证分析是否正确
    x, result = FFT(Fs, y)

    # 绘图
    fig1 = plt.figure(figsize=(16, 9))
    plt.title('original data')
    plt.plot(t, y)
    plt.xlabel('time/s')
    plt.ylabel('Amplitude')
    plt.grid()

    fig2 = plt.figure(figsize=(16, 9))
    plt.title('FFT')
    plt.plot(x, result)
    plt.xlabel('Frequency/Hz')
    plt.ylabel('Amplitude')
    plt.grid()
    plt.show()

输出结果如下:

图1. 原始信号
图2. 对信号进行FFT

结果分析:从FFT的结果我们可以看到,在频率为390到2000Hz的地方有两个凸起的“高峰”,说明原信号在此处的频率占据比重较大,而且FFT的结果也满足原信号的幅值大小,即390Hz信号的幅度为2,2000Hz的信号幅度为5。


如果我们对信号进行加噪,可以查看FFT分析结果。

噪声信号选择随机信号,注意,点数要和被分析的信号的点数保持一致。

# 0-1 之间的随机噪声乘以10倍,也即0-10之间的噪声
noise1 = 10*np.random.random(10000)  
图3. 没有加入噪声的FFT频谱
图4.加入噪声之后的FFT频谱

 从图中可以看到的是,加入的噪声信号有直流成分(0Hz),其他横坐标上密密麻麻分布的是一些零碎的频率信号。

参考文章:

Python 中 FFT 快速傅里叶分析 - 知乎

有关Python中利用FFT(快速傅里叶变换)进行频谱分析的更多相关文章

  1. ruby-on-rails - 使用 Ruby on Rails 进行自动化测试 - 最佳实践 - 2

    很好奇,就使用ruby​​onrails自动化单元测试而言,你们正在做什么?您是否创建了一个脚本来在cron中运行rake作业并将结果邮寄给您?git中的预提交Hook?只是手动调用?我完全理解测试,但想知道在错误发生之前捕获错误的最佳实践是什么。让我们理所当然地认为测试本身是完美无缺的,并且可以正常工作。下一步是什么以确保他们在正确的时间将可能有害的结果传达给您? 最佳答案 不确定您到底想听什么,但是有几个级别的自动代码库控制:在处理某项功能时,您可以使用类似autotest的内容获得关于哪些有效,哪些无效的即时反馈。要确保您的提

  2. python - 如何使用 Ruby 或 Python 创建一系列高音调和低音调的蜂鸣声? - 2

    关闭。这个问题是opinion-based.它目前不接受答案。想要改进这个问题?更新问题,以便editingthispost可以用事实和引用来回答它.关闭4年前。Improvethisquestion我想在固定时间创建一系列低音和高音调的哔哔声。例如:在150毫秒时发出高音调的蜂鸣声在151毫秒时发出低音调的蜂鸣声200毫秒时发出低音调的蜂鸣声250毫秒的高音调蜂鸣声有没有办法在Ruby或Python中做到这一点?我真的不在乎输出编码是什么(.wav、.mp3、.ogg等等),但我确实想创建一个输出文件。

  3. ruby-on-rails - 按天对 Mongoid 对象进行分组 - 2

    在控制台中反复尝试之后,我想到了这种方法,可以按发生日期对类似activerecord的(Mongoid)对象进行分组。我不确定这是完成此任务的最佳方法,但它确实有效。有没有人有更好的建议,或者这是一个很好的方法?#eventsisanarrayofactiverecord-likeobjectsthatincludeatimeattributeevents.map{|event|#converteventsarrayintoanarrayofhasheswiththedayofthemonthandtheevent{:number=>event.time.day,:event=>ev

  4. ruby - 使用 C 扩展开发 ruby​​gem 时,如何使用 Rspec 在本地进行测试? - 2

    我正在编写一个包含C扩展的gem。通常当我写一个gem时,我会遵循TDD的过程,我会写一个失败的规范,然后处理代码直到它通过,等等......在“ext/mygem/mygem.c”中我的C扩展和在gemspec的“扩展”中配置的有效extconf.rb,如何运行我的规范并仍然加载我的C扩展?当我更改C代码时,我需要采取哪些步骤来重新编译代码?这可能是个愚蠢的问题,但是从我的gem的开发源代码树中输入“bundleinstall”不会构建任何native扩展。当我手动运行rubyext/mygem/extconf.rb时,我确实得到了一个Makefile(在整个项目的根目录中),然后当

  5. ruby - 如何进行排列以有效地定制输出 - 2

    这是一道面试题,我没有答对,但还是很好奇怎么解。你有N个人的大家庭,分别是1,2,3,...,N岁。你想给你的大家庭拍张照片。所有的家庭成员都排成一排。“我是家里的friend,建议家庭成员安排如下:”1岁的家庭成员坐在这一排的最左边。每两个坐在一起的家庭成员的年龄相差不得超过2岁。输入:整数N,1≤N≤55。输出:摄影师可以拍摄的照片数量。示例->输入:4,输出:4符合条件的数组:[1,2,3,4][1,2,4,3][1,3,2,4][1,3,4,2]另一个例子:输入:5输出:6符合条件的数组:[1,2,3,4,5][1,2,3,5,4][1,2,4,3,5][1,2,4,5,3][

  6. ruby - 即使失败也继续进行多主机测试 - 2

    我已经构建了一些serverspec代码来在多个主机上运行一组测试。问题是当任何测试失败时,测试会在当前主机停止。即使测试失败,我也希望它继续在所有主机上运行。Rakefile:namespace:specdotask:all=>hosts.map{|h|'spec:'+h.split('.')[0]}hosts.eachdo|host|begindesc"Runserverspecto#{host}"RSpec::Core::RakeTask.new(host)do|t|ENV['TARGET_HOST']=hostt.pattern="spec/cfengine3/*_spec.r

  7. ruby - 是否可以覆盖 gemfile 进行本地开发? - 2

    我们的git存储库中目前有一个Gemfile。但是,有一个gem我只在我的环境中本地使用(我的团队不使用它)。为了使用它,我必须将它添加到我们的Gemfile中,但每次我checkout到我们的master/dev主分支时,由于与跟踪的gemfile冲突,我必须删除它。我想要的是类似Gemfile.local的东西,它将继承从Gemfile导入的gems,但也允许在那里导入新的gems以供使用只有我的机器。此文件将在.gitignore中被忽略。这可能吗? 最佳答案 设置BUNDLE_GEMFILE环境变量:BUNDLE_GEMFI

  8. ruby - 在 Windows 机器上使用 Ruby 进行开发是否会适得其反? - 2

    这似乎非常适得其反,因为太多的gem会在window上破裂。我一直在处理很多mysql和ruby​​-mysqlgem问题(gem本身发生段错误,一个名为UnixSocket的类显然在Windows机器上不能正常工作,等等)。我只是在浪费时间吗?我应该转向不同的脚本语言吗? 最佳答案 我在Windows上使用Ruby的经验很少,但是当我开始使用Ruby时,我是在Windows上,我的总体印象是它不是Windows原生系统。因此,在主要使用Windows多年之后,开始使用Ruby促使我切换回原来的系统Unix,这次是Linux。Rub

  9. Python 相当于 Perl/Ruby ||= - 2

    这个问题在这里已经有了答案:关闭10年前。PossibleDuplicate:Pythonconditionalassignmentoperator对于这样一个简单的问题表示歉意,但是谷歌搜索||=并不是很有帮助;)Python中是否有与Ruby和Perl中的||=语句等效的语句?例如:foo="hey"foo||="what"#assignfooifit'sundefined#fooisstill"hey"bar||="yeah"#baris"yeah"另外,类似这样的东西的通用术语是什么?条件分配是我的第一个猜测,但Wikipediapage跟我想的不太一样。

  10. java - 什么相当于 ruby​​ 的 rack 或 python 的 Java wsgi? - 2

    什么是ruby​​的rack或python的Java的wsgi?还有一个路由库。 最佳答案 来自Python标准PEP333:Bycontrast,althoughJavahasjustasmanywebapplicationframeworksavailable,Java's"servlet"APImakesitpossibleforapplicationswrittenwithanyJavawebapplicationframeworktoruninanywebserverthatsupportstheservletAPI.ht

随机推荐