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

本文链接：https://blog.csdn.net/qq_46234417/article/details/146019039

代码来源及详情下载网址：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 方法的方差进一步增加。

关键观察

周期图法 vs. Welch 方法：
- 周期图法的频谱估计结果波动较大，尤其是在高频段，这表明其方差较高。
- Welch 方法通过分段和平均操作，显著降低了频谱估计的方差，结果更加平滑和稳定。
Welch 方法的参数影响：
- 窗口长度（M）：
  - 较长的窗口长度（如 M=160）可以提供更好的频率分辨率，但计算量较大。
  - 较短的窗口长度（如 M=80）会降低频率分辨率，但计算量较小。
- 重叠长度（D）：
  - 较大的重叠长度（如 D=80）可以进一步降低频谱估计的方差，结果更加平滑。
  - 无重叠（D=0）时，Welch 方法的方差较高，尤其是在高频段。
频率范围：
- 在低频段（0-10Hz），三种 Welch 方法的结果较为接近，且与周期图法的结果也较为接近。
- 在中频段（10-20Hz），Welch 方法的结果开始显示出平滑的优势。
- 在高频段（20-50Hz），Welch 方法的结果显著优于周期图法，尤其是当窗口长度较长且重叠长度较大时。