在Python中绘制脑电图数据的频谱图很困难

Question

我正在做一些任务，模仿以下链接的一些情节：

“https://www.nature.com/articles/s41598-023-37829-z”，

但是当我尝试使用 scipy 库在 python 中绘制脑电图频谱图时遇到问题。

我正在尝试绘制从 Muse 头带数据集记录的脑电图数据的频谱图。以下是我根据提供的说明执行的步骤：

1. 预处理：我使用 signal.butter 对原始 EEG 信号应用了 1 Hz 以上的高通滤波器和 50 Hz 以下的低通滤波器。然后，我尝试使用 signal.resample 和 signal.decimate 将信号从 256 Hz 下采样到 128 Hz，但遇到了 ValueError：“值的长度与索引的长度不匹配。”

2. 绘制频谱图：我使用 signal.spectrogram 来计算通道 TP9 和 TP10 的频谱图。然而，绘图时，只有轴可见，并且没有脑电图信号。

    # Extract the desired channels (TP9 and TP10)

    channels = ['TP9', 'TP10']  
    eeg_data = eeg_df[channels]
    
    # Convert the timestamps to time in seconds

    sampling_rate = 256  # Sampling rate of the EEG data in Hz
    time_seconds = eeg_df['timestamps'] / sampling_rate
    
    # Preprocess the EEG signals

    for channel in channels:
        # Apply high-pass filter above 1 Hz and low-pass filter below 50 Hz
        b, a = signal.butter(4, [1, 50], btype='bandpass', fs=sampling_rate)
        eeg_data[channel] = signal.filtfilt(b, a, eeg_data[channel])
        
        # Downsample the signal from 256 to 128 Hz

        new_length = len(eeg_data[channel]) // 2
        eeg_data[channel] = signal.resample(eeg_data[channel], new_length)
    
    # Define the parameters for the spectrogram

    window = 'hann'  # Windowing function

    nperseg = 128    # Length of each segment

    noverlap = 64    # Overlap between segments

    nfft = 128       # Number of points for the FFT
    
    # Compute the spectrogram for each channel

    frequencies, times, spectrogram_TP9 = signal.spectrogram(eeg_data['TP9'], fs=sampling_rate / 2, window=window,nperseg=nperseg, noverlap=noverlap, nfft=nfft)
    frequencies, times, spectrogram_TP10 = signal.spectrogram(eeg_data['TP10'], fs=sampling_rate / 2, window=window,nperseg=nperseg, noverlap=noverlap, nfft=nfft)
    
    # Plot the spectrogram for TP9

    plt.figure(figsize=(12, 6))
    plt.subplot(2, 1, 1)
    plt.pcolormesh(times, frequencies, 10 * np.log10(spectrogram_TP9), shading='gouraud', cmap='viridis')
    plt.title('Spectrogram for TP9')
    plt.xlabel('Time (seconds)')
    plt.ylabel('Frequency (Hz)')
    plt.colorbar(label='Power/Frequency (dB/Hz)')
    plt.ylim(0, 30)  # Limit the y-axis to frequencies up to 30 Hz
    
    # Plot the spectrogram for TP10

    plt.subplot(2, 1, 2)
    plt.pcolormesh(times, frequencies, 10 * np.log10(spectrogram_TP10), shading='gouraud', cmap='viridis')
    plt.title('Spectrogram for TP10')
    plt.xlabel('Time (seconds)')
    plt.ylabel('Frequency (Hz)')
    plt.colorbar(label='Power/Frequency (dB/Hz)')
    plt.ylim(0, 30)  # Limit the y-axis to frequencies up to 30 Hz
    
    plt.tight_layout()
    plt.show()

   

我怀疑该问题可能与下采样后的长度不匹配或信号处理函数的错误使用有关。有人可以帮我解决这个问题，并提供有关如何正确绘制脑电图数据频谱图的指导吗？我在图像中附上了我的错误，以及我的数据集信息和值。

如有任何帮助或建议，我们将不胜感激！谢谢！

相关图片：

Answer 1

我试图使用类似的数据创建所需的图。在创建频谱图之前，我首先可视化原始信号和滤波信号。

我认为您看到的错误是因为您将下采样信号分配给原始数据帧，这将不起作用，因为它们的长度不同。为了解决这个问题，我将下采样的数据分配给一个新的数据框

eeg_downsampled

。

原始和加工渠道：

频谱图：

