基于汉宁窗、多窗口和小波的时频分析

有态度 有深度 有温度

你的科研小伙伴

如何使用汉宁窗、多窗口和小波对单个对象的脑电信号进行时频分析？（本文以MEG数据为例）相信你看完这篇文章，会有所收获。

我们先来认识一下汉宁窗、多窗口和小波分别是什么：

（1）汉宁窗：其实汉宁窗的英文写法有两种：hann/hanningwindow。目前这两种表述都可以，而且MATLAB中也存在hann和hanning两个函数。汉宁窗是窗函数之一，是升余弦窗的一个特例。汉宁窗可以看作是3个矩形时间窗的频谱之和，或者说是3个sinc（t）型函数之和，而括号中的两项相对于第一个谱窗向左、右各移动了π/T，从而使旁瓣互相抵消，消去高频干扰和漏能，适用于非周期性的连续信号。

（2）多窗口：多窗口法（multitaper）通过平均加窗数据的周期图以降低频谱估计的方差，多窗口法使用了不同的窗口函数，但数据还是用的整段信号。多窗口法使用由一系列正交锥形窗（tapers）组成的离散长球序列（discrete prolate spheroidalsequence）来产生一系列的加窗数据，并计算这些加窗数据周期图的平均值，将其作为信号的频谱估计。除了正交性之外，这些锥形窗还具有最佳的时频集中特性。因此，多窗口法得到的频谱估计方差小、频率分辨率较高。

本图来源于书籍《脑电信号处理与特征提取》

（3）小波：我们往期文章完美通俗解读小波变换，终于懂了小波是什么中有关于小波的详细介绍，感兴趣的小伙伴可以点此链接进行阅读。

某种小波的示意图

了解完汉宁窗、多窗口和小波后，接下来看如何运用它们进行时频分析。总的分析步骤如下：

（1）使用ft_definetrial和ft_preprocessing将数据读入MATLAB中；

（2）使用ft_selectdata将试次与每个条件分开；

（3）使用ft_freqanalysis计算每个频率箱和每个时间箱的功率值；

（4）可视化结果。通过单通道（ft_singleplotTFR）或多通道（ft_multiplotTFR）创建时间-频率图，或者以指定的时间-频率间隔创建地形图（ft_topoplotTFR）来实现。

时频分析的步骤示意图

本文为大家详细介绍了4种时频分析步骤。标题分别为时频分析①，时频分析②，时频分析③，时频分析④。

先用ft_preprocessing读取数据。该脚本是从-1.0秒到2.0秒读取数据，为了减少在试次开始和结束时发生的边界效应，建议设置比原定计划稍大一些的时间间隔。在MATLAB的命令窗口中输入Load dataFIC加载已经预处理过的数据。选择一个条件进行时频分析：

cfg = [];
cfg.trials = data_all.trialinfo == 3;
dataFIC = ft_redefinetrial(cfg, data_all);

时频分析①

汉宁锥度，固定窗长

如何使用汉宁锥计算时频，选择一个固定的时间窗长，根据时间窗长(δ T)确定频率分辨率。频率分辨率（δ f）=1 /时间窗长（秒）（δ T）。因此，500ms的时间窗则等于2Hz的频率分辨率（1/0.5s = 2 Hz），这说明可以用2Hz、4Hz、6Hz等频率来计算功率。

load dataFIC

下面的代码是以500ms的时间窗为例。

cfg              = [];
cfg.output       = 'pow';
cfg.channel      = 'MEG';
cfg.method       = 'mtmconvol';
cfg.taper        = 'hanning';
cfg.foi          = 2:2:30;                           % analysis 2 to 30 Hz in steps of 2 Hz
cfg.t_ftimwin    = ones(length(cfg.foi),1).*0.5;   % length of time window = 0.5 sec
cfg.toi          = -0.5:0.05:1.5;                    % time window "slides" from -0.5 to 1.5 sec in steps of 0.05 sec (50 ms)
TFRhann = ft_freqanalysis(cfg, dataFIC);

不管用什么方法来计算TFR，输出的格式都是相同的。输出格式是一个具有以下字段的结构：

