Сложность построения спектрограммы данных ЭЭГ в Python

Я стремился создать желаемый график, используя аналогичные данные. Я начал с визуализации исходных и отфильтрованных сигналов, прежде чем создавать спектрограммы.

Я думаю, что ошибка, которую вы видели, заключалась в том, что вы назначали исходному кадру данных сигнал с пониженной дискретизацией, который не будет работать, поскольку они имеют разную длину. Чтобы обойти эту проблему, я присвоил уменьшенные данные новому фрейму данных 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)')

Это функция, которую я создал для повторной выборки сигналов ЭЭГ, поскольку у меня были проблемы с повторной выборкой Scipy.

def resample_filter(index, data, rate, target, interpolation='linear'):
    if target > rate:
        updata = []
        size = len(index)
        time = size // rate
        upindex = np.linspace(0, time, int(time * target))
        for i in range(len(data)):
            updata.append(scipy.interpolate.interp1d(np.linspace(0, time, size), data[i], kind=interpolation)(upindex))
        return upindex, updata, target
    factor = int(round(rate / target))
    return index[::1 if factor == 0 else factor], [data[i][::1 if factor == 0 else factor] for i in range(len(data))], rate // factor

Он возвращает новый индекс, новый сигнал и новую частоту дискретизации.

Редактировать :

Версия ресэмплера scipy.resample

def resample_filter(index, data, rate, target):
    time = len(index) // rate
    size = int(time * target)
    return np.linspace(0, time, size), [scipy.signal.resample(data[i], size) for i in range(len(data))], target

И функция, которую я использую для спектрограммы с окном в секундах, но она такая же, как у вас:

def stft_spectrogram(index, data, rate, window=1.0, type='blackman'):
    size = len(data)
    n = int(window * rate)
    f, t, sxx = signal.spectrogram(
        x=np.array(data), mode = "complex", window=type, fs=rate, nperseg=n if n <= size else size)
    return np.linspace(index[0], index[-1], len(sxx[0])), f, sxx

Он возвращает индекс времени, частоты и комплексы stft, где np.abs(sxx) — амплитуда частоты, а np.angle(sxx) — фаза.

Спектрограмма STFT позволяет увидеть, что происходит, но вы столкнетесь с принципом неопределенности между временной и частотной областями, и вам придется поиграть с размером окна, чтобы выбрать, на какой стороне вы предпочитаете видеть явно теряющие данные на другой стороне. Чтобы убедиться, что ваши фильтры эффективны, вам следует начать с построения простого графика спектральной плотности мощности, используя signal.welch. Это что-то вроде STFT. Вот пример:

STFT-спектрограмма сигнала ЭЭГ с европейским пиком переменного тока частотой 50 Гц из моей работы :(

Уэлч PSD сигнала

После простого Notch-фильтра:

После проверки эффективности ваших фильтров вы можете использовать спектрограмму STFT, но вам также следует подумать о построении некоторых вейвлетов Морлета или Гильберта Хуанга.