脑电信号时频分析仿真

zxdnjp

已于 2025-10-29 16:23:16 修改

阅读量140

点赞数 5

CC 4.0 BY-SA版权

分类专栏：脑电信号处理学习笔记文章标签： matlab 信号处理傅里叶分析

于 2025-10-29 16:22:12 首次发布

本文链接：https://blog.youkuaiyun.com/zxdnjp/article/details/154073197

脑电信号处理学习笔记专栏收录该内容

5 篇文章

订阅专栏

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

文章目录

前言
一、时频分析是什么？
二、代码脚本
总结

前言

提示：仅作个人学习记录用：

EEG（脑电信号）处理过程中，肯定避免不了对采集到的信号进行时频分析并绘制频谱图，所以本片笔记就记录一下博主自己写的一个简单脚本，用于绘制信号时频分析结果。

提示：以下是本篇文章正文内容，下面案例可供参考

一、时频分析是什么？

脑电信号的时频分析是观察大脑电活动随时间变化的频率成分。把脑电波拆解成不同频段的成分（如α波、β波），然后看这些成分如何随着时间变化。它能同时展示 “什么时候” 发生了 “什么样的” 脑活动，比单独看波形或频谱更能捕捉大脑动态的变化过程。这种方法可以揭示认知任务中或疾病状态下脑电活动的细微波动，是研究大脑功能的重要工具。

话不多说！！你能看到这篇帖子肯定是已经有基础认知了！！上代码！！！！

二、代码脚本

代码如下（示例）：

clear;
close all;
fclose('all');
clc;

%% 信号频谱分析 - FFT实现

% 参数设置
Fs = 128;               % 采样频率 128 Hz
total_points = 120000;  % 总点数 120000
window_duration = 2;    % 窗长 2秒
step_duration = 0.5;    % 步长 0.5秒

% 计算窗口参数
window_length = Fs * window_duration;      % 窗长点数: 256点
step_length = Fs * step_duration;          % 步长点数: 64点
overlap_length = window_length - step_length; % 重叠点数: 192点

% 生成示例信号
% 假设信号包含10Hz和25Hz成分，加上一些噪声
t = (0:total_points-1)/Fs;
signal = 2*sin(2*pi*10*t) + 1.5*sin(2*pi*25*t) + 0.5*randn(size(t));

num_windows = floor((total_points - overlap_length) / step_length);
freq_axis = (0:window_length/2) * Fs / window_length;  % 单边频谱，y轴的索引

% 预分配存储时频结果的矩阵
time_freq_matrix = zeros(length(freq_axis), num_windows);
time_axis = zeros(1, num_windows);

% 创建汉明窗
window_func = hamming(window_length);

% 执行滑动窗口FFT（其实就是每次取一段窗口进行功率谱分析，然后画出来对应的功率谱图）
for i = 1:num_windows
    % 计算当前窗口的起始和结束索引
    start_index = (i-1) * step_length + 1;
    end_index = start_index + window_length - 1;
    
    % 提取当前窗口信号并加窗
    if end_index <= total_points
        windowed_signal = signal(start_index:end_index) .* window_func';
        
        % 执行FFT
        NFFT = 2^nextpow2(window_length);  % 使用2的幂次方提高FFT效率
        Y = fft(windowed_signal, NFFT);
        
        % 计算单边幅度谱 (取前NFFT/2+1个点)
        P2 = abs(Y/NFFT);       % 能量归一化，使频谱幅度反映信号真实幅度，也是为了补偿FFT求和操作导致的幅度放大
        P1 = P2(1:NFFT/2+1);    % 实信号的频谱是共轭对称的，只需要一半，还能减少数据量，提高后续处理效率
        P1(2:end-1) = 2*P1(2:end-1);    % 能量补偿，因为单边谱只包含了一半的能量信息，这样可以使单边谱的幅度与原始信号幅度对应
        
        % 确保P1长度与freq_axis匹配
        if length(P1) > length(freq_axis)
            P1 = P1(1:length(freq_axis));
        elseif length(P1) < length(freq_axis)
            P1(end+1:length(freq_axis)) = 0;
        end
        
        % 存储结果
        time_freq_matrix(:, i) = P1;        % 每列是一组结果
        time_axis(i) = (start_index + window_length/2) / Fs;  % 窗口中心时间
    end
end

%% 绘制时频分析结果
figure('Position', [100, 100, 1200, 800]);  % 控制画布位置和画布大小

% 子图1: 原始信号
subplot(3,1,1);
plot(t, signal);
title('原始信号 (120000点, 128Hz采样)');
xlabel('时间 (s)');
ylabel('幅度');
grid on;

% 子图2: 时频图 (手动FFT结果)
subplot(3,1,2);
imagesc(time_axis, freq_axis, 20*log10(time_freq_matrix + eps));  % 转换为dB
axis xy;  % 确保频率轴方向正确
colorbar;
title('时频分析 - 手动FFT实现 (窗长2s, 步长0.5s)');
xlabel('时间 (s)');
ylabel('频率 (Hz)');
ylim([0, Fs/2]);  % 显示奈奎斯特频率以下

% 子图3: 使用spectrogram函数对比
subplot(3,1,3);
[S, F, T] = spectrogram(signal, hamming(window_length), overlap_length, 1024, Fs);
surf(T, F, 10*log10(abs(S)), 'EdgeColor', 'none');
axis tight;
view(0, 90);
colorbar;
title('时频分析 - spectrogram函数结果');
xlabel('时间 (s)');
ylabel('频率 (Hz)');