研一脑电小白日记 之 频谱与时频分析代码学习2 Welch法频谱估计

代码来源及详情下载网址:https://github.com/zhangzg78/eegbook

demo_welch.m Welch法频谱估计:使用不同参数计算EEG信号的Welch谱,对比分段重叠对频谱平滑度的影响
【脑电信号处理与特征提取】P6-张治国:频谱分析和时频分析

代码分析

demo_welch.m代码:

%% load data and define parameters
clear all; clc;

load data_eeg.mat
% data_eeg.mat contains 2 variables
%   - x: the EEG signal (with eyes-closed)
%   - Fs: the sampling rate (Fs = 160Hz)

N = length(x); % the number of samples (N=480)
x = detrend(x); % remove the low-frequency trend from EEG

%% spectral estimation (the Welch's method)
nfft = 2^nextpow2(N); % the number of FFT points
[P_per, f] = periodogram(x,[],nfft,Fs); % periodogram is also estimated for comparison
% check the help file to learn how to specify parameters in "pwelch.m"
% three parameter settings are used below
[P_wel_1, f] = pwelch(x,Fs,Fs/2,nfft,Fs);
[P_wel_2, f] = pwelch(x,Fs,0,nfft,Fs);
[P_wel_3, f] = pwelch(x,Fs/2,0,nfft,Fs);

%% display spectral estimates
f_lim = f((f>0)&(f<=50)); % specify the frequency range to be shown

figure('units','normalized','position',[0.1    0.3    0.8    0.5])
subplot(1,2,1)
hold on; box on;
plot(f,10*log10(P_per),'k','linewidth',0.5)
plot(f,10*log10(P_wel_1),'r','linewidth',2)
xlabel('Frequency (Hz)'); ylabel('Power (dB)')
hl = legend('Periodogram','Welch''s method (M=160, D=80)');
set(hl,'box','off','location','southwest')
set(gca,'xlim',[min(f_lim),max(f_lim)])

subplot(1,2,2)
hold on; box on;
plot(f,10*log10(P_wel_1),'r','linewidth',2)
plot(f,10*log10(P_wel_2),'g','linewidth',1)
plot(f,10*log10(P_wel_3),'b','linewidth',1)
xlabel('Frequency (Hz)'); ylabel('Power (dB)')
hl = legend('Welch''s (M=160, D=80)','Welch''s (M=160, D=0)','Welch''s (M=80, D=0)');
set(hl,'box','off','location','southwest')
set(gca,'xlim',[min(f_lim),max(f_lim)])
1. 参数定义与频谱估计
%% spectral estimation (the Welch's method)
nfft = 2^nextpow2(N); % 计算FFT点数(不小于信号长度的最小2的幂)
[P_per, f] = periodogram(x,[],nfft,Fs); % 周期图法估计频谱(用于对比)

% 三种不同的Welch方法参数设置
[P_wel_1, f] = pwelch(x, Fs, Fs/2, nfft, Fs);    % 设置1:窗口长度=160,重叠=80
[P_wel_2, f] = pwelch(x, Fs, 0, nfft, Fs);       % 设置2:窗口长度=160,重叠=0
[P_wel_3, f] = pwelch(x, Fs/2, 0, nfft, Fs);     % 设置3:窗口长度=80,重叠=0
  • 关键参数

    • nfft:FFT点数,由2^nextpow2(N)计算(例如N=480,则nfft=512)。

    • pwelch参数说明

      • 输入顺序pwelch(x, window, noverlap, nfft, Fs)

      • 窗口长度(window):以样本点为单位,但用户代码中错误地将Fs(采样率160Hz)作为窗口长度,实际应为样本点数(如160样本对应1秒窗口)。

      • 重叠点数(noverlap):需小于窗口长度,但用户代码中的Fs/2(80)可能表示重叠80样本。

      • 正确参数示例

        window_length = 160; % 窗口长度(样本点数)
        overlap = 80;        % 重叠样本点数
        [P_wel_1, f] = pwelch(x, window_length, overlap, nfft, Fs);


2. 频谱可视化
%% display spectral estimates
f_lim = f((f>0)&(f<=50)); % 限定显示频率范围(0-50Hz)

figure('units','normalized','position',[0.1 0.3 0.8 0.5]) % 创建全屏图形窗口

% 左子图:周期图 vs Welch方法(设置1)
subplot(1,2,1)
hold on; box on;
plot(f, 10*log10(P_per), 'k', 'linewidth', 0.5) % 周期图(对数刻度)
plot(f, 10*log10(P_wel_1), 'r', 'linewidth', 2) % Welch设置1(对数刻度)
xlabel('Frequency (Hz)'); ylabel('Power (dB)')
hl = legend('Periodogram','Welch''s method (M=160, D=80)');
set(hl, 'box', 'off', 'location', 'southwest') % 图例设置
set(gca, 'xlim', [min(f_lim), max(f_lim)])     % 频率范围限制

% 右子图:不同Welch参数对比
subplot(1,2,2)
hold on; box on;
plot(f, 10*log10(P_wel_1), 'r', 'linewidth', 2) % 设置1
plot(f, 10*log10(P_wel_2), 'g', 'linewidth', 1) % 设置2
plot(f, 10*log10(P_wel_3), 'b', 'linewidth', 1) % 设置3
xlabel('Frequency (Hz)'); ylabel('Power (dB)')
hl = legend('Welch''s (M=160, D=80)','Welch''s (M=160, D=0)','Welch''s (M=80, D=0)');
set(hl, 'box', 'off', 'location', 'southwest')
set(gca, 'xlim', [min(f_lim), max(f_lim)])
  • 左子图功能

    • 黑线(周期图):方差较大,频谱波动明显,尤其在低频段。

    • 红线(Welch设置1):通过分段和重叠平均,显著降低方差,曲线更平滑。

  • 右子图功能

    • 对比不同参数的影响

      • M=160, D=80(红色):长窗口+50%重叠,平滑度高但分辨率较低。

      • M=160, D=0(绿色):长窗口+无重叠,平滑度次之,计算效率更高。

      • M=80, D=0(蓝色):短窗口+无重叠,分辨率高但方差更大。

总结

该代码通过对比周期图法和不同参数下的Welch方法,展示了频谱估计中平滑度与分辨率的权衡

  • 左子图表明Welch方法(红色)显著优于周期图(黑色),但需注意参数设置的合理性。

  • 右子图揭示了窗口长度和重叠对频谱平滑度的影响,为实际应用中的参数选择提供了直观参考。

  • 改进方向:修正参数单位、优化窗函数、添加生理频段标注。

运行结果

结果分析

这段代码的运行结果图展示了两种频谱估计方法(周期图法和 Welch 方法)的比较,以及 Welch 方法在不同参数设置下的结果。以下是对运行结果图的详细分析:

左子图:周期图法与 Welch 方法的比较

  • 横轴:频率(Frequency,单位为 Hz)。

  • 纵轴:功率(Power,单位为 dB)。

  • 曲线说明

    • 黑色曲线:周期图法(Periodogram)的结果。

      • 周期图法的频谱估计结果波动较大,尤其是在高频段(如 20Hz 以上),这表明其方差较高。

    • 红色曲线:Welch 方法(M=160, D=80)的结果。

      • Welch 方法通过分段和平均操作,显著降低了频谱估计的方差,曲线更加平滑。

      • 在低频段(如 0-10Hz)和中频段(10-20Hz),Welch 方法的结果与周期图法的结果较为接近,但在高频段(20-50Hz),Welch 方法的结果更加稳定。

右子图:Welch 方法在不同参数设置下的比较

  • 横轴:频率(Frequency,单位为 Hz)。

  • 纵轴:功率(Power,单位为 dB)。

  • 曲线说明

    • 红色曲线:Welch 方法(M=160, D=80)的结果。

      • 窗口长度为 160,重叠长度为 80,频谱估计结果较为平滑,方差较低。

    • 绿色曲线:Welch 方法(M=160, D=0)的结果。

      • 窗口长度为 160,无重叠(D=0),频谱估计结果的波动较大,尤其是在高频段(20-50Hz),这表明无重叠时 Welch 方法的方差较高。

    • 蓝色曲线:Welch 方法(M=80, D=0)的结果。

      • 窗口长度为 80,无重叠(D=0),频谱估计结果的波动更大,尤其是在高频段(20-50Hz),这表明窗口长度较短且无重叠时,Welch 方法的方差进一步增加。

关键观察

  1. 周期图法 vs. Welch 方法

    • 周期图法的频谱估计结果波动较大,尤其是在高频段,这表明其方差较高。

    • Welch 方法通过分段和平均操作,显著降低了频谱估计的方差,结果更加平滑和稳定。

  2. Welch 方法的参数影响

    • 窗口长度(M

      • 较长的窗口长度(如 M=160)可以提供更好的频率分辨率,但计算量较大。

      • 较短的窗口长度(如 M=80)会降低频率分辨率,但计算量较小。

    • 重叠长度(D

      • 较大的重叠长度(如 D=80)可以进一步降低频谱估计的方差,结果更加平滑。

      • 无重叠(D=0)时,Welch 方法的方差较高,尤其是在高频段。

  3. 频率范围

    • 在低频段(0-10Hz),三种 Welch 方法的结果较为接近,且与周期图法的结果也较为接近。

    • 在中频段(10-20Hz),Welch 方法的结果开始显示出平滑的优势。

    • 在高频段(20-50Hz),Welch 方法的结果显著优于周期图法,尤其是当窗口长度较长且重叠长度较大时。

总结

  • Welch 方法的优势

    • 通过分段和平均操作,显著降低了频谱估计的方差,结果更加平滑和稳定。

    • 适用于需要高精度频谱估计的场景,尤其是在高频段。

  • 参数选择的影响

    • 较长的窗口长度和较大的重叠长度可以进一步降低方差,但会增加计算量。

    • 无重叠时,Welch 方法的方差较高,尤其是在高频段。

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值