TFRhann =

        label: {149x1 cell}                % Channel names
       dimord: 'chan_freq_time'            % Dimensions contained in powspctrm, channels X frequencies X time
         freq: [2 4 6 8 10 12 14 16 18 20 22 24 26 28 30]  % Array of frequencies of interest (the elements of freq may be different from your cfg.foi input depending on your trial length)
         time: [1x41 double]               % Array of time points considered
    powspctrm: [149x15x41 double]          % 3-D matrix containing the power values
         elec: [1x1 struct]                % Electrode positions etc
         grad: [1x1 struct]                % Gradiometer positions etc
          cfg: [1x1 struct]                % Settings used in computing this frequency decomposition

TFRhann功率谱包含每个特定频率的原始功率值的时间进程。

可视化

为了将事件相关电位的变化可视化，将基线间隔标准化。

使用ft_multiplotTFR函数绘制所有传感器的TFR，在cfg结构中可以调整设置参数。例如：

cfg = [];
cfg.baseline     = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.zlim         = [-2.5e-27 2.5e-27];
cfg.showlabels   = 'yes';
cfg.layout       = 'CTF151_helmet.mat';
figure
ft_multiplotTFR(cfg, TFRhann);

注意，使用选项cfg基线和cfg基线类型对数据进行基线校正。基线校正也可以通过ft_freqbaseline脚本实现。

使用ft_singleplotTFR绘制单个channel的图，以传感器MRC15的TFR为例，脚本设置如下。

cfg = [];
cfg.baseline     = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.maskstyle    = 'saturation';
cfg.zlim         = [-2.5e-27 2.5e-27];
cfg.channel      = 'MRC15';
cfg.layout       = 'CTF151_helmet.mat';
figure
ft_singleplotTFR(cfg, TFRhann);

使用ft_singleplotTFR获得的单个传感器的时频图

从上图中可以看出，在刺激开始后的0.9到1.3s的时间间隔内，功率在15-20 Hz左右增加。使用ft_topoplotTFR来绘制beta功率增加的地形图。

cfg = [];
cfg.baseline     = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.xlim         = [0.9 1.3];
cfg.zlim         = [-1e-27 1e-27];
cfg.ylim         = [15 20];
cfg.marker       = 'on';
cfg.layout       = 'CTF151_helmet.mat';
cfg.colorbar     = 'yes';
figure
ft_topoplotTFR(cfg, TFRhann);

使用ft_topoplotTFR获得时频表征（15-20Hz，刺激后0.9-1.3s）的地形图

时频分析②

汉宁锥度，频率相关窗长

根据频率变化的时间窗计算TFR。通常，时间窗随着频率的增加而变短。这种方法的主要优点是，时间平滑随着频率的增加而减小，从而增加了对短期效应的敏感性。然而，时间分辨率的提高是以频率分辨率为代价的。接下来主要演示如何使用滑动窗口汉宁锥度的方法实现频率相关的时间窗分析。该方法与小波分析非常相似。用Morlet小波进行的小波分析不同于高斯形状的锥度（见时频分析④）。

最好方法是首先选择每个时间窗口的周期数，这对所有频率都是相同的。例如，如果每个窗口的周期是7，那么7Hz的时间窗口是1000ms（1/7x 7个周期）；10Hz和时间窗是700ms  (1/10 x 7周期)和20Hz的时间窗是350 ms (1/20 x 7周期)。然而，太细的频率分辨率不仅不会提供新的信息，反倒会造成冗余。

以下是一个时间窗周期为7的具体设置参数。这里只对一个传感器（MRC15）进行计算，但当然也可以扩展到所有传感器。

cfg              = [];
cfg.output       = 'pow';
cfg.channel      = 'MRC15';
cfg.method       = 'mtmconvol';
cfg.taper        = 'hanning';
cfg.foi          = 2:1:30;
cfg.t_ftimwin    = 7./cfg.foi;  % 7 cycles per time window
cfg.toi          = -0.5:0.05:1.5;
TFRhann7 = ft_freqanalysis(cfg, dataFIC);

使用ft_singleplotTFR绘制结果：

cfg              = [];
cfg.baseline     = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.maskstyle    = 'saturation';
cfg.zlim         = [-2e-27 2e-27];
cfg.channel      = 'MRC15';
cfg.interactive  = 'no';
cfg.layout       = 'CTF151_helmet.mat';
figure
ft_singleplotTFR(cfg, TFRhann7);

使用ft_singleplotTFR获得的MRC15的时频图

时频分析③

多窗口