---

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt

from scipy import signal

#
# Load some test data and modify it to look more like the OP's data
#
eeg_df = pd.read_csv('../musedata_meditation.csv')
eeg_df = df.rename(columns={'RAW_TP9': 'TP9', 'RAW_TP10': 'TP10'})
eeg_df['sample_number'] = range(len(eeg_df))

# Sampling rate of the EEG data in Hz
sampling_rate = 256  

# Convert the sample index to time in seconds
# This means you need to divide the sample *number* by the sampling rate
eeg_df['time_seconds'] = eeg_df['sample_number'] / sampling_rate

# Extract the desired channels (TP9 and TP10)
keep_channels = ['TP9', 'TP10']
eeg_data = eeg_df[keep_channels + ['time_seconds']].copy()

#
# Preprocess the EEG signals
#

#Define the pre-preprocessing parameters
downsample_factor = 2
bp_order = 4
bp_passband = [1, 50]
downsample_factor = 2

eeg_downsampled = pd.DataFrame()
downsampled_length = len(eeg_data) // downsample_factor

for channel in keep_channels:
    filtered_name = channel + '_bp'
    # Apply high-pass filter above 1 Hz and low-pass filter below 50 Hz
    b, a = signal.butter(bp_order, bp_passband, btype='bandpass', fs=sampling_rate)
    eeg_data[filtered_name] = signal.filtfilt(b, a, eeg_data[channel])
    
    # Downsample the signal from 256 to 128 Hz
    # Put the shorter signals in a new dataframe "eeg_downsampled"
    eeg_downsampled[channel] = signal.resample(eeg_data[filtered_name], downsampled_length)
eeg_downsampled['time_seconds'] = np.arange(new_length) / (sampling_rate / downsample_factor)

#
# View the original and processed data
#
plot_height = 2
n_plots = 3
f, axs = plt.subplots(nrows=n_plots, ncols=1,
                      figsize=(11, plot_height * n_plots),
                      layout='constrained', sharex=True)

plot_slice = slice(0, 2500) #window of data to plot

#Original data
eeg_data[plot_slice].plot(
    y=keep_channels, x='time_seconds',
    ylabel='signal', xlabel='time (s)',
    linewidth=1, ax=axs[0],
    title=f'original data (samples {plot_slice.start} to {plot_slice.stop})'
)

#Filtered
eeg_data[plot_slice].plot(
    y=[chan + '_bp' for chan in keep_channels], x='time_seconds',
    ylabel='signal', xlabel='time (s)', legend=False,
    linewidth=1, ax=axs[1], title=f'bandpass {bp_passband}Hz'
)

#Downsampled
eeg_downsampled[plot_slice.start:plot_slice.stop//downsample_factor].plot(
    y=keep_channels, x='time_seconds',
    ylabel='signal', xlabel='time (s)', legend=False,
    linewidth=1, ax=axs[2], title=f'downsampled {downsample_factor}x'
)


# Define the parameters for the spectrogram
window = 'hann'  # Windowing function
nperseg = 128    # Length of each segment
noverlap = 64    # Overlap between segments
nfft = 128       # Number of points for the FFT

# Compute the spectrogram for each channel
spectrograms = {}
for channel in keep_channels:
    frequencies, times, spectrogram = signal.spectrogram(
        eeg_downsampled[channel],
        fs=sampling_rate / downsample_factor,
        window=window, nperseg=nperseg, noverlap=noverlap, nfft=nfft
    )
    spectrograms[channel] = spectrogram

# Plot the spectrograms
f, axs = plt.subplots(nrows=2, ncols=1, figsize=(11, 5),
                      layout='constrained', sharex=True)
f.suptitle('Spectrograms')
vmin, vmax = -10, 30

for channel, ax in zip(keep_channels, axs):
    im = ax.pcolormesh(
        times, frequencies, 10 * np.log10(spectrograms[channel]),
        vmin=vmin, vmax=vmax, #comment out for auto-range
        shading='nearest', cmap='viridis'
    )
    
    ax.set_ylim(0, 30)  # Limit the y-axis to frequencies up to 30 Hz
    ax.set_ylabel(('Left' if channel=='TP9' else 'Right') + f' ({channel})\nfrequency (Hz)')
    
    cbar = f.colorbar(im, aspect=8, pad=0.01, label='power/frequency (dB/Hz)')
    
ax.set_xlabel('time (s)')