代码来源及详情下载网址: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 方法的结果显著优于周期图法,尤其是当窗口长度较长且重叠长度较大时。
-
总结
-
Welch 方法的优势:
-
通过分段和平均操作,显著降低了频谱估计的方差,结果更加平滑和稳定。
-
适用于需要高精度频谱估计的场景,尤其是在高频段。
-
-
参数选择的影响:
-
较长的窗口长度和较大的重叠长度可以进一步降低方差,但会增加计算量。
-
无重叠时,Welch 方法的方差较高,尤其是在高频段。
-