研一脑电小白日记之频谱与时频分析代码学习4 连续小波变换（CWT）时频分析

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

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

demo_mwt.m 连续小波变换（CWT）时频分析：使用Morlet小波分析VEP信号的时频特征，需调用子函数subfunc_mwt.m

subfunc_mwt.m CWT计算子函数：基于Morlet小波基实现时频分解，支持调节小波中心频率与带宽参数。

代码解释 demo_mwt.m

这段代码用于对视觉诱发电位（VEP）信号进行多尺度小波变换（MWT）的时间频率分析，并与短时傅里叶变换（STFT）的结果进行比较。以下是代码的逐步解释：

1. 数据预处理

N = numel(x); % 时间点的数量
x = x - mean(x(t<0)); % 基线校正

N：信号的时间点数量。
x = x - mean(x(t<0))：通过减去刺激前信号的平均值进行基线校正，以消除基线漂移。

2. 短时傅里叶变换（STFT）

nfft = 2^nextpow2(N); % FFT 点数
[P_stft,f] = subfunc_stft(x,100, nfft, Fs); % 计算 STFT，窗口大小为 100
N_F = length(f); % 频率 bin 的数量

nfft：计算大于等于信号长度的最小 2 的幂次，作为 FFT 点数。
subfunc_stft：调用自定义函数计算 STFT，窗口大小为 100。
P_stft：STFT 的功率谱密度。
f：频率轴。

3. 多尺度小波变换（MWT）

omega = [0.5 1];    % 测试两个 omega 值
sigma = [0.5 .2];   % 测试两个 sigma 值

ff = f/Fs;     % 归一化频率
L_hw = N; % 滤波器长度
for fi=1:N_F
    scaling_factor = omega(1)./ff(fi);
    u = (-[-L_hw:L_hw])./scaling_factor;
    w(:,fi) = sqrt(1/scaling_factor)*exp(-(u.^2)/(2*sigma(1).^2));
    hw(:,fi) = w(:,fi).'.* exp(1i*2*pi*omega(1)*u);
end

omega：测试两个尺度参数。
sigma：测试两个形状参数。
ff：归一化频率。
L_hw：滤波器长度。
scaling_factor：尺度因子。
u：时间索引。
w：小波核的权重。
hw：小波核的复数表示。

for n_omega=1:numel(omega)
    for n_sigma=1:numel(sigma)
        [P_mwt(:,:,n_omega,n_sigma)] = subfunc_mwt(x, f, Fs, omega(n_omega), sigma(n_sigma));
    end
end

subfunc_mwt：调用自定义函数计算 MWT，返回功率谱密度 P_mwt。

4. 显示结果

f_lim = [min(f(f>0)), 30]; % 频率范围（去除 0Hz）
f_idx = find((f<=f_lim(2))&(f>=f_lim(1)));
t_lim = [-0.2, 1]; % 时间范围
t_idx = find((t<=t_lim(2))&(t>=t_lim(1)));

figure('units','normalized','position',[0.1    0.15    0.8    0.7])
for n_omega=1:numel(omega)
    for n_sigma=1:numel(sigma)
        n = (n_omega-1)*2+n_sigma;
        subplot(2,2,n)
        imagesc(t(t_idx),f(f_idx),P_mwt(f_idx,t_idx,n_omega,n_sigma))
        xlabel('Time (s)')
        ylabel('Frequency (Hz)')
        set(gca,'xlim',t_lim,'ylim',f_lim)
        axis xy; hold on;
        plot([0 0],[0 Fs/2],'w--')
        title(['Scaleogram (\omega = ',num2str(omega(n_omega),'%1.2g'),', \sigma = ',num2str(sigma(n_sigma),'%1.2g'),')'],'fontsize',12)
        text(t_lim(2),f_lim(2)/2,'Power (\muV^2/Hz)','rotation',90,'horizontalalignment','center','verticalalignment','top')
        colorbar 
    end
end

f_lim：限定频率范围为 0Hz 到 30Hz。
t_lim：限定时间范围为 -0.2 秒到 1 秒。
imagesc：绘制 MWT 的结果（尺度图），显示不同时间和频率下的功率。
xlabel 和 ylabel：设置坐标轴标签。
set(gca,'xlim',t_lim,'ylim',f_lim)：设置时间和频率范围。
axis xy：设置坐标轴方向。
plot([0 0],[0 Fs/2],'w--')：绘制刺激点的参考线。
title：设置子图标题，显示参数值。
colorbar：添加颜色条，表示功率（dB）。

总结

这段代码通过 MWT 方法对 VEP 信号进行时间频率分析，并与 STFT 的结果进行比较。主要步骤包括：

数据预处理：基线校正。
计算 STFT：用于对比参考。
计算 MWT：通过不同的尺度参数和形状参数，得到不同时间频率分布。
显示结果：绘制不同参数下的尺度图，展示信号在不同时间和频率下的功率分布。

代码解释 subfunc_mwt.m

这段代码实现了一个基于 Morlet 小波变换（MWT）的时间频率分析，用于分析信号在不同时间和频率下的功率分布。以下是代码的逐步解释：

1. 预处理和参数设置

if size(x,2)==1; x = x.'; end
x = detrend(x,'linear'); % 去除线性趋势

N_F = length(f); % 频率 bin 的数量
N_T = length(x); % 时间样本数量
f = f/Fs;     % 归一化频率

S = single(zeros(N_F,N_T)); % 定义输出的大小