为了更好地控制频率平滑，通常使用多窗口法。在一个给定的时间窗口中，更多的锥形窗将导致更强的平滑。振荡伽玛活动（30-100Hz）是相当宽的频带了，分析这一信号成分得益于多窗口。对于小于30Hz的信号，建议只使用单锥，例如上面所示的Hanning锥（注意，在下面的例子中，使用多锥来分析低频，因为这个数据集的伽玛波段不受影响）。

基于多锥的时频分析可以通过ft_freqanalysis来实现。该函数使用滑动时间窗，在给定的频率下计算其功率。在通过离散傅里叶变换计算功率之前，数据是“锥形的”。每个时间窗可以使用几个正交锥形。计算每个锥形数据段的功率，然后进行合并。参数设置如下：

（1）Cfg.foi，感兴趣的频率，从1Hz到30Hz，以2Hz为梯度。

（2）Cfg.toi，兴趣的时间间隔。决定了计算功率值的时间窗的中心时间。Cfg.toi= -0.5:0.05:1.5表示功率值从-0.5到1.5s，以50ms为梯度。

（3）Cfg.t_ftimwin为滑动时间窗的长度，单位为s（=tw）。选择cfg.t_ftimwin = 5. / cfg.foi，表示每个时间窗有5个周期。

（4）Cfg.tapsmofrq确定频率平滑的宽度，单位为Hz（=fw）。选择cfg.tapsmofrq = cfg.foi*0.4，即平滑会随着频率的增加而增加。

设置完这些参数后，呈现的图如下：

a)使用多窗口法的TFRs设置的特征图；b)由以上设置参数产生的时间-频率图示例。

cfg = [];
cfg.output     = 'pow';
cfg.channel    = 'MEG';
cfg.method     = 'mtmconvol';
cfg.foi        = 1:2:30;
cfg.t_ftimwin  = 5./cfg.foi;
cfg.tapsmofrq  = 0.4 *cfg.foi;
cfg.toi        = -0.5:0.05:1.5;
TFRmult = ft_freqanalysis(cfg, dataFIC);

绘制结果

cfg = [];
cfg.baseline     = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.zlim         = [-2e-27 2e-27];
cfg.showlabels   = 'yes';
cfg.layout       = 'CTF151_helmet.mat';
cfg.colorbar     = 'yes';
figure
ft_multiplotTFR(cfg, TFRmult)

使用多窗口法计算功率的时频图。

时频分析④

Morlet小波

用多窗口法计算TFR的另一种方法是使用Morlet小波。该方法等效于使用高斯形状的锥体计算时间窗随频率变化的TFR。这里有一个关键参数是cfg.width，它以周期数确定小波的宽度。给定频率F的频谱带宽等于F/width*2（因此，在30 Hz和宽度为7时，频谱带宽为30/7*2= 8.6 Hz），而小波持续时间等于width/F/pi（在这种情况下，7/30/pi = 0.074s = 74ms）。

利用Morlet小波计算TFR。

cfg = [];
cfg.channel    = 'MEG';
cfg.method     = 'wavelet';
cfg.width      = 7;
cfg.output     = 'pow';
cfg.foi        = 1:2:30;
cfg.toi        = -0.5:0.05:1.5;
TFRwave = ft_freqanalysis(cfg, dataFIC);

绘制结果

cfg = [];
cfg.baseline     = [-0.5 -0.1];
cfg.baselinetype = 'absolute';
cfg.zlim         = [-2e-25 2e-25];
cfg.showlabels   = 'yes';
cfg.layout       = 'CTF151_helmet.mat';
cfg.colorbar     = 'yes';
figure
ft_multiplotTFR(cfg, TFRwave)

用Morlet小波计算功率的时频图。

参考资料：

1.胡理.脑电信号处理与特征提取[M].科学出版社.2020.

2.Percival, D. B. , &  Walden, A. T. . (1993). Spectral analysis for physical applications (multitaper and conventional univariate techniques) || references. , 10.1017/CBO9780511622762, 546-561.

3.Tallon-Baudry, C. , &  Bertrand, O. . (1999). Oscillatory gamma activity in humans and its role in object representation. Trends in Cognitive Sciences.

4.Mitra, P. P. , &  Pesaran, B. . (1998). Analysis of dynamic brain imaging data. Biophysical Journal, 76(2), 691-708.

5.https://www.fieldtriptoolbox.org/tutorial/timefrequencyanalysis/

MC_Brain

茗创科技工作室

觉得有帮助，欢迎转发收藏或者点个在看哦~

听说点在看的人SCI接收率都提升了18%呢！