专栏首页脑机接口信号处理之频谱原理与python实现

信号处理之频谱原理与python实现

EEG信号是大脑神经元电活动的直接反应,包含着丰富的信息,但EEG信号幅值小,其中又混杂有噪声干扰,如何从EEG信号中抽取我们所感兴趣的信号是一个极为重要的问题。自1932年Dietch首先提出用傅里叶变换方法来分析EEG信号,该领域相继引入了频域分析、时域分析等脑电分析的经典方法。

频谱分析

下面是一组用于描述和解释信号属性的常用量(matlab的常见形式,python中的常见形式也类似):

x: 采样的数据;

n=length(x): 样本数量;

fs: 采样频率(每单位时间或空间的样本数)(单位常用:赫兹Hz);

dt=1/fs :每样本的时间或空间增量(如果是时间上的增量,则又称:采样间隔或采样步长,单位常用:s);

t=(0:n-1)/fs : 数据的时间或空间范围;

y=fft(x) : 数据的离散傅里叶变换(DFT);

abs(y) :DFT的振幅;

(abs(y).^2)/n :DFT的幂;

fs/n : 频率增量;

f=(0:n-1) * (fs/n) : 频率范围;

fs/2 :Nyquist频率(频率范围的中点);

频谱分析是一种将复噪声号分解为较简单信号的技术。真实世界中的信号可能由多种简单信号叠加而成。找出一个信号在不同频率下的信息(可能是幅度、功率、强度或相位等)的作法就是频谱分析。

采样定理:采样频率要大于信号频率的两倍。

N个采样点经过FFT变换后得到N个点的以复数形式记录的FFT结果。

假设采样频率为Fs,采样点数为N。那么FFT运算的结果就是N个复数(或N个点),每一个复数就对应着一个频率值以及该频率信号的幅值和相位。

第一个点对应的频率为0Hz(即直流分量),最后一个点N的下一个点对应采样频率Fs。其中任意一个采样点n所代表的信号频率:

这表明,频谱分析得到的信号频率最大为 (N-1)*Fs/N,对频率的分辨能力是Fs/N。采样频率和采样时间制约着通过FFT运算能分析得到的信号频率上限,同时也限定了分析得到的信号频率的分辨率。

每一个复数的模值对应该点所对应的频率值的幅度特性,具体的定量关系如下:

假设信号由以下周期的原始信号叠加而成:

那么,在经过FFT分析后得到的第一个点的模值是A1的N倍,而且只有在FFT结果点对应的频率在ω2,ω3时,其模值才明显放大,在其他频率点,模值接近于0。在这些模值明显放大的点中,除第一个点之外的其它点模值是相应信号幅值的N/2倍。

每个复数的相位就是在该频率值下信号的相位:φ2,φ3。

FFT结果有对称性,通常我们只是用前半部分的结果,也就是小于采样频率一半的结果。同时也只有采样频率一半以内、具有一定幅值的信号频率才是真正的信号频率。

下面就用python案例进行说明

案例1

import numpy as np
import pylab as pl
import math

# 采样频率
fs=1048
# 采样步长
t = [x/1048.0 for x in range(1048)]
"""
设计的采样值
假设信号y由4个周期信号叠加所得,如下所示
"""
y = [ 3.0 * np.cos(2.0 * np.pi * 50 * t0 - np.pi * 30/180)
     + 1.5 * np.cos(2.0 * np.pi * 75 * t0 + np.pi * 90/180)
     +  1.0 * np.cos(2.0 * np.pi * 150 * t0 + np.pi * 120/180)
     +  2.0 * np.cos(2.0 * np.pi * 220 * t0 + np.pi * 30/180)
     for t0 in t ]
pl.plot(t,y)
pl.xlabel('time(s)')
pl.title("original signal")
pl.show()
"""
现在对上述信号y在0-1秒时间内进行频谱分析,

本案例中采样频率为1048Hz,即单位时间内采样点数为1048
"""
# 采样点数
N=len(t)
# 采样频率
fs=1048.0
# 分辨率
df = fs/(N-1)
# 构建频率数组
f = [df*n for n in range(0,N)]
Y = np.fft.fft(y)*2/N  #*2/N 反映了FFT变换的结果与实际信号幅值之间的关系
absY = [np.abs(x) for x in Y]      #求傅里叶变换结果的模

pl.plot(f,absY)
pl.xlabel('freq(Hz)')
pl.title("fft")
pl.show()

案例2

from scipy.fftpack import fft, fftshift, ifft
from scipy.fftpack import fftfreq
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

"""
t_s:采样周期
t_start:起始时间
t_end:结束时间
"""
t_s = 0.01
t_start = 0.5
t_end = 5
t = np.arange(t_start, t_end, t_s)

f0 = 5
f1 = 20

# 绘制图表
plt.figure(figsize=(10, 12))

# 构建原始信号序列
y = 1.5*np.sin(2*np.pi*f0*t) + 3*np.sin(2*np.pi*20*t) + np.random.randn(t.size)
ax=plt.subplot(511)
ax.set_title('original signal')
plt.tight_layout()
plt.plot(y)

"""
FFT(Fast Fourier Transformation)快速傅里叶变换
"""
Y = fft(y)
ax=plt.subplot(512)
ax.set_title('fft transform')
plt.plot(np.abs(Y))