x：输入信号，假设是一个一维数组。
detrend：去除信号中的线性趋势，以消除基线漂移。
N_F：频率 bin 的数量。
N_T：时间样本数量。
f：归一化频率。
S：初始化一个二维数组，用于存储 MWT 的复数结果。

2. Morlet 小波变换

L_hw = N_T; % 滤波器长度
for fi=1:N_F
    scaling_factor = omega./f(fi); % 尺度因子
    u = (-[-L_hw:L_hw])./scaling_factor;
    hw = sqrt(1/scaling_factor)*exp(-(u.^2)/(2*sigma.^2)).* exp(1i*2*pi*omega*u);
    S_full = conv(x,conj(hw));
    S(fi,:) = S_full(L_hw+1:L_hw+N_T); % 复数值
end
P = abs(S).^2;  % 功率值

L_hw：滤波器长度，设置为时间样本数量。
scaling_factor：尺度因子，用于调整小波的尺度。
u：时间索引，用于构造小波核。
hw：Morlet 小波核，由高斯窗和复数正弦波组成。
- sqrt(1/scaling_factor)*exp(-(u.^2)/(2*sigma.^2))：高斯窗部分。
- exp(1i*2*pi*omega*u)：复数正弦波部分。
conv：对信号 x 和小波核 hw 进行卷积，得到复数结果 S_full。
S(fi,:)：提取卷积结果中对应的时间段，存储在 S 中。
P：计算功率值，P = abs(S).^2，单位为 μV2/Hz。

3. 参数说明

omega：Morlet 小波的中心频率。
sigma：Morlet 小波的带宽参数，控制小波的时间和频率分辨率。
f：频率轴，归一化到采样率 Fs。

4. 总结

这段代码通过 Morlet 小波变换（MWT）对信号进行时间频率分析，主要步骤包括：

预处理：去除信号中的线性趋势。
参数设置：定义频率 bin 的数量、时间样本数量和归一化频率。
小波核构造：构造 Morlet 小波核，由高斯窗和复数正弦波组成。
卷积计算：对信号和小波核进行卷积，得到复数结果。
功率计算：计算功率值，用于显示信号在不同时间和频率下的功率分布。

通过这些步骤，可以得到信号的时间频率表示（TFR），用于分析信号在不同时间和频率下的特性。

运行结果

结果分析

这段代码的运行结果图展示了使用 Morlet 小波变换（MWT）对视觉诱发电位（VEP）信号进行时间频率分析的结果。图中显示了不同参数设置下的尺度图（Scaleogram）。以下是对运行结果图的详细分析：

1. 尺度图（Scaleogram）

横轴：时间（Time，单位为秒，s）。
纵轴：频率（Frequency，单位为赫兹，Hz）。
颜色条：功率（Power，单位为 μV2/Hz），颜色越亮表示功率越高。

2. 参数设置

omega：Morlet 小波的中心频率，测试了两个值：0.5 和 1。
sigma：Morlet 小波的带宽参数，测试了两个值：0.5 和 0.2。

3. 结果分析

3.1 尺度图（ω=0.5,σ=0.5）

频率范围：0-30Hz。
时间范围：-0.2 秒到 1 秒。
主要特征：
- 在时间 0 之后，频率成分主要集中在 5Hz 左右。
- 功率较高的频率成分出现在 5Hz 左右，表明视觉刺激引发了特定频率成分的变化。
- 颜色条的范围较宽，表明功率变化较大。

3.2 尺度图（ω=0.5,σ=0.2）

频率范围：0-30Hz。
时间范围：-0.2 秒到 1 秒。
主要特征：
- 在时间 0 之后，频率成分主要集中在 5Hz 左右。
- 功率较高的频率成分仍然集中在 5Hz 左右，但颜色条的范围较窄，表明功率变化较小。
- 由于 σ 较小，小波的带宽较窄，频率分辨率较高。

3.3 尺度图（ω=1,σ=0.5）

频率范围：0-30Hz。
时间范围：-0.2 秒到 1 秒。
主要特征：
- 在时间 0 之后，频率成分主要集中在 5Hz 左右。
- 功率较高的频率成分仍然集中在 5Hz 左右，但颜色条的范围较宽，表明功率变化较大。
- 由于 ω 较大，小波的中心频率较高，时间分辨率较高。

3.4 尺度图（ω=1,σ=0.2）

频率范围：0-30Hz。
时间范围：-0.2 秒到 1 秒。
主要特征：
- 在时间 0 之后，频率成分主要集中在 5Hz 左右。
- 功率较高的频率成分仍然集中在 5Hz 左右，但颜色条的范围较窄，表明功率变化较小。
- 由于 ω 较大，σ 较小，小波的中心频率较高且带宽较窄，时间分辨率和频率分辨率都较高。

4. 关键观察

频率成分：
- 在时间 0 之后，频率成分主要集中在 5Hz 左右，表明视觉刺激引发了特定频率成分的变化。
- 无论参数如何设置，这一特征都较为明显。
参数影响：
- omega：较大的 ω 值提高了时间分辨率，但可能会降低频率分辨率。
- sigma：较小的 σ 值提高了频率分辨率，但可能会降低时间分辨率。
- 选择合适的参数需要根据具体应用需求进行权衡。
功率分布：
- 功率较高的频率成分主要集中在 5Hz 左右，表明这一频率成分在 VEP 信号中占主导地位。
- 不同参数设置下的功率分布略有差异，但总体趋势一致。