"""
Y = fftshift(X) 通过将零频分量移动到数组中心,重新排列傅里叶变换 X。
"""
shift_Y = fftshift(Y)
ax=plt.subplot(513)
ax.set_title('shift fft transform')
plt.plot(np.abs(shift_Y))

"""
得到正频率部分
"""
pos_Y_from_fft = Y[:Y.size//2]
ax=plt.subplot(514)
ax.set_title('fft transform')
plt.tight_layout()
plt.plot(np.abs(pos_Y_from_fft))

"""
直接截取 shift fft结果的前半部分
"""
pos_Y_from_shift = shift_Y[shift_Y.size//2:]
ax=plt.subplot(515)
ax.set_title('shift fft cut')
plt.plot(np.abs(pos_Y_from_shift))
plt.show()

代码来源于网络,本文对代码进行注释并整理

众所周知,傅里叶变换的快速算法FFT可以用来对信号的频域特征进行分析,然而,FFT仅能用于平稳信号的分析,对于非平稳信号,则需要采用短时傅里叶变换(STFT)进行分析。

对于非平稳信号,短时傅里叶变换所采用的策略是在信号上面加窗,一般是hamming窗,当然也可以是其它类型的窗函数,加窗之后的信号被分割为一组短长度子序列,子序列可以近似的看为平稳序列,可以用傅里叶变换的方式去进行分析,这也是短时傅里叶变换的精髓所在。

用STFT进行脑电信号分析一般有两种思路,第一种是用STFT来分离EEG信号的波段,从而求得每一个波段的能量作为特征(alpha、beta、theta、gamma、deta)。第二种是利用STFT计算功率谱密度作为特征,功率谱密度(PSD)特征可以针对整个信号子序列也可以针对子序列中特定的波段来计算。这两种思路中,第二种思路用的比较广,下面对其进行说明。

matlab中进行STFT的函数为spectrogram,计算功率谱密度(PSD)时使用如下格式:

[S,F,T,P]=spectrogram(x,window,noverlap,nfft,fs)

其中,S为输入信号x的短时傅里叶变换,F为频率向量,T为时间向量,P为功率谱密度矩阵,x为输入信号,window为时间窗,noverlap为overlap的点数,如果为0就是没有overlap,nfft为DFT的点数,fs为采样频率。其中,F向量的维度和P的行数一致,可以根据F向量来选取特定波段的PSD,还可以将alpha、beta、theta、gamma、deta这几个波段分别分为几个窄波段,提取窄带PSD。

参考:

https://blog.csdn.net/lbaihao/article/details/81280030

https://blog.csdn.net/weiwei9363/article/details/86545691

https://blog.csdn.net/qrlhl/article/details/72155395

本文分享自微信公众号 - 脑机接口社区(Brain_Computer),作者:Rose

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-06-25

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 信号处理之功率谱原理与python实现

    功率谱是功率谱密度函数的简称,它定义为单位频带内的信号功率。它表示了信号功率随着频率的变化情况,即信号功率在频域的分布状况。

    脑机接口社区
  • EMD算法原理与实现

    SSVEP信号中含有自发脑电和大量外界干扰信号,属于典型的非线性非平稳信号。传统的滤波方法通常不满足对非线性非平稳分析的条件,1998年黄鄂提出希尔伯特黄变换(...

    脑机接口社区
  • 脑波的分类

    人类的脑电图中脑波频率一般在0.5~30Hz,通常按照频率进行分类以表示各种成分。下面是国际上的分类标准。一般将比α波慢的δ波与θ波统称为慢波;而将比α波快的β...

    脑机接口社区
  • Python机器学习的练习八:异常检测和推荐系统

    在这篇文章中,将会涉及两个话题——异常检测和推荐系统,我们将使用高斯模型实现异常检测算法并且应用它检测网络上的故障服务器。我们还将看到如何使用协同过滤创建推荐系...

    AiTechYun
  • 吴恩达机器学习笔记-2

    逻辑回归 (Logistic Regression)是分类问题的一个代表算法,这是目前最流行使用最广泛的一种学习算法。

    happy123.me
  • 【5分钟玩转Lighthouse】Python绘制图表

    本文将讲解如何在Lighthouse等云服务器上通过display、Python、Matplotlib等工具查看和绘制各类图表。

    溪歪歪
  • python 对矩阵进行复制操作 np.repeat 与 np.tile区别

    python 对矩阵进行复制操作 np.repeat 与 np.tile区别 二者区别 二者执行的是均是复制操作; np.repeat:复制的是多维数组的...

    学到老
  • Python的numpy库使用

    plt.plot(x, y, color='green', linewidth=2)

    用户7886150
  • [Python3] NumPy基础

    一、创建数组二、数组操作类型1. 数组属性2. 数组索引:获取单个元素3. 切片4. 数组的变形5. 数组拼接和分裂

    用户7886150
  • Python | 100道测试题,带你玩转Numpy模块!

    Numpy是Python做数据分析所必须要掌握的基础库之一。以下为入门Numpy的100题小练习,原为github上的开源项目,由和鲸社区的小科翻译并整理(保留...

    咸鱼学Python

扫码关注云+社区

领取腾讯云代